ABSTRACT

Most numerical transient flow models that consider dynamic friction employ a finite differences approach or the method of characteristics. These models assume a single fluid (water only) with constant density and pressure wave velocity. But when transient flow modeling attempts to integrate the presence of air, which produces a variable density and pressure-wave velocity, the resolution scheme becomes increasingly complex. Techniques such as finite volumes are often used to improve the quality of results because of their conservative form. This paper focuses on a resolution technique for unsteady friction using the Godunov approach in a finite volume method employing single-equivalent twophase flow equations. The unsteady friction component is determined by taking into account local and convective instantaneous accelerations and the sign of both convective acceleration and velocity values. The approach is used to reproduce a set of transient flow experiments reported in the literature, and good agreement between simulated and experimental results is found.

1. Introduction

Friction in pressurized flows can be decomposed into two components: static friction in steady flows and dynamic friction in transient flows. Transient flows present fairly important dynamic frictions that can significantly modify system behavior. Several types of model are used to calculate the dynamic friction component. One of the first types is the convolution-based model developed by Zielke [1], which uses the local acceleration and weight functions to calculate the unsteady friction component for laminar flow. The integration procedure in the Zielke model requires lots of memory and is very time-consuming. Others such as Suzuki et al. [2] and Schohl [3] suggested improvements in the determination of the unsteady friction component. Later, Vardy and Brown [4,5], Prashanth Reddy et al. [6] and Vitkovsky et al. [7] proposed an extension of the Zielke convolution model for smooth and rough pipes with turbulent flow. Another type of model called the instantaneous accelerations-based model is noted in the literature. It attributes the attenuation of flow amplitudes to local and convective instant acceleration [6,8,9]. In instantaneous accelerations-based model, the friction coefficient is decomposed into a static friction component and an unsteady friction component (Equation (1)).

(1)

with:

(2)

where is the Brunone coefficient, is the pipe diameter, is the water velocity, is the pressure wave celerity, is the time and is the abscissa.

Bergant et al. [10] improved this formulation by introducing instead of, to better reproduce the attenuation of flow acceleration and deceleration. For producing good accuracy in case of an upstream valve closing, second improvement was considered by the sign of the velocity and its acceleration or deceleration [11-13].

For solving transient flows with the dynamic friction, several numerical methods are proposed in the literature and two approaches are generally used. The first approach involves calculating the local and convective acceleration in the source term (see Equation (4)) as treated by Bergant et al. [10], Brunone et al. [14] and Bughazem and Anderson [15]. The second one transfers acceleration or one part of acceleration in the variables and the flux terms (see Equation (4)), as illustrated by Bughazem and Anderson [16], Wylie [17] and Vítkovský et al. [11]. The first approach uses finite differences, while the second employs the method of characteristics that is a graphical procedure for the integration of partial differential equations (PDEs) [18]. Prashanth Reddy et al. [6] used the second type of resolution with the characteristic equations developed by Vitkovsky et al. [11] to calculate the dynamic friction component. The formulation of the dynamic friction component in the instantaneous accelerations-based model has been tested with finite difference techniques and the method of characteristics. However, there is a lack of literature on the use of the finite volume technique in instantaneous accelerationsbased models.

Finite volume techniques, particularly those using the Godunov method (as presented by Toro [19] and Guinot [20]) are now widely used for transient flow in sewers (as in León et al. [21], Vasconcelos and Wright [22] and Sanders and Bradford [23]) but without unsteady friction. This paper presents and applies the instantaneous accelerations-based model including a numerical method based on Godunov approach. Four specific objectives are targeted: i) to modify the single equivalent two-phase flow equations, ii) to consider air content, i.e. some compressibility, in the equations, iii) to propose a numerical method, and iv) to compare numerical and experimental results.

2. Methodology

The first section of the paper proposes the governing equations and the second a numerical method adapted from Guinot’s [20] traffic flow model. The third section presents two case studies in which numerical and experimental results are analyzed. The proposed model is governed by the single-equivalent two-phase flow equations. The dynamic friction component refers to the instantaneous accelerations-based model formulation and the resolution scheme uses the Godunov approach. This formulation of equations considers some compressibility of the fluid with variable density and pressure wave celerity. The model results are compared to experimental measurements from Adamkowski and Lewandowski [24] and Bergant et al. [25].

3. Governing Equations

The pressurized flow Equations (3) and (4) are those presented by Guinot [20] for the calculation of waterhammer in the presence of a two-phase fluid with an air content.

(3)

Vectors, and correspond, respectively, to the variables, flux and source term defined by Equation (4).

(4)

where: and represent the time and the abscissa, the flow cross-section, the flow discharge, the fluid density, the mass of fluid per unit length of the pipe, the mass discharge, the pressure, the pipe slope, the friction slope calculated by Dary-Weisbach formula as following, acceleration of gravity, and the water velocity.

