Investigation of Air Pocket Behavior in Pipelines Using Rigid Column Model and Contributions of Time Integration Schemes

Abstract

This paper studies the air pressurization problem caused by a partially pressurized transient flow in a reservoir-pipe system. The purpose of this study is to analyze the performance of the rigid column model in predicting the attenuation of the air pressure distribution. In this regard, an analytic formula for the amplitude and frequency will be derived, in which the influential parameters, particularly, the driving pressure and the air and water lengths, on the damping can be seen. The direct effect of the driving pressure and inverse effect of the product of the air and water lengths on the damping will be numerically examined. In addition, these numerical observations will be examined by solving different test cases and by comparing to available experimental data to show that the rigid column model is able to predict the damping. However, due to simplified assumptions associated with the rigid column model, the energy dissipation, as well as the damping, is underestimated. In this regard, using the backward Euler implicit time integration scheme, instead of the classical fourth order explicit Runge–Kutta scheme, will be proposed so that the numerical dissipation of the backward Euler implicit scheme represents the physical dissipation. In addition, a formula will be derived to calculate the appropriate time step size, by which the dissipation of the heat transfer can be compensated.

1. Introduction

The operational problems in sewer network systems, either due to overflowing or blockage at pipeline ends, disturb the boundary conditions and result in the transition of the gravity flow to a partially pressurized flow, and then the high pressurized air pocket could be released to the atmosphere either at upstream or downstream boundaries. For numerical simulation, among various models, the rigid column (RC) model is less complicated, while incorporates the fundamental features of this type of transient flows. The RC model has been studied in extensive research exemplified by the following studies [1,2,3,4,5,6,7,8,9,10]: Cabrera et al. [2], Liou and Hunt [3], Li and McCorquodale [4], Zhou et al. [5,6], Lee [7], Fuertes-Miquel et al. [8], Coronado-Hernández et al. [9,10].

Li and McCorquodale [4] used the RC model proposed in [1] to calculate the transient pressure which causes blowing off of storm manholes. In the proposed RC model, the water column is assumed as an incompressible fluid with uniform but unsteady velocity. The transition from free-surface flow to pressurized flow was described in six stages, and the trapped air pocket was assumed to undergo pseudo adiabatic expansion and compression. It was found that the numerical transient pressure overpredicts the experimental data by a factor of two, and the attenuation of the pressure oscillation is underestimated, while the frequency of pressure surge oscillations predicts the measured data accurately. Li and McCorquodale [4] linked this poor behavior to the superposition of various air release processes and the steady-state friction factor used in the mathematical model.

Liou and Hunt [3] proposed an RC model to simulate the transient flow in an empty pipeline with an undulating elevation profile. They found that the calculated velocity distribution is comparable with the experimental data.

Lee [7] analytically and experimentally studied the behavior of air pressurization in horizontal pipelines and derived a closed form of the solutions calculated by the RC model, in which a vertical air-water interface is assumed. It was found that the maximum pressure increases when the reservoir pressure increases, and the RC model overpredicts the experimental data. In addition, Lee showed that the damping behavior of the air pressure is contributed to the heat transfer between the air, water, and the pipe wall. Additionally, he showed that the maximum pressure calculated by both variable and constant water length assumptions is similar, while he indicated that these two assumptions lead to calculating different times at which the maximum pressure occurs. Lee demonstrated that the maximum pressure in a frictionless pipe is independent of initial air and water lengths. Moreover, Lee, by using the small wave theory, linearized the governing equations and provided a formula for the frequency of the air pressure oscillation. This formula showed that the frequency depends on the reservoir pressure and the initial length of the air and water phases.

Hou et al. [11] experimentally and numerically studied the rapid filling of a large-scale pipeline. For the numerical calculation, they applied the RC model proposed in [3] along with the modification proposed by Axworthy and Karney [12]. Note that the modification is to eliminate the pressure head at junctions from the momentum equation. They showed that the numerical result of the RC model is in good agreement with the experimental data.

Fuertes-Miquel et al. [8], using the RC model, analyzed the effects of isothermal and adiabatic behaviors of the air pocket, which is expelled from a pipe through an air valve by admitting water flow. It was found that numerical results can be significantly different depending on the heat transfer mechanism. In particular, in contrary to the isothermal assumption, the adiabatic assumption results in smaller air pressure and larger water velocity so that the induced water hammer is greater.

Coronado-Hernández et al. [9] studied the sub-atmospheric pressure occurring during an emptying process in a pipeline with irregular elevation. They proposed an RC model to simulate the water flow, while the air pocket was simulated using the thermodynamic formulation. By comparing to experimental data, they showed that the RC model can predict the flow variables (water velocity and air pressure) accurately.

The air entrapment following rapidly filling storm water systems (SWSs) has great complexity and challenges. This complexity, for example, includes the transition from gravity to pressurized flow and their interaction and significant head loss at the interface, the turbulent nature of the flow, lumped nature of the air pocket [5], and the heat transfer between water, air, and the pipe wall [7]. It is known that the available mathematical models, which practitioners are interested in, including the RC model, do not take the effects of all of these parameters into account so that these models overestimate the numerical solutions including the maximum pressure. This overestimation implies that these mathematical models, due to their assumptions, are not able to sufficiently predict the dissipation of the driving energy. Thus, to dissipate this energy, or in other words, to compensate for the effects of the aforementioned neglected terms, considering a numerical dissipation term seems to be useful. Note that in the literature, the classical fourth order Runge–Kutta explicit scheme has mostly been employed for the time integration of the governing equations. Recently, the backward Euler implicit scheme has been used in [13] to simulate the dynamics of the entrapped air pocket using the 3D Navier–Stokes equation and the volume of fluid (VOF) method. In addition, Rokhzadi and Fuamba [14] proposed using the backward Euler implicit scheme to solve a similar problem by a shock-fitting approach.

