Simulations and measurements of friction in oscillating flow

Rapid changes in the operating point of a hydropower plant cause transients in the hydraulic system. The flow will oscillate between the surge shaft and the reservoir until the oscillations are damped out by the friction. This paper investigates different friction models for oscillating flow, such as a quasi-steady model, Vítkovský’s model, and a one-term model. The governing equations for transient flow were solved by using rigid water column theory. The simulations were compared with experimental results from a small-scale test rig. The test rig consisted of a small reservoir, a horizontal pipe, and a surge shaft. The oscillations were induced by closing a valve downstream of the surge shaft. Measurements were obtained with pressure transducers.


Introduction
Changes in the operating point of a hydropower plant induce transients in the hydraulic system. Transients are often split into fast transients such as the water hammer and slow transients such as mass oscillations [1]. The friction in such transitional, unsteady flow is different than the friction in steady flow [2]. But few models are able to accurately account for the friction loss in transient flow, while at the same time being computationally efficient and easy to use. This paper investigates a few different friction models for the slow transient mass oscillations and compares them with experimental results.
Fluid transients have been studied since the 17th century [2]. But the equations governing transient flow were not fully established until the 1960s [3].
For slow transients, the governing equations are often simplified. Rigid water column theory in which the pipe walls are considered rigid and the water is considered incompressible, is often used for mass oscillations since it is simpler than taking the elasticity into account and still captures the physics quite well [4].
Conventionally, the friction relations for steady flow have been expected to hold at every instant in time in transient flow analysis [3]. The expression for the head loss in steady flow was proposed by Weisbach in 1850 and is well-known today as the Darcy-Weisbach equation [5]. The head loss is caused by viscosity, and is related to the wall shear stress [6]. A quasi-steady approach is sometimes used, in which the steady friction relation is updated at each instant in time. However, the Darcy-Weisbach head loss relation is usually not satisfactory for transient flow, since flow reversal will occur in transient flow and due to the no-slip condition at the pipe wall, this will yield larger wall shear stresses than in steady flow [3]. Indeed, research has shown that using steady or quasi-steady relations results in discrepancies between experimental and numerical results [3]. An unsteady component is therefore often included in the friction expression for transient flow.
Daily et al. [7] found that accelerating flow had a slight increase in frictional resistance compared with steady flow, while decelerating flow had a slightly smaller frictional resistance. They argued that this was because during acceleration the velocity profile steepens and gives higher shear. They developed an expression for the unsteady friction, in which it depends on the instantaneous mean velocity, V , and the instantaneous local acceleration, ∂V /∂t. Brunone modified the model proposed by Daily et al. and included the instantaneous convective acceleration, ∂V /∂x [8].
Vítkovský introduced a sign term and an absolute value around the instantaneous convective acceleration in Brunone's formulation of the friction model in order to ensure that the convective term has the correct sign for all cases [8]. Since his investigations had shown that the original formulation gave the wrong sign of the convective term, −a∂V /∂x for some cases, such as closure of upstream valve where initial flow in the pipe is in positive x-direction [8]. Bergant et al. [8] showed that Vítkovský's model gave a better fit than the quasi-steady model when compared with measurements of an upstream valve closure when looking at the first second following the valve closure. Even though Brunone's friction model was developed for the fast transient water hammer, it is expected to perform well for mass oscillations as well, since many of the flow characteristics which occurs in water hammer flow such as flow reversal, also occurs in mass oscillations.
Bergset [9] proposed a new one-term model which builds on the quasi-steady friction model. In order to correct the quasi-steady friction model, he proposed to multiply the friction term by a constant which had the nondimensionalized convective acceleration and instantaneous local acceleration terms in the exponent.

Governing equations
The equations governing transient flow in pipes with circular cross-section can be derived from conservation of mass and continuity as shown in for example [10]. The resulting continuity and momentum equation for transient flow, expressed in terms of the piezometric head, H, are equation (1) and equation (2), respectively.
These equations are valid for slightly compressible fluids flowing at low Mach numbers. H is as mentioned the piezometric head, a is the speed of sound, g is the gravitational acceleration, V is the fluid velocity, and h f is the head loss per unit length. For transient flow, the head loss is often split into a steady and unsteady part, h f = h f,s + f h,u . The speed of sound, a, for a thick-walled or rigid pipe is a = K/ρ where K is the bulk modulus of elasticity of the fluid, and ρ is the density of the fluid [10].