The friction coefficient is determined by considering local and convective instantaneous accelerations as in the model by Brunone et al. [8,9] that was modified by Bergant et al. [10] as in Equation (5).

(5)

where the Brunone coefficient depends on expressed in Equation (6).

(6)

Friction is therefore first decomposed into static and dynamic components as shown in Equation (5). The dynamic friction component is transformed according to the variables flow and in Equation (7):

(7)

In which:

(8)

The Brunone coefficient is determined as in Bergant et al. [10] by:

(9)

where (Equation ) depends more on the flow regime (Reynolds number).

(10)

The contributions of the dynamic friction acceleration in the vector variable and in the flux vector (Equation (4)) allows to obtain Equation (11), in which is the friction slope due to the dynamic friction component.

(11)

In vectorial form Equation (11) become Equation (12), where the sought variables are and. The other variables will be determined by considering the Equations (24), (25) and.

(12)

with:

(13)

and

(14)

Equations (12) and (13) are equivalent to Equations (3) and (4) if the dynamic friction component is not considered (i.e.,).

The next proposes a numerical method for solving these equations for each internal cell. The exterior cells will be determined by the boundary conditions formalized in each case. Since a variable is known for each pipe end, one of the Riemann invariants will be used to calculate the second variable. For example, for a known discharge or pressure to the left boundary, the second Riemann invariant is used to calculate the second variable (discharge or pressure). When it is the right boundary, the first Riemann invariant is used. Details of the calculation will be discussed in the next section.

4. Procedure for Numerical Solution of Partial Differential Equations (PDEs)

Equations (12) and (13) resemble the Guinot’s [20] traffic flow model, in which the acceleration/deceleration of a vehicle influences the speed of the preceding vehicle. Solving this non-conservative system of balance equations with the Godunov approach by using the solution of the Riemann problem, involves four steps:

The first step is to make the solution of the Riemann problem of Equation (15) to express the conservation component of Equation (12):

(15)

The second is to analyze the solution of the Riemann problem of the second hyperbolic part (Equation (16)):

(16)

The third is to obtain the solution of the Riemann problem of the source term corresponding to Equation (17):

(17)

The fourth step is to obtain all parameters (as velocity and pressure in each cell) needed to determine the flow. The next time step is calculated by taking into account the Courant condition:, corresponding to the spatial step.

4.1. Step 1: Solution of the Riemann Problem for the Conservation Part

The conservative part is close to the two-phase flow of the equation presented by Guinot [20] and León et al. [21]. The only difference is the factor in the second line of the flux regarding Equation (13). If the dynamic friction component (i.e.,) is not considered, the equations are those usually used when only the static friction component is considered. Considering, the Jacobian matrix of respecting the matrix, can be found in Equation (18):

(18)

After processing and arrangement, the eigenvalues of the matrix are given by Equation (19).

(19)

With:

(20)

The eigenvalues, corresponding to the eigenvectors for the matrix, are and as seen in Equation (21).

(21)

4.1.1. Riemann invariants of the Conservative Solution

The Riemann invariants along each characteristic are given by Equation (22).

(22)

Equation (22) can be integrated respectively between and following approximation according to the trapezoidal rule, as with static friction. represents the variables to the left of the region of the constant state, the region of the constant state and those to the right of the constant state region. The resulting equation (Equation (23)) is solved by the Newton Raphson method for determining i.e. and.

(23)

Note that the celerity, which depends on, is related to other parameters through Equation (24), as presented by Guinot [26].

(24)

With fluid density, the volume fraction of air at the reference pressure, the polytropic coefficient of air (for isothermal conditions and for adiabatic conditions). The pressure waves celerity in water is calculated as suggested by Wylie and Streeter [27].

The solution of Equation (23) yields: and i.e. and used for calculating the mass discharge and the vector variable. The pressure force is obtained by solving Equation (25) [20].

(25)

4.1.2. Boundary Conditions Calculation of the Conservative Part

The approach assumes half-virtual cells at and as suggested by Guinot [20], using a reference pressure/velocity depending on the boundary conditions type.

Prescribed pressure at the left boundary

If the pressure is prescribed, i.e. known at the left boundary, then Equation (25) is used to determine. Therefore Equation (24) can be used to determine from. The Riemann invariant along (Equation (23)) yields the velocity(Equation (26)) with corresponds to flow parameter at cell 1 at time.

(26)

Prescribed discharge at the left boundary