In this paper, the damping and frequency of the pressure distribution of an entrapped air caused by a rapid filling of a horizontal reservoir-pipe system will be studied using some available experimental data as well as the numerical results of the RC model. In this regard, the effects of reservoir pressure and the air and water initial lengths will be discussed. To analytically show the relation between these parameters and the damping and frequency, the governing equations are linearized around a point, in which the air pressure is equal to the reservoir pressure. This linearization implies that the air pressure does not significantly deviate from the reservoir pressure. Although this assumption is not valid for every case, it can provide insightful information about the damping and frequency. In addition, using the backward Euler implicit scheme will be introduced to discuss the numerical dissipation as a representative for the physical dissipation, which the RC model cannot predict properly for many reasons. In this regard, the modified equation will be derived to provide a formula, by which the size of the time step and the amount of numerical dissipation can be controlled so that excessive numerical dissipation imposed on the solutions is avoided. It will be shown that this formula adjusts the time step in proportion to the values of the effective physical parameters on the damping behavior.

It is worth mentioning that for problems with large initial air pocket sizes, the elasticity of entrapped air pocket is much higher than the elasticity of the water [15]. Thus, it can be expected that for these cases, which are more problematic in SWSs, the RC model can calculate the solutions appropriately. However, as indicated in the literature, e.g., [16], the RC model cannot capture the interactions between different waves exist in pipelines.

2. Numerical Analysis

2.1. Governing Equation

In this paper, the air entrapment caused by transient flow in a horizontal reservoir-pipe system is studied. As shown in Figure 1, a pipe segment, with the diameter $\mathrm{D}$ and length ${\mathrm{L}}_{\mathrm{t}}$, is connected at upstream to a reservoir, in which the absolute pressure is ${\mathrm{p}}_{\mathrm{R}}$, and is blocked at downstream by a valve. The water in a stationary pressurized flow regime occupies an upstream part of the pipe with the initial length $\mathrm{L}$. The entrapped air occupies the rest of the pipe at downstream, which is separated by a valve from the water column, with the initial length ${\mathrm{L}}_{\mathrm{a}}$, which ${\mathrm{L}}_{\mathrm{a}}={\mathrm{L}}_{\mathrm{t}}-\mathrm{L}$. Note that the interface of air–water phases is assumed vertical so that the gravity flow forming during the transient flow is neglected.

The momentum equation for a horizontal water column is,

$$\frac{\mathrm{dV}}{\mathrm{dt}}=\frac{\left({\mathrm{p}}_{\mathrm{R}}-\mathrm{p}\right)}{\mathsf{\rho }\mathrm{L}}-\frac{\mathrm{f}}{2\mathrm{D}}\mathrm{V}\left|\mathrm{V}\right|,$$

(1)

where $\mathrm{t}$ is the time variable, $\mathsf{\rho }$ and $\mathrm{V}$ are the water density and velocity, respectively, and $\mathrm{g}$ is the gravitational constant ($\mathrm{g}=9.81{\mathrm{ms}}^{-2}$), $\mathrm{p}$ is the air absolute pressure, and $\mathrm{f}$ is the friction factor. Note that the length of the water column is assumed constant as Lee [7] demonstrated that the variable length does not provide significantly different results compared to the constant length assumption.

The continuity equation gives a relationship between the air volume rate and the water flow rate as,

$$\frac{\mathrm{d}\forall }{\mathrm{dt}}=-\mathrm{VA},$$

(2)

where $\mathrm{A}$ is the pipe cross-section area and $\forall$ is the air volume. Note that the velocity is assumed positive in the downstream direction. Following Lee [7], applying the first thermodynamic law to the air pocket and assuming the air as a perfect gas yield,

$$\frac{\mathrm{dp}\forall }{\mathrm{dt}}=-\left(\mathrm{k}-1\right)\mathrm{p}\frac{\mathrm{d}\forall }{\mathrm{dt}},$$

(3)

where $\mathrm{k}$ is the polytropic exponent. Note that Equation (3) implies a balance between the internal energy and the work done through the air compression and expansion. However, Lee [7] attributed the damping of the air pressure distribution to the heat transfer. Thus, to consider the effect of heat transfer, the term $\left(\mathrm{k}-1\right)\mathrm{q}$ must be added to the right hand side (RHS) of Equation (3). In addition, in the present study, the polytropic exponent is set to $\mathrm{k}=1.4$, as indicated in [7], meaning that this value gives more appropriate results. Substituting Equation (2) into Equation (3) yields,

$$\frac{\mathrm{dp}\forall }{\mathrm{dt}}=\left(\mathrm{k}-1\right)\mathrm{pVA}.$$

(4)

2.2. Linearized Governing Equations

