Influence of Second Viscosity on Pressure Pulsation

A mathematical model of pulsating flow is proposed in the paper. The model includes more accurate description of energy dissipation, so it allows, for example, better stability analysis of water power plant control and more effective operation. Flow in a pipeline system is usually treated as a one-dimensional flow. This is also applied for more difficult cases of the Newtonian and non-Newtonian liquids simulations in the rigid or flexible pipes. Computational simulations of pressure pulsations in pipelines often predict lower damping than what the experimental results show. This discrepancy can be caused by neglecting one of the important damping mechanisms. The second viscosity describes the energy losses due to the compressibility of the liquid. Its existence and use in the computations specifies the real pulsations damping descriptions and predictions. A frequency dependent model of pressure pulsations including second viscosity is introduced. The second viscosity is determined from the system eigenvalue. The experiments were performed with water for low frequencies (from 0.1 to 1 kHz). This area is not fully covered by the current available research results.


Introduction
Fluid flow in pipeline systems is usually treated as a one-dimensional flow because one component of velocity is dominant; thus, other components can be neglected. Considering the one-dimensional flow assumption, computational simulations of transient events, pressure pulsations or water hammer often suffer from inaccuracy of the mathematical description of the dissipation mechanism, which is usually derived from the energy loss of steady flow and is proportional to velocity squared.
Usually, real energy dissipation is stronger than numerical models predict. This is obvious when the method of characteristic is used. Some authors compensate this difference by introducing the unsteady friction model, which is a type of turbulent model. Daily et al. in [1], wrote a mathematical description of the basic turbulent model and performed an experimental investigation of the energy loss of an accelerating or decelerating flow, and this work is still often cited.
Caulk [2] presented a one-dimensional nonlinear equations model derived for axisymmetric motion inside a slender surface of revolution. The general equations were specialized to reflect specific conditions on the fluid, such as the presence of surface tension, a confining elastic membrane or a rigid tube. Robertson [3] applied the nine-director theory on the modeling of the blood flow in walls of arteries, that are typical case of flexible wall and non-Newtonian liquid fluid structure interaction case. Carapau [4] used the director theory approach related to fluid dynamics for reduction of the nonlinear three-dimensional equations governing the axisymmetric unsteady motion of a non-Newtonian incompressible fluid to a one-dimensional system of ordinary differential equations depending on time and on a single spatial variable. Unfortunately, the model of unsteady friction is not as well described as the model of steady friction, and investigation in this area is still ongoing [5][6][7]. In contrast, the inclusion of second viscosity in the model instead of unsteady friction makes results more accurate and it has an advantage: second viscosity is associated with fluid matter, and it is possible to use it in the frequency domain.
The energy dissipation due to second viscosity occurs when the fluid is being compressed or when the fluid is expanding. Stokes formulated a hypothesis that second viscosity equals negative two thirds of the shear viscosity η 2 = − 2 3 η; thus the fluid is supposed to be in thermodynamic equilibrium [8]. It has already been proven that Stokes' hypothesis is not valid for polyatomic gases and liquids [9]. The author showed the influence of second viscosity on the absorption and dispersion of sound in fluid and discussed the frequency dependence of the second viscosity.
Examination of the theoretical model of liquid was performed in [10], where authors noted the role of second viscosity in the spreading of wave fronts, dilatation and compression of fluid or attenuation of sound waves. Numerical determination of the bulk viscosity of water was also conducted in [11,12].
When we look for experimental data of second viscosity, we find out that they are very rare. It is very difficult to find any experimental data within the last 150 years. Regarding water, possibly the oldest experiments were performed by Karim [13]. Foldyna published the dependence of second viscosity on pressure [14,15]. In addition, [16] provided measured data for water, and [17,18] measured the temperature dependence of the second viscosity.
All published experiments were performed at high frequencies (on the order of MHz or GHz) in a narrow range.
A method of the second viscosity determination from the damping of pressure pulsations is described within this paper. The suggested procedure allows determining of the second viscosity at low frequencies.

Model Derivation
This chapter describes the derivation of a mathematical model of pulsating flow in a pipeline, when second viscosity is assumed. The influence of an unsteady velocity profile is considered as well. The solution is obtained in the frequency domain.
The pulsations can be simulated using a transfer matrix that describes system properties, which include the fluid and pipe. The derivation of the transfer matrix originates from the continuity equation and momentum equation for a one-dimensional flow, which can be assumed for the pipe.