Methods of Solutions
The partial differential equations governing transient flow, equation (1) and (2), can be solved by different numerical methods. The method of characteristics is the most common method used when investigating water hammers [10]. For slow transients, such as mass oscillations, it is common to solve the equations using rigid liquid column theory. Under the assumption 3 that mass oscillations act as a rigid liquid column, the partial differential equations reduce to a set of ordinary differential equations. With rigid liquid column theory, the fluid is considered incompressible, and the pipe walls are considered rigid. This implies an infinite speed of sound [10]. Consequently, at any point in the system, a flow change is transmitted instantaneously throughout the system, and the water will have a common fluid particle velocity along the pipe [2]. With the assumptions above, the partial differential equations reduce to ordinary differential equations since the flow variables do not vary with distance, only with time [2]. The continuity equation becomes and the momentum equation is reduced to where k = f L 2gA 2 D with L being the length of the pipe. A is the cross sectional area of the pipe, Q is the flow rate, and z is the height of the free surface in the surge shaft. Equations (1) and (2) are thus reduced to an initial value problem consisting of two ordinary differential equations, equation (3) and (4) with two unknowns, z and Q, and two known initial values, z 0 and Q 0 . Solutions of this can be found by different numerical methods, such as the Euler method, a modified Euler method or a Runge-Kutta method. In this paper, the Euler method was chosen. For equation (3) and equation (4), the Euler method yields with f (Q n ) = gA L (z − kQ n |Q n |), and This method is first-order accurate in time, meaning that the error is of order ∆t. Thus, ∆t should not be chosen too large.

Friction models
The steady component in the expression for the head loss is simply taken as the head loss in steady flow. This is well-known as the Darcy-Weisbach equation. Expressed per unit length, the steady part of the head loss in a pipe with diameter, D, is where the friction factor, f , can be found from Moody's diagram or a formula approximating Moody's diagram such as Haaland's or Colebrook's equation. For the simulations presented in this paper, Haaland's equation was used. The quasi-steady friction model uses head loss relation for steady flow, but updates the flow variables with each time step. The unsteady term is zero in the quasi-steady approach [3]. Different expressions for the unsteady friction component have been proposed. Vítkovský's formulation of the unsteady friction term takes the sign of the convection term in Brunone's model into account, and is expressed per unit length as follows [8] Bergset [9] developed a one-term model for the friction loss based on the quasi-steady model. The friction expression was modified by a term taking the velocity change into account, as seen in equation (10). The model was then implemented in the same way as the quasi-steady model, by updating the friction factor to correspond with the current Reynolds number. Per unit length the proposed head loss was as follows The proposed expression for the constant B was where the subscripts 1 and 2 denote sections of a pipe with different diameter.

Experiment
A drawing of the test rig which was used for the experiments can be seen in figure 1. The rig consists of an upper reservoir, a 21 m long headrace tunnel, a surge shaft and siphon system. The reservoir has a spillway which ensures that the water level in the reservoir stays constant Measurements were done with three UNIK 5000 pressure transducers from GE Druck, located at three different places of the rig as seen in figure 1. An electromagnetic flowmeter from Krohne, called OPTIFLUX 2300C, was used for measurements of volume flow. Data acquisition from the measuring instruments was done with a NI USB-6211 device, and sampling was done at a frequency of 10 Hz.
The fluid in the experiment was water with an assumed dynamic viscosity of 1.138e-3 kg/(ms). Initial flow rates ranging from 0.0033 m 3 /s to 0.0015 m 3 /s were tested. This corresponds to flow with initial Reynolds number ranging from 24592 to 11178. Valve 1 through 3 were open initially, while valve 4 was adjusted to control the flow rate. The mass oscillations were induced in the system by manually closing a butterfly valve downstream of the surge shaft, valve 3 in figure 1.