As indicated in [7], the variable water column length does not have significant effects on the solutions. Thus, the only nonlinearities in the governing equations are related to the friction force in the momentum equation and the polytropic process equation $\left(\mathrm{p}{⩝}^{\mathrm{k}}=\mathrm{cte}\right)$. As shown in [7], for small ranges of reservoir pressure $({\mathrm{p}}_{\mathrm{R}}<2.0{\mathrm{p}}_{\mathrm{a}})$, in which ${\mathrm{p}}_{\mathrm{a}}$ is the atmospheric absolute pressure $\left({\mathrm{p}}_{\mathrm{a}}=101\mathrm{kPa}\right)$, the maximum air pressure, so as the minimum air pressure, slightly changes with the friction loss coefficient $\left({\mathrm{fL}}_{\mathrm{a}}/\mathrm{D}\right)$ and it is almost constant for both polytropic coefficients $(\mathrm{k}=1.2$ and $1.4)$. Thus, for simplicity of the discussion presented only in Section 2.2 and Section 2.3, the friction and local loss terms are neglected, meanwhile, the contributions of these terms to the damping behavior of the entrapped air pressure have been discussed in the literature.

Thus, the linearized governing equations are expected to appropriately approximate the solutions of the nonlinear equations when ${\mathrm{p}}_{\mathrm{R}}/{\mathrm{p}}_{\mathrm{a}}<2.0$. Therefore, the governing equations are linearized around a point, in which $\mathrm{p}={\mathrm{p}}_{\mathrm{R}}$, so that the effects of different variables on the damping and frequency of the air pressure distribution can be analyzed.

A continuous function of two variables $\mathrm{f}\left(\mathrm{u},\mathrm{v}\right)$, in which $\mathrm{u}$ and $\mathrm{v}$ denote the dependent variables, could be linearized using the Taylor series expansion around any point $\left({\mathrm{u}}_{\mathrm{R}},{\mathrm{v}}_{\mathrm{R}}\right)$ as,

$$\mathrm{f}\left(\mathrm{u},\mathrm{v}\right)=\mathrm{f}\left({\mathrm{u}}_{\mathrm{R}},{\mathrm{v}}_{\mathrm{R}}\right)+{\frac{\partial \mathrm{f}}{\partial \mathrm{u}}|}_{\left({\mathrm{u}}_{\mathrm{R}}\right)}\left(\mathrm{u}-{\mathrm{u}}_{\mathrm{R}}\right)+{\frac{\partial \mathrm{f}}{\partial \mathrm{v}}|}_{\left({\mathrm{v}}_{\mathrm{R}}\right)}\left(\mathrm{v}-{\mathrm{v}}_{\mathrm{R}}\right).$$

(5)

Note that the momentum equation, Equation (1), and the continuity equation, Equation (2), for a frictionless horizontal pipe are already linear. The energy conservation equation, Equation (4), can be linearized using Equation (5) as,

$$\frac{\mathrm{dp}}{\mathrm{dt}}=\mathrm{kA}\frac{{\mathrm{p}}_{\mathrm{R}}}{{\forall }_{\mathrm{R}}}\mathrm{V}+\left(\mathrm{k}-1\right)\frac{{\mathrm{AV}}_{\mathrm{R}}}{{\forall }_{\mathrm{R}}}\left(\mathrm{p}-{\mathrm{p}}_{\mathrm{R}}\right),$$

(6)

where the subscript $\left(\mathrm{R}\right)$ implies the value of a variable when the air pressure equals to the reservoir pressure.

2.3. Pressure Distribution Analysis

The linearized governing equations, Equations (1) and (6), can be rearranged in a vector form as,

$$\frac{\mathrm{d}\overrightarrow{\mathrm{F}}}{\mathrm{dt}}=G\overrightarrow{\mathrm{F}},$$

(7)

where,

$$\overrightarrow{\mathrm{F}}=\left\{\begin{array}{c}\mathrm{V}\\ {\mathrm{p}}^{*}\end{array}\right\},G=\left(\begin{array}{cc}0& -1/\mathsf{\rho }\mathrm{L}\\ {\mathrm{kAp}}_{\mathrm{R}}/{⩝}_{\mathrm{R}}& \left(\mathrm{k}-1\right){\mathrm{AV}}_{\mathrm{R}}/{⩝}_{\mathrm{R}}\end{array}\right),$$

(8)

where ${\mathrm{p}}^{*}=\mathrm{p}-{\mathrm{p}}_{\mathrm{R}}$. Two time-dependent decoupled equations can be derived using the Eigen-decomposition method, by which the numerical solutions, which have an exponential form in terms of the time, can be found. The amplitude and the frequency of the solutions are related to the eigenvalues of the matrix $\mathrm{G}$, which can be calculated as,

$${\mathsf{\lambda }}_{1,2}=\frac{1}{2}\left\{-\left(\mathrm{k}-1\right){\left(\frac{{\mathrm{p}}_{\mathrm{R}}}{{\mathrm{p}}_{\mathrm{a}}}\right)}^{1/\mathrm{k}}\frac{{\mathrm{V}}_{\mathrm{R}}}{{\mathrm{L}}_{\mathrm{a}}}\pm \mathrm{i}\sqrt{4\mathrm{k}\frac{1}{{\mathrm{L}}_{\mathrm{a}}\mathrm{L}}\frac{{\mathrm{p}}_{\mathrm{R}}}{\mathsf{\rho }}{\left(\frac{{\mathrm{p}}_{\mathrm{R}}}{{\mathrm{p}}_{\mathrm{a}}}\right)}^{1/\mathrm{k}}-{\left(\left(\mathrm{k}-1\right){\left(\frac{{\mathrm{p}}_{\mathrm{R}}}{{\mathrm{p}}_{\mathrm{a}}}\right)}^{1/\mathrm{k}}\frac{{\mathrm{V}}_{\mathrm{R}}}{{\mathrm{L}}_{\mathrm{a}}}\right)}^2}\right\}$$

(9)