Continuity Equation
The continuity equation for liquid can be written as follows: where relation between density and pressure for barotropic fluid dρ = dp a 2 is included [19], which is a common practice for liquids. Integration over the liquid volume dV (see Figure 1) yields: Flow rates Q 1 and Q 2 are associated with the corresponding cross-sections S 1 and S 2 . Pressure p is now the mean value over the entire volume V of liquid at one time instant. Considering Q 2 = Q 1 + ∂Q ∂x · dx, it is possible to write Equation (2) in the form of Equation (3), where the volume dV = S · dx: Because the flow in a pipe is being solved, only the one dimensional case can be considered. Thus, the space coordinate x has only one dimension that is associated with the pipe axis. The integral over the pipe inner surface P equals zero when the pipe wall is rigid and leak-proof. However, usually, pipe elasticity allows pipe wall deformations and reduces the sound speed in the system.
The behavior of the pipe wall can be described using the Kelvin-Voigt model; see Figure 2. Using force balance and deforming properties of the model, it is possible to derive: where subscripts refer to the corresponding parts of the body in Figure 2. For a polymeric material pipe a more complex Kelvin-Voigt model must be used. For this case it is complicated to find the appropriate material constants. A detailed discussion of this problem is written in [20,21].
When the Laplace transform with zero initial conditions is applied to Equations (4) to (6), one can obtain:σ where the expression in parenthesis is complex Young's modulus E c .
In the case of a thin wall, the stress and velocity of the wall can be computed according to the following equations: and after the Laplace transform:σ The combination of Equations (7), (10) and (11) yields a relationship between the wall velocity and pressure in the pipe after the Laplace transform: The wall velocity has a radial direction and is non-zero when the pressure changes (pipe is expanding or contracting).
When the Laplace transform is applied to Equation (3), it is possible to obtain a continuity equation in the form of Equation (13). Relationship (12) is substituted for velocity in the integral over the surface P.
When the complex speed of sound Equation (14) is introduced, one obtains the final shape of the continuity equation for the liquid in a pipe Equation (15).

Momentum Equation
The momentum equation expresses the equilibrium of forces that originate from inertia, pressure and viscosity. It is possible to write it as follows: where the stress tensor is a function of shear viscosity and second viscosity: Substituting Equation (17) into Equation (16) yields: where the convective acceleration is neglected because the liquid velocity is much lower than the speed of sound. When Equation (18) is integrated over the volume defined in Figure 1, the momentum equation for the axial direction is obtained: Here, the function W represents energy loss due to an unsteady velocity profile. The function is well described in [22,23]. The unsteady velocity profile has a significant impact on energy loss in the case of viscous liquids (e.g., oil).
Applying the Laplace transform to Equation (19) and using Equation (15) yields the final shape of the momentum equation:

Transfer Matrix
Transfer matrices determine the transfer of the state vector from the coordinate x = 0 to any coordinate x on the pipe axis. Equations (15) and (20) can be written in matrix form: where: J 0 and J 1 are Bessel functions. The coefficient B includes the Laplace transform of the function W: Equation (21) can be rearranged: Now, the Laplace transform with respect to x is applied to Equation (26), and after rearranging and applying the inverse Laplace transform, Equation (27) is obtained.
where s * is the parameter of the Laplace transform with respect to x. The transfer matrix P appears after computing the inverse matrix and executing the inverse Laplace transform in Equation (27): Thus, Equation (27) can be written in the final form: whereū represents the state vector:ū When the pressure and flow rate at the beginning of the pipe are known, it is possible to compute the pressure and flow rate in any place along the pipe using Equation (29). All information about the pipe and liquid is included in the transfer matrix P f .

Transfer Matrix of the Pipe
Oscillation of the fluid in the pipe is sufficiently described by Equation (29) only for the case when the ends (flanges) of the pipe do not move in the axial direction. The experiment, which is described in the next section, does not fulfil this condition because the excitation is performed only via axial displacement of the flange. Thus, it is also necessary to derive the transfer matrix of the pipe. Generally, the procedure is the same as the abovementioned procedure for the fluid. The transfer matrix of the pipe includes the material properties and mass at both ends of the pipe. The mass represents flanges, valves and sensors because they influence the dynamic behavior of the pipe the most. The scheme is shown in Figure 3.  The transfer equation between places 0 and 3, when there is no liquid in the pipe, is: where the index refers to the corresponding place in Figure 3. Matrices P m and P p are transfer matrices of the mass at the end of pipe and the pipe, respectively. Their derivation follows. The transfer matrix P m comes from the equation of motion: When the flange is ideally stiff, it is possible to consider the same displacement of both flange faces u 0 = u 1 . When the Laplace transform is applied on Equation (32) , the transfer matrix of the flange can be written in the form: The same procedure can be applied on the flange indexed 2 and 3. The transfer matrix of the pipe P P comes from a longitudinal oscillation: Equations (4) to (6) describing the Kelvin-Voigt model of the body can be implied into the Equation (34). Then it is possible after the Laplace transform to write the Equation (35) as: where Equations (36) and (37) can be obtained from the stressσ (0) and the displacementū (0) taken for the pipe beginning. Knowing that ε = ∂u ∂x and σ = F S p . Speed of sound a 0 in the pipe is a complex number: The transfer matrix of the pipe comes from the combination of Equations (35) and (7): where: In this case, the variable x means the length of the pipe. Boundary conditionsF 0 = 0 andF 3 = 0 can be implied to Equation (31), then it is possible to write: When the determinant of the matrix in Equation (41) equals zero, it is possible to find constants E and b for the measured eigenvalues s.