Results
The simulations were done in MATLAB. The Euler method was used to solve the initial value problem, with initial conditions for the flow rate as given by the measured flow rate before valve closure. Different time steps were tried, and the results presented here are from the simulations with time step ∆t = 0.01.
The experimental results for the different initial flow rates showed the same tendencies when compared with the existing friction models, so results for only one flow rate is presented here, flow rate 0.0033 m 3 /s, corresponding to Reynolds number, Re = 24592.
The simulations with the quasi-steady friction model, Vítkovský's friction model and the oneterm model is plotted together with measurements of the pressure from one of the experimental runs, in figure 2a, 2b, and 2c, respectively. The pressure measurements are from the third pressure transducer, PT3, in figure 1, situated in the surge shaft.
The three friction models performs quite similarly. They all model the first peak accurately. But with time, none of the models provide enough damping compared with the damping measured in the experiment. Vítkovský's model performs marginally better than the quasi-steady model. This is difficult to see from plots in figure 2a and 2b, however when the difference between the simulated results and the measured result is studied at each local extrema, one can see that Vítkovský's model provides a slightly larger frictional loss than the quasi-steady model. This can be seen in table 1 where the mean and largest absolute difference between the measurements and the simulations are presented. The quasi-steady model as seen in figure 2a and the one-term model as seen in figure 2c turns out to be very similar since the correction term in the one-term model becomes approximately 1. The coefficient B has an exponent which is on average of the order e-6 during the time period for the oscillations. In this project, B, was found in the same way as Bergset proposed in [9]. With this method, the B will be completely dependent on the diameter change. For this rig there is no diameter change, and the diameter D 2 was arbitrarily chosen to be the same as Bergset used.  Simply choosing a different D 2 will therefore yield completely different results. The D 2 value was set to 0.325 m, and with the flow rate from this experiment, 0.0033 m 3 /s, this resulted in a B equal to 8.1250e+05. However, as the exponent is so small the value of B doesn't really matter much in the case tested, as the whole correction term will tend towards 1.
In an attempt to better match the one-term model to the experimental results, different modifications were tested. As the exponent in the correction term is so small, simply increasing the value of B did not help much. The best way that was found was to multiply the friction loss by a factor larger than one. Through trial and error it was found that multiplication by 3.2 yielded the best results. This modification of the one-term model can be seen in figure 2d.
After this modification, the model yields too much damping in the beginning and for about the first 140 s, while it does not yield enough damping towards the end. As for the other three simulations, the pressure waves are not fully damped out.
The modification of the one-term model reduces the mean absolute difference between the simulations and the experimental results by nearly 74 % compared with the original one-term model. The largest difference between the simulation and experimental results in a single extrema is reduced by 47 %. The modification improved the fit of the one-term model compared with the experimental results, however the fact that it yields too much damping in the beginning is undesirable. Thus, further investigations of how the friction models can be modified to accurately predict the friction loss for oscillating flow should be done.

Conclusion
Simulations with the quasi-steady friction model, Vítkovský's formulation of Brunone's friction model and the one-term friction model have been compared with experimental measurements. The models perform quite similarly. The model by Vítkovský performs only slightly better than the other two. They all match well for the first peak, but for the following peaks they do not yield enough damping of the oscillations. A modification of the one-term model was also investigated, in which the friction term was multiplied by a constant. This resulted in more damping, however the fit for the first peaks are then not satisfactory as the model yields too much damping for the first peaks.

Further work
Further investigations of the characteristics of transient flow should be done. More measurements of mass oscillations should therefore be performed, and finding measurement techniques that can give more detailed insight into the flow characteristics is essential. For example, particle image velocimetry, laser Doppler anemometry or constant temperature anemometry could be considered as measuring techniques which could increase the knowledge of the velocity profile and the shear stress in slow transient flow.
In this paper, the governing equations for transient flow were simplified by assuming that the mass oscillations act as a rigid liquid column, thus neglecting the compressibility of the water and the elastic effects of the pipe. Other methods of solving the partial differential equations such as the well-known method of characteristics could therefore be used instead, as this takes the compressibility and elasticity into account.
The unsteady friction models do not yield enough damping of mass oscillations, and further research must be done to find good modifications of the models so that they match better with experimental measurements. In this case, modifications of the friction expressions were only investigated for the one-term model. Modifications of the quasi-steady friction model and the Vítkovský friction model should also be investigated. The one-term friction model was so far only modified to match the experiment with Re = 24592. The modification that was found should also be tested and compared with experimental results from different flow conditions, in order to see how the constant value is affected by the flow conditions in the pipe, and hopefully find some correlation between the value and the changes in flow conditions.