where $\mathrm{i}$ denotes the complex unit number. Note that the polytropic process relationship, i.e., $\mathrm{p}{⩝}^{\mathrm{k}}=\mathrm{cte}$, and the definition of $⩝={\mathrm{AL}}_{\mathrm{a}}$ were used to further rearrange Equation (9). Thus, the amplitude and the frequency of the air pressure distribution can be calculated as,

$$\mathrm{M}={\mathrm{e}}^{-\frac{1}{2}\left(\left(\mathrm{k}-1\right){\left(\frac{{\mathrm{p}}_{\mathrm{R}}}{{\mathrm{p}}_{\mathrm{a}}}\right)}^{1/\mathrm{k}}\frac{{\mathrm{V}}_{\mathrm{R}}}{{\mathrm{L}}_{\mathrm{a}}}\right)},$$

(10)

$$\mathsf{\phi }=\sqrt{4\mathrm{k}\frac{1}{{\mathrm{L}}_{\mathrm{a}}\mathrm{L}}\frac{{\mathrm{p}}_{\mathrm{R}}}{\mathsf{\rho }}{\left(\frac{{\mathrm{p}}_{\mathrm{R}}}{{\mathrm{p}}_{\mathrm{a}}}\right)}^{1/\mathrm{k}}-{\left(\left(\mathrm{k}-1\right){\left(\frac{{\mathrm{p}}_{\mathrm{R}}}{{\mathrm{p}}_{\mathrm{a}}}\right)}^{1/\mathrm{k}}\frac{{\mathrm{V}}_{\mathrm{R}}}{{\mathrm{L}}_{\mathrm{a}}}\right)}^2},$$

(11)

where $\mathsf{\phi }$ is the frequency and $\mathrm{M}$ is the amplitude. Note that the damping behavior can be realized from the amplitude $\left(\mathrm{M}\right)$ in Equation (10). As explained, ${\mathrm{V}}_{\mathrm{R}}$ denotes the water column velocity at a time when the air pressure is equal to reservoir pressure. Malekpour et al. [17] analyzed the energy conversion of a transient flow in a reservoir-pipe system and derived an equation for the water column velocity in terms of initial air and water column lengths, reservoir pressure, and the air pressure. However, a simple equation for the cases with constant water column length assumption is derived here.

Multiplying both sides of Equation (1), in which the friction term is neglected, with $\mathrm{VA}$ and summing up with Equation (4) yields an energy equation as,

$$\frac{\mathrm{d}}{\mathrm{dt}}\left(\mathsf{\rho }\mathrm{LA}\frac{{\mathrm{V}}^2}{2}+\frac{\mathrm{p}⩝}{\mathrm{k}-1}\right)={\mathrm{p}}_{\mathrm{R}}\mathrm{AV}.$$

(12)

For an initially stationary water column, taking integral of both side of Equation (12) from initial condition to the time when the pressure is equal to the air pressure, so called ${\mathrm{t}}_{\mathrm{R}}$, yields,

$$\mathsf{\rho }\mathrm{LA}\frac{{\mathrm{V}}_{\mathrm{R}}^2}{2}+\frac{{\mathrm{p}}_{\mathrm{R}}{⩝}_{\mathrm{R}}}{\mathrm{k}-1}-\frac{{\mathrm{p}}_0{⩝}_0}{\mathrm{k}-1}={{\displaystyle \int }}_0^{{\mathrm{t}}_{\mathrm{R}}}{\mathrm{p}}_{\mathrm{R}}\mathrm{AVdt}.$$

(13)

Since the time taken by the air pressure to reach the reservoir pressure is not so long, for initially stationary flow, the integral on RHS of Equation (13) is negligible. Thus, the following relation for ${\mathrm{V}}_{\mathrm{R}}$ can be derived,

$${\mathrm{V}}_{\mathrm{R}}^2=\frac{2{\mathrm{p}}_{\mathrm{R}}{\mathrm{L}}_{\mathrm{a}}\left({\mathrm{p}}_0-{\left({\mathrm{p}}_0/{\mathrm{p}}_{\mathrm{R}}\right)}^{1/\mathrm{k}}\right)}{\mathsf{\rho }\mathrm{L}\left(\mathrm{k}-1\right)}.$$

(14)

As can be seen, ${\mathrm{V}}_{\mathrm{R}}$ is a function of $\sqrt{{\mathrm{p}}_{\mathrm{R}}{\mathrm{L}}_{\mathrm{a}}/\mathrm{L}}$. Note that a similar function can be found in [17] (Equation (20)). Therefore, Equations (10) and (11) show that the frequency and the damping of the air pressure distribution increase as the reservoir pressure increases and as the product of the air and water lengths $\left({\mathrm{L}}_{\mathrm{a}}\mathrm{L}\right)$ decreases. Later, these relationships will be confirmed by explaining the results of some available experiments. Therefore, Equation (10) shows that the RC model is able to predict the damping behavior of the air pressure distribution. However, in the next section, it will be shown that this ability depends on the time integration scheme.

2.4. Time Integration Scheme

From the literature, it can be realized that commonly the classical fourth order Runge–Kutta scheme, hereafter called the 4th-RK scheme, has been mostly used in the literature to integrate the governing equations of the RC model [6,11]. In addition, it is known that the numerical results are always overestimated, caused by underestimating the physical energy damping for many reasons including the assumption of lumped air pocket process [5], and neglecting the heat transfer between the air, water, and pipe wall [7]. Therefore, this paper proposes using the backward Euler implicit time integration scheme, hereafter called the BE scheme, which has more dissipative properties compared to 4th-RK scheme. In other words, numerical dissipation can represent physical dissipation.