Transfer Matrix of the Entire System
The transfer matrix, which is derived in Section 2.3, allows computation of the fluid dynamics. The transfer matrix, which is derived in Section 2.4, allows computation of the pipe dynamics (in the axial direction) but without fluid inside. Now, the last step is to convert these two matrices into one transfer matrix that describes the pipe + fluid system. Therefore, the following boundary conditions at the ends of the pipe are used: fluid velocity equals flange velocity, and force, which acts on the flange, equals the force that acts on the fluid.
After the Laplace transform:ū Now, Equations (29), (31), (44) and (45) give: where the matrix Λ originates from Equations (44) and (45): When the eigenvalue from the experiment is substituted into the matrix in Equation (46), it is possible to find fluid properties (speed of sound and second viscosity) from the zero determinant of this matrix.

Experiment
The theoretical analysis described in the previous chapter was used to evaluate the experimental data. The experiment was performed on a simple pipe and gives data that are necessary to calibrate the numerical model. The scheme of the experiment is shown in Figure 4. Steel pipes (0.8, 1 and 2 m long) were used. The length is important because it determines the eigenfrequency of the response. The pipe was hung on ropes to ensure that the system mounting had a very low eigenfrequency and did not affect the measurement.
The pipe was filled with water, which was delivered by a pump through a hose and a valve at one end. The air can exit through the valve at the other end. As soon as the pipe was full of water, the valves were shut, and the pipe was ready to be excited with a hammer. The static pressure in the pipe was kept at 1.9 MPa to suppress the influence of air because there was always some amount of air left.
Response to the excitement was measured using pressure transducers, which were mounted to the pipe wall at both ends. The sampling frequency of the transducers was set to the value of ten times higher than the expected frequency of the response (e.g., 26 kHz for the pipe 1 m long), and the uncertainty was 0.25% from the range, which was 10 MPa (one sensor) and 25 MPa (the other sensor). In addition there were accelerometers on the flanges for measuring the acceleration of both ends in the axial direction (range ± 500 m s −2 ).
The response was measured with water and without water. Thus, it was possible to determine four material constants, which are necessary for calibrating the numerical model. First, the modulus of elasticity E (see Figure 2) and damping b were determined from the measurement without water. The damping is supposed to follow a hyperbolic law: where ζ is the damping ratio. The shear viscosity was 0.001 Pa s. Then, the speed of sound and the second coefficient of viscosity were determined from the measurement with water. In both cases, the coefficients were computed numerically to make the determinant of the appropriate matrix equal zero in Equations (41) and (46). The eigenvalue of the system can be computed from a discrete Fourier transform of the measured response according to Figure 5 and Equation (49).

Amplitude [Pa]
Frequency [Hz] f p/2 p w Figure 5. Determination of the real part of the eigenvalue from the discrete Fourier transform.
where the frequency f corresponds to the first eigen frequency of the system.

Results
The experiment without water in the pipe was evaluated based on the records of accelerometers. The influence of air was neglected. Thus, it was possible to determine material properties of the pipes. The values of damping ratio ζ and modulus of elasticity E (see Figure 2) were used for evaluating experiments with water in the pipe. The results are listed in Table 1. The speed of sound is different for different pipe dimensions (diameter and wall thickness) and does not depend on the frequency. The second viscosity is shown in Figure 6 as a function of frequency. When the second viscosity is supposed to be zero, the shear viscosity has to be changed to obtain the correct eigenvalue. Figure 8 shows the sensitivity on the shear viscosity. The damping corresponds to the experiment when the shear viscosity is 0.19 Pa s. This value is almost 200× greater than the physically correct number for water at given conditions (0.001 Pa s).

Conclusions
The mathematical model of pulsating flow in a pipe is derived in the paper. The obtained description in the frequency domain can be useful for stability analysis of a hydraulic system, typically, penstock of a hydro power plant. Beside the stability analysis it is possible to find the liquid and pipeline material properties by using measured eigenvalues of the system. Young's modulus and damping of the pipe as well as second viscosity and speed of sound of water were evaluated. Obtained result comes from the simple experiment with the system that is excited by an external force (hammer or water hammer in the case of very large pipelines).
Despite the inclusion of elastic behavior of the pipe, it is still possible to extend the model by including the influence of the surrounding air, which is important in the case of measuring systems transporting gas. The reader should also take into account that the solution is obtained in the frequency domain, and it is not easy to use the results in the time domain.