If the discharge is known, then the velocityRiemann invariant (Equation (27) alongis combined with Equations (24) and (25) yields and.

(27)

Prescribed pressure at the right boundary

If the pressure is prescribed at the right boundary, Equations (25) and (24), and Riemann invariant along (Equation (28) or (29)) provide and with corresponding to flow parameter at cell at time.

(28)

(29)

Prescribed discharge at the right boundary

When discharge is prescribed, Riemann invariant (Equation (28)) along, combined with Equations (24) to (25), yields, and as shown in Equation (30):

(30)

4.1.3. Balance over the Cells

The conservative part yields the homogeneous solution (Equation (31)) at each cell without the source term as presented by Guinot [20]:

(31)

The flux will be calculated by Equation (32):

(32)

The provisional values will be corrected in two successive stages, integrating the second hyperbolic part and the source term.

4.2. Step 2: Solution of the Riemann Problem of the Second Hyperbolic Part

Equation (16) is solved by seeking the eigenvalues and vectors of (Equations (33)).

(33)

Given that eigenvalues and vectors of differ depending on the sign of, each case of Equations (33) is treated separately in sections 4.2.1 for and in section 4.2.2 for.

4.2.1. Resolution of Case 1

Eigenvalues and eigenvectors of for Case 1 are depicted in Equations (34) and (35), respectively:

(34)

(35)

Considering the validity of the eigenvectors in Equation (35), along each characteristic line, solution of the Riemann problem (Equation (36)) is obtained.

(36)

Using the trapezoidal rule, as in Equation (37), the velocity of the constant state is calculated by Equation (38).

(37)

(38)

After the determination of and, the vector variable is determined for each inter face. The new variable considering the second hyperbolic part is then calculated by applying the Riemann invariants principle [19,20]:

(39)

4.2.2. Balance over the Cells

The following formulation (Equation (40)) is obtained when is integrated over the cell for.

(40)

4.2.3. Resolution of Case 2

In Case 2, eigenvalues (and) and eigenvectors (and) of are respectively Equations (41) and (42):

(41)

(42)

As in Case 1, the Riemann invariant is expressed by Equation (43).

(43)

The second part of the solution of the Riemann problem (Equation (43)) is used to calculate the velocity of the constant state (Equation (44) or (45)).

(44)

(45)

The quantities and yield the vector variable

.

The new variable taking into account the second hyperbolic part is then calculated by applying the principle of Riemann invariants:

(46)

4.2.4. Balance over the Cells

Integrating Equation (46) over the cell yields for:

(47)

4.2.5. Boundary Conditions

The last cell is calculated by considering either the first component of Equation (37) for or the second component of Equation (44) for. Each of the resulting Equations (48) or (49) is combined with the known flow condition:

(48)

(49)

with dependant respectively on and.

4.3. Step 3: Solution of the Riemann Problem of the Source Term

The source term is solved by Equation (50) to estimate the variable vector of each cell depending on, as presented by Guinot [20] and Toro [19].

(50)

where indicates that is used to evaluate the source term of the Equation (13).

4.4. Step 4: Flow Parameter Calculation

After determining the variables and, velocity, density, pressurewave celerity by Equation (24), and pressure by Equation (25) are determined for each cell.

5. Results and Discussion

Two case studies ((1) a closed downstream valve and (2) a closed upstream valve) are performed, and numerical results compared to the experimental results.

5.1. Case Study 1: Closed Downstream Valve Analysis

5.1.1. Experimental Setup of the Case Study 1

The used experimental results are from the study of Adamkowski and Lewandowski [24]. The experiments were conducted at a test rig composed of a 98.11 m long copper pipe with an internal diameter of 0.016 m and a wall thickness of The pipe slope is over the horizontal plan. A steel cylinder with a diameter of about 1.6 m is used as an upstream reservoir. A quick-closing ball valve is installed at the downstream end of the pipe. Four absolute pressure semiconductor transducers are mounted on the pipe at 0.25 L, 0.5 L, 0.75 L and L from the reservoir. The test procedure consists in an instant closing of the valve. During the tests, the steady head-water level in the upstream reservoir is. The water temperature is 22.6˚C and the kinematic viscosity coefficient is m²/s. The pressure wave celerity calculated according to the pipe characteristics is m/s. The initial velocity in the pipe before closing of the valve is.

5.1.2. Discussion of Case Study 1 Results

Comparison of numerical results with measured results needs to consider several parameters such as the rate of air (Equation (24)) and the static friction. As shown by Wylie and Streeter [27], the rate of air reduces pressure wave celerity. This reduction is reflected in the numerical model by an increase in the pressure-oscillation period. This analysis was used to estimate the rate of air. Simulations were performed to find the value of the rate of air for better correspondence between measured and calculated oscillations. The rate of air is then considered to be 0.02%. The static friction is selected from Axworthy et al. [28], who compared their numerical results to the experimental results of Adamkowski and Lewandowski [24]. The polytropic coefficient is considered as. After several simulations, this value, corresponding to adiabatic conditions, seemed to yield best results.

Comparison between simulated and measured results (Figure 1) shows good correlation with the frequency of oscillations. The period of frequency is slightly larger in the case without taking into account the dynamical friction component. There is a slight overestimation of pressure when the dynamical friction is not taken into account, which may reach a maximum of 3 m at the downstream end of the pipe. For pressure amplitudes, the difference between measured and simulated results increases over time. This indicates that the dissipation of the numerical model is not sufficient; in other words, the friction coefficient is underestimated. This may be due either to the static friction component or the dynamic friction component. Moreover, the increased gap between numerical results and measurements may be due to the Brunone coefficient calculation by Equations (9) and (10). Indeed the difference becomes important beyond 2.5 s, i.e. when the Reynolds number in the peak pressure falls below 6000. This shows that when the flow dissipates or becomes less turbulent, the coefficient becomes less accurate.

5.2. Case Study 2: Closed Upstream Valve Analysis

5.2.1. Experimental Setup of the Case Study 2

The experimental tests used to examine the unsteadyflow, in case of a closed valve installed upstream just after the pump are taken from Pezzinga and Scandura [29]. These results are extracted from the paper by Prashanth Reddy et al. [6]. The tests were carried out on experimental installation composed of a zinc-plated steel pipeline (length 143.7 m, diameter 53.2 mm, thickness 3.35 mm, modulus of elasticity 2.06 × 10¹¹ N/m², roughness 0.1 mm) fed by a centrifugal electropump. A 1 m³ pressure tank is located at the downstream end of the pipe. Closing the upstream valve generates an interesting transient flow useful for analyzing dynamic friction.

The pressure variation is measured by pressure transducers located at the upstream valve and at the middle of the pipe. The average temperature of the water during the tests was 15˚C; the values of the kinematic viscosity and of the elastic modulus K were assumed to be 1.14 × 10⁻⁶ m²/s and 2.14 × 10⁹ N/m², respectively. The theoretical pressure wave celerity considered by the authors [29] is equal to 1360 m/s. As with the downstream valve in Case Study 1, a low rate of air into the water is assumed. Several simulations are performed with different rates of air to find out which allows for a better comparison between calculated and measured pressure oscillations. The initial velocity in the pipe, before closing of the valve, is, and the static friction component is considered equal to.

5.2.2. Discussion of Case Study 2 Results

Figure 2 shows the comparison between calculated and measured pressure values. The best superposition of oscillations is obtained with a rate of air of 0.01% and a polytropic coefficient considered equal to. The calculated frequencies of oscillations coincide perfectly with those measured at the valve and at the middle of the pipe. Consideration of the static friction component or dynamic friction component shows that they tend to dissipate energy. The inclusion of the dynamic friction

component enables a better agreement between calculated and measured pressure values. The maximum difference obtained is 6%, and is greatest at the valve position. This difference may be due to the boundary condition calculation. The difference between measured and calculated pressure values is very low in the first pressure peaks, i.e., just after the valve closure. These pressure peaks are the most damaging to water plants. Overall, the numerical results tend to overestimate pressure values, especially at the valve position, which may be due to the choice of static friction, the initial flow condition before the valve closure, the boundary condition calculation, or the early stage of computing. Another factor that might affect the quality of results is the choice of the rate of air. The choice of the polytropic coefficient may also influence the quality of results. Experiments are rarely realized under ideal adiabatic or isothermal conditions. Thus, tests with a uniform measure of the rate of air are needed to better calibrate the model.

6. Conclusions and Recommendations

This paper focuses on the resolution of unsteady friction using the Godunov approach in a finite volume method with single-equivalent two-phase flow equations. The calculated results show good agreement with experimental measurements, especially for the first pressure peaks, which are the most dangerous for pipe safety. However, dissipation is found to be lower in the calculated than in the measured results. The differences between simulated results with a dynamic friction component and those with a static friction component are more important in the case of a closed upstream valve than a closed downstream valve. It would appear that taking the dynamic friction component into account is more relevant in the case of upstream valve closure. This could be due partly to the formalization of the dynamic friction component, and in particular, the inclusion of acceleration and deceleration, as well as direction of flow. The difference could also be due to the calculation of boundary conditions, which are different in the two case studies. These results demonstrate that the proposed approach allows dynamic friction to be taken into account in finite-volume models, using the increasingly popular Godunov approach. The results also point to the possibility of considering the effect of air in order to improve the quality of simulation models. The model requires further improvements to more accurately reproduce the pressure values in cases of upstream valve closure. Searching for an adequate solution should address the determination of boundary conditions, the formalization of the dynamic friction component and the choice of the polytropic coefficient and experiments in which all parameters are adequately defined.

REFERENCES