The governing equations of the RC model for a frictionless horizontal pipe-reservoir system are Equation (1), without the last term in its RHS, and Equations (2) and (3). However, Equation (3) can be further simplified, using the continuity equation in Equation (2), as,

$$\frac{\mathrm{dp}}{\mathrm{dt}}=\mathrm{k}\frac{\mathrm{p}}{⩝}\mathrm{AV}.$$

(15)

By applying the BE scheme to the governing equations, the discretized form of the equations become,

$${\mathrm{V}}^{\mathrm{n}+1}={\mathrm{V}}^{\mathrm{n}}+∆\mathrm{t}\frac{\left({\mathrm{p}}_{\mathrm{R}}-{\mathrm{p}}^{\mathrm{n}+1}\right)}{\mathsf{\rho }\mathrm{L}},$$

(16)

$${\forall }^{\mathrm{n}+1}={\forall }^{\mathrm{n}}-∆{\mathrm{tAV}}^{\mathrm{n}+1},$$

(17)

$${\mathrm{p}}^{\mathrm{n}+1}={\mathrm{p}}^{\mathrm{n}}+∆\mathrm{tk}\frac{{\mathrm{p}}^{\mathrm{n}}}{{⩝}^{\mathrm{n}}}{\mathrm{AV}}^{\mathrm{n}+1},$$

(18)

where superscript $\mathrm{n}$ implies the current time step. Note that in Equation (18) the term $\mathrm{p}/⩝$ can also be calculated implicitly. However, it was realized that it does not have significant effects on the results. The details of the classical fourth order Runge–Kutta scheme can be found in other references, e.g., [18].

2.5. Effective Time Step Size

Rokhzadi and Fuamba [14], for problems with small ratios of the reservoir pressure to the atmospheric pressure, used a shock-fitting approach, which is a combination of the RC model and the Saint-Venant equations. Rokhzadi and Fuamba proposed using the BE scheme instead of the 4th-RK scheme. In this shock-fitting approach, the time step size is controlled by the Courant Friedrichs Lewy (CFL) condition associated with the gravity flow. Since by using the RC model there is no such criterion to control the time step size, so as the numerical dissipation, the question is how much numerical dissipation can be allowed on the calculation to avoid excessively imposing numerical dissipation and to prevent spoiling the solutions. In this section, a formula will be derived, by which it will be shown that the time step size is controlled by the physical variables, which affect the damping.

In order to find the numerical dissipation terms imposed by the BE scheme, the modified equation is needed. The modified equation can be found using Equations (16)–(18) and by substituting the Taylor series of the variables at the time step $\mathrm{n}+1$. Note that here, only those terms with $\mathrm{O}\left(∆\mathrm{t}\right)$ are kept and the terms with the higher order of $∆\mathrm{t}$ are neglected.

$$\frac{{\mathrm{dV}}^{\mathrm{n}}}{\mathrm{dt}}=\frac{{\mathrm{p}}_{\mathrm{R}}-{\mathrm{p}}^{\mathrm{n}}}{\mathsf{\rho }\mathrm{L}}-∆\mathrm{t}\frac{1}{\mathsf{\rho }\mathrm{L}}\frac{{\mathrm{dp}}^{\mathrm{n}}}{\mathrm{dt}}+\mathrm{O}\left(∆{\mathrm{t}}^2\right),$$

(19)

$$\frac{{\mathrm{dp}}^{\mathrm{n}}}{\mathrm{dt}}=\mathrm{k}\frac{{\mathrm{p}}^{\mathrm{n}}}{{⩝}^{\mathrm{n}}}{\mathrm{AV}}^{\mathrm{n}}+∆\mathrm{tkA}\frac{{\mathrm{p}}^{\mathrm{n}}}{{⩝}^{\mathrm{n}}}\frac{{\mathrm{dV}}^{\mathrm{n}}}{\mathrm{dt}}+\mathrm{O}\left(∆{\mathrm{t}}^2\right).$$

(20)

By comparing Equations (19), and (20) to Equation (1), without the friction term, and Equation (15), it can be realized that the first terms in RHS of Equations (19) and (20) are the physical terms involved in the governing equations and the second terms are the truncation error terms, which are responsible for imposing numerical dissipation. Substituting Equation (19) into Equation (20) yields,

$$\frac{{\mathrm{dp}}^{\mathrm{n}}}{\mathrm{dt}}=\mathrm{k}\frac{{\mathrm{p}}^{\mathrm{n}}}{{⩝}^{\mathrm{n}}}{\mathrm{AV}}^{\mathrm{n}}+∆\mathrm{tkA}\frac{{\mathrm{p}}^{\mathrm{n}}}{{⩝}^{\mathrm{n}}}\left(\frac{{\mathrm{p}}_{\mathrm{R}}-{\mathrm{p}}^{\mathrm{n}}}{\mathsf{\rho }\mathrm{L}}-∆\mathrm{t}\frac{1}{\mathsf{\rho }\mathrm{L}}\frac{{\mathrm{dp}}^{\mathrm{n}}}{\mathrm{dt}}+\dots \right).$$

(21)

As mentioned before, the energy equation, in which the effect of heat transfer is included, has a form as,

$$\frac{\mathrm{dp}}{\mathrm{dt}}=\mathrm{k}\frac{\mathrm{p}}{⩝}\mathrm{AV}+\left(\mathrm{k}-1\right)\frac{\mathrm{q}}{⩝}.$$

(22)

By comparing Equation (21) to Equation (22), an equation can be derived for the time step as,

$$∆\mathrm{t}=\frac{\mathsf{\rho }\mathrm{L}\left(\mathrm{k}-1\right)\mathrm{q}}{\mathrm{kpA}\left({\mathrm{p}}_{\mathrm{R}}-\mathrm{p}\right)}$$

(23)

Note that since Equation (23) was found by equating the heat transfer term to the numerical dissipation term, thus, these analyses imply that by applying the BE scheme to the governing equations of the RC model and using the time step size as in Equation (23), the physical dissipation of the heat transfer can be compensated by the numerical dissipation of the BE scheme. In addition, later, this compensation will be further explained. To clarify Equation (23), the heat transfer term needs to be further explained.

Following Lee [7], by neglecting the conduction and radiation heat transfers, the convection heat transfer can be calculated as,

$$\mathrm{q}={\mathrm{HA}}_{\mathrm{q}}∆\mathrm{T}={\mathrm{HA}}_{\mathrm{q}}\left({\mathrm{T}}_0-\mathrm{T}\right),$$

(24)

where $\mathrm{H}$ is the convection heat transfer coefficient, ${\mathrm{A}}_{\mathrm{q}}$ is the thermal area, ${\mathrm{T}}_0$ is the initial air pocket temperature, and $\mathrm{T}$ is the variable air temperature, both in the Kelvin scale. Following [7,19], the convection coefficient is calculated as,

$$\mathrm{H}=3.5{\left|∆\mathrm{T}\right|}^{1/3}.$$

(25)

The thermal area for the typical examples depicted in Figure 1 can be calculated as,

$${\mathrm{A}}_{\mathrm{q}}=2\mathrm{A}+{\mathsf{\pi }\mathrm{DL}}_{\mathrm{a}},$$

(26)

so that the heat transfer between the air, water, and pipe wall can be taken into account.

The air temperature is calculated using the perfect gas state relationship as,

$$\mathrm{T}=\frac{\mathrm{P}⩝}{{\mathrm{Rm}}_{\mathrm{g}}},$$

(27)

where $\mathrm{R}$ is the air constant $\left(\mathrm{R}=287\mathrm{J}/\mathrm{kg}·\mathrm{K}\right)$ and ${\mathrm{m}}_{\mathrm{g}}$ is the air mass, which is constant since there is no air release.

Equation (23) can be further expanded, using Equations (24)–(28), as,

$$∆\mathrm{t}=\frac{\mathsf{\rho }\mathrm{L}\left(\mathrm{k}-1\right){\mathrm{HA}}_{\mathrm{q}}∆\mathrm{T}}{\mathrm{kpA}\left({\mathrm{p}}_{\mathrm{R}}-\mathrm{p}\right)}.$$

(28)

As can be seen in Equation (28), the time step size depends on physical variables including the air and water initial lengths, the reservoir pressure, and heat transfer.

3. Results

To further discuss the effects of time integration schemes on the air pressure distribution and to clarify the effects of the reservoir pressure $\left({\mathrm{p}}_{\mathrm{R}}\right)$ and the product of the air and water initial lengths $\left({\mathrm{L}}_{\mathrm{a}}\mathrm{L}\right)$ as well as to compare the performances of the BE and 4th-RK schemes, some test cases, for which the experimental data exist, and a hypothetical example have been solved by the RC model, for which the governing equations are Equations (1), (2) and (15).

3.1. Hypothetical Test Case

A horizontal frictionless pipe-reservoir system with $\mathrm{D}=1\mathrm{m}$ and different air and water lengths and the reservoir pressure ratio ${\mathrm{p}}_{\mathrm{R}}/{\mathrm{p}}_{\mathrm{a}}=2.0$ has been solved and the results are shown in Figure 2 and Figure 3. Note that for both BE and 4th-RK schemes, the time step is set to $∆\mathrm{t}=0.001\mathrm{s}$.

As can be seen in Figure 2, as the product of the air and water lengths decreases, either by decreasing the air length or the water length, the rate of damping and the frequency increase. In addition, as indicated in [7] as well, Figure 2 shows that the first maximum pressure, for a frictionless pipe, is independent of the air and water initial lengths. Note that in Figure 2, it was attempted to examine large lengths of the air and columns as well (second row). Figure 3 shows the results of the same test case, in which the 4th-RK scheme is used instead of the BE scheme. As can be seen, by using the 4th-RK scheme, the RC model is unable to simulate the relationship between the damping and the product of the air and water initial lengths. However, the frequency and its relation with ${\mathrm{L}}_{\mathrm{a}}\mathrm{L}$ are properly predicted.

Similarly, to illustrate the effect of the reservoir pressure on the damping and the frequency and to compare the effect of time integration schemes, Figure 4 shows the non-dimensional air pressure distribution $\left(\mathrm{p}/{\mathrm{p}}_{\mathrm{a}}\right)$ in terms of time for a specific example with pipe diameter $\mathrm{D}=1.0\mathrm{m}$, and the air and water initial lengths ${\mathrm{L}}_{\mathrm{a}}=2.0$, and $\mathrm{L}=1.0\mathrm{m}$. The left graph illustrates the results calculated by the BE scheme, and the right graph is calculated by the 4th-RK scheme. As can be seen, the frequency increases with increasing the reservoir pressure, the feature that appropriately is predicted by both time integration schemes. However, from the right graph, it can be seen that the 4th-RK scheme is unable to appropriately calculate the damping and its direct relation with the reservoir pressure, while the BE scheme in the left graph effectively calculates the damping behavior. In addition, as shown in [7], the first maximum pressure, for initially stationary water flows, is only a function of the reservoir pressure and it increases with increasing the reservoir pressure, a fact that is calculated by both time integration schemes.

Therefore, it can be claimed that implicit time integration schemes, represented by the BE scheme, is more appropriate time integration scheme for using with the RC model.

3.2. Experiment of Zhou

Note that, hereafter, in every test case, the time step is calculated using Equation (28). One of the test cases is of Zhou [20], which is a reservoir-pipe system, in which $\mathrm{D}=0.035\mathrm{m}$, ${\mathrm{L}}_{\mathrm{t}}=8.96\mathrm{m}$, ${\mathrm{p}}_{\mathrm{R}}/{\mathrm{p}}_{\mathrm{a}}=2.43$, and $\mathsf{\lambda }=0.56$, and $0.89$, in which $\mathsf{\lambda }=\mathrm{L}/{\mathrm{L}}_{\mathrm{t}}$. In addition, the pipe is horizontal with a friction factor $\mathrm{f}=0.035$ and a local head loss coefficient ${\mathrm{K}}_{\mathrm{loss}}=0.093$. Figure 5 shows the air pressure distribution calculated by the 4th-RK and BE schemes and the corresponding experimental data for both ${\mathsf{\lambda }}^{\prime }$s.

First, it is shown that the physical solution approves the inverse relation between the damping and the product of the air and water initial lengths. The experimental data in Figure 5a shows that the first peak with the value of $230$ occurs at $\mathrm{t}=0.45\mathrm{s}$, and the second peak value of $164$ occurs at $\mathrm{t}=1.83\mathrm{s}$. Thus, the rate of damping in this test case is around $47.83{\mathrm{s}}^{-1}$. The same calculation in Figure 5b shows that the rate of damping is $102.5{\mathrm{s}}^{-1}$. These values show that when ${\mathrm{L}}_{\mathrm{a}}\mathrm{L}$ is smaller, as in Figure 5b, the rate of damping increases. Therefore, it can be claimed that the RC model and the relation derived for the amplitude of the oscillations in Equation (10) provide appropriate information about the physical behavior of the air pressurization.

In addition, the numerical solutions in Figure 5 show that the BE scheme outperforms the 4th-RK scheme. It is worth mentioning that the phase shift between numerical and experimental data was reported in other references, e.g., [5,6]. As Zhou et al. [5] explained, this phase shift could be due to the assumption that the air pocket remains intact, while the air roll ups and splits into several pockets with smaller sizes. However, for the test case with a smaller initial air length (Figure 5b), in the last period, the BE scheme causes a more obvious phase shift, compared to the 4th-RK scheme.

Here, it is shown that the numerical dissipation is imposed in proportion to physical conditions so that the performance of the RC model is improved. As demonstrated in [7], the dimensional analysis shows that the governing equations of the RC model have a similarity form with a scale factor ${\mathrm{fL}}_{\mathrm{a}}/\mathrm{D}$, which appears as a coefficient of the friction loss term. The details of the non-dimensional variables can be found in Equations (3)–(85) in reference [7]. Thus, since the only differences between the test cases of Zhou [13], as in Figure 5, are the initial air and water lengths, if the coefficient ${\mathrm{fL}}_{\mathrm{a}}/\mathrm{D}$ set equal for both cases, then the relation between the numerical dissipation imposed on the solutions and the product of the air and water initial lengths can be analyzed quantitatively. Therefore, the friction factor in the test case with $\mathsf{\lambda }=0.56$ is changed to $\mathrm{f}=0.035/4.0$, and the result of this test case along with Figure 5b, are shown in Figure 6.

The first peak calculated by the 4th-RK and BE schemes, for the test with $\mathsf{\lambda }=0.56$, shown in Figure 6a, are approximately $440$, and $331$, with a difference of around $109$. The same calculation for the case with $\mathsf{\lambda }=0.89$, in Figure 6b, shows that the first peak calculated by the 4th-RK and BE schemes are approximately $445$ and $362$, with a difference around $83$. These values show that the numerical dissipation imposed on the solutions of the test case with $\mathsf{\lambda }=0.56$, in which ${\mathrm{L}}_{\mathrm{a}}\mathrm{L}$ is larger, is larger than the test case with $\mathsf{\lambda }=0.89$, in which ${\mathrm{L}}_{\mathrm{a}}\mathrm{L}$ is smaller. Referring to Equation (15), it can be claimed that the BE scheme can improve the performance of the RC model because, according to Equation (15), when ${\mathrm{L}}_{\mathrm{a}}\mathrm{L}$ is large the damping of the RC model is small and the overestimation is large and, as seen in Figure 6a, the BE scheme imposes more dissipation to improve the overestimation.

3.3. Experiment of Zhou and Liu

Another experiment presented in this paper was carried out in [21], in which the air pressurization in a horizontal reservoir-pipe system was experimentally analyzed for different tailwater depths including the initially dry-bed condition. The experiment setup is similar to the schematic provided in Figure 1, with $\mathrm{D}=0.04\mathrm{m}$, ${\mathrm{L}}_{\mathrm{t}}=8.824\mathrm{m}$, ${\mathrm{L}}_{\mathrm{a}}=3.25\mathrm{m}$, and ${\mathrm{p}}_{\mathrm{R}}={\mathrm{p}}_{\mathrm{a}}+120\mathrm{kPa}$. The friction factor was calculated as $\mathrm{f}=0.075$, and, in the present paper, an additional factor representing the local head loss is considered as ${\mathrm{K}}_{\mathrm{loss}}=0.093$. The air pressure distribution of this test case, calculated by both BE and 4th-RK schemes, are shown in Figure 7a. Note that for this test case, $\mathrm{k}=1.2$ was used since it provides more accurate results than $\mathrm{k}=1.4$.

The experimental data in Figure 7a shows that the first peak with the value of $172$ occurs at $\mathrm{t}=0.8\mathrm{s}$, and the second peak value of $142$ occurs at $\mathrm{t}=2.0\mathrm{s}$. Thus, the rate of damping in this test case is around $25.0{\mathrm{s}}^{-1}$. Note that the product of air and water initial lengths in the test case is almost the same as the test case of Zhou [20], in which $\mathsf{\lambda }=0.56$ and the rate of damping was calculated as $47.83{\mathrm{s}}^{-1}$. Comparing the rate of damping in these two test cases shows that when the reservoir pressure is large, as in the test case of Zhou [20], the rate of damping is large as well. Therefore, the physical solution proves a direct relation between the damping and the reservoir pressure. In addition, it can be claimed that the RC model and the relation derived for the amplitude of the oscillations in Equation (10) provide appropriate information about the physical behavior of the air pressurization.

Moreover, it can be seen that the BE scheme improves the overestimated peak values compared to the 4th-RK scheme and no obvious phase shift is seen.

In order to further analyze the relation between numerical dissipation and the reservoir pressure, in Figure 7b, the reservoir pressure is increased to ${\mathrm{p}}_{\mathrm{R}}={\mathrm{p}}_{\mathrm{a}}+200\mathrm{kPa}$. The difference between the first peak value calculated by BE and 4th-RK schemes in Figure 7a, in which the reservoir relative pressure is $120\mathrm{kPa}$, is around seven. In Figure 7b, in which the reservoir relative pressure is around $200\mathrm{kPa}$, the difference between numerical peak values is around $18$. As the only difference between Figure 7a,b, is the reservoir pressure, it can be claimed that for cases with larger reservoir pressure, in which the overestimation of the RC model is larger, the BE scheme imposes more dissipation to improve the overestimation.

Therefore, considering the results of both test cases, it can be expected that implicit time integration schemes, represented by the BE scheme, are more helpful in using with the RC model. In addition, by using the derived formula for the time step, it can be ensured that excessive numerical dissipation is avoided.

3.4. Test Case of Lee

In this section, to show that the BE scheme by using the time step derived in Equation (28) can compensate for the effect of heat transfer, one of the test cases solved in [7] by 4th-RK scheme, in which the effect of heat transfer is included is compared to the result calculated by BE without the heat transfer. The test case is a frictionless horizontal reservoir-pipe system with $\mathrm{D}=1.0$ m, the initial water length $\mathrm{L}=40.0$ m and the initial air length ${\mathrm{L}}_{\mathrm{a}}=10.0$ m. The ratio of the reservoir absolute pressure to the atmospheric pressure is ${\mathrm{p}}_{\mathrm{r}}/{\mathrm{p}}_{\mathrm{a}}=2.0$. Note that following Lee [7], the variable water column length is considered in the RC model. As can be seen in Figure 8, the result calculated by BE scheme is almost the same as the result calculated in [7]. Thus, it can be claimed that the numerical dissipation of the BE scheme, by using the time step in Equation (28), can compensate for the effect of heat transfer.

4. Conclusions

This study analyzes the rigid column model to further explain that this model is able to produce the fundamental features of the physics of the problem, in particular, the damping behavior of the surge pressure distribution. However, this model underestimates the damping behavior and overestimates the peak values. The reasons are linked to simplified approximations associated with the rigid column model and the complexity of the transient flow, which cause this model to underestimate the prediction of the physical energy dissipation. To address this poor behavior of the rigid column model, using implicit schemes has been proposed, by which the numerical dissipation can represent the physical dissipation. In this regard, a criterion has been provided to control the amount of numerical dissipation to ensure that the solutions are not spoiled and, in the meantime, the performance of the rigid column model is improved.

The contribution of the friction loss term to the damping behavior of the surge pressure is known and discussed in the literature. Therefore, to provide an analytic formula by which the relation between the frequency and damping and physical parameters can be shown, the governing equations, applied to a frictionless horizontal reservoir-pipe system, were linearized around the point where the air pressure is equal to the reservoir pressure. This linearization implies that the air pressure does not deviate significantly from the reservoir pressure. Although this assumption is not valid for every example, it was shown that it can provide insightful information about the air pocket behavior. It should be noted that the linearized governing equations were just used to analyze the influential factors on the damping and frequency, while for solving the test cases, the nonlinear governing equations were used.

It was found that the damping and frequency of the air pressure, calculated by the RC model, has a direct relation with the reservoir pressure and inverse relation with the product of the air and water initial lengths. The same relationship was found in some available experimental data, meaning that it confirms that the RC model is able to predict the physical behavior of the air pressure distribution. In addition, it was shown that, despite the abilities of the RC model, the performance of this model depends on the type of the temporal scheme used for time integration, which is because of many simplified approximations associated with the RC model, including neglecting the heat transfer and lumped assumption of the air pocket. In this regard, it was shown that implicit schemes, represented by the backward Euler scheme, are more effective than explicit schemes, represented by the classical fourth order Runge–Kutta scheme. The reason is that the numerical dissipation associated with implicit schemes can compensate for the physical dissipation associated with the neglected factors in the RC model. Moreover, to avoid imposing excessive numerical dissipation on the solutions, a formula for calculating the time step size was derived, which is a function of the effective parameters on the physical damping and frequency. It was shown that, by using this time step, the backward Euler scheme can help the RC model to more effectively predict the air pressure distribution.