Analyses of 2D transient models for vaporous cavitating flows in reservoir-pipeline-valve systems

Water hammer causes not only pressure rise but also pressure drop in pipes. When the pressure drops below the vapor pressure of pure liquid, cavitation will occur. The occurrence, development and collapse of the cavitation bubble are processes of phase changes, which have the characteristics of complex physical mechanisms and difficult cognition. Therefore, it is important to analyze the problem of transient vaporous cavitation.

A lot of experiments and theoretical research have been carried out and many mathematical models of transient cavitating flows have been proposed in recent years. The Discrete Vapor Cavity Model (DVCM) with constant wave-speed proposed by Streeter (1983) assumed that the whole cross section of the pipe is full of vapor when the pressure of the section is lower than the vapor pressure and its pressure holds at vapor pressure until the collapse of the vapor. The evolution of the cavitation cavity is calculated by mass conservation of the liquid column. The model shows the basic characteristic of the transient cavitating flow and is easily programmed, so it is widely applied in the existing commercial water hammer software (Flowmaster, Hammer, etc., Marcinkiewicz et al. 2008). However, the DVCM may predict unrealistic pressure peaks in the case of collapse of multi-cavities (Simpson & Bergant 1995). The Discrete Gas Cavity Model (DGCM) (Simpson & Bergant 1994), which is similar to the DVCM model, has the hypothesis that there is always a bubble (i.e. gas nucleus) at the computational node and the size of the bubble changes with the change of ideal gas (isothermal or adiabatic). The calculated results predicted by the DGCM model are in good agreement with the experimental results and less affected by the input parameters than those of the DVCM model. The Generalized Interface Vaporous Cavitation Model (GIVCM) carried out by Simpson (1986) included water hammer equations for the liquid phase, two-phase equations for a distributed vaporous cavitation region, shock equations for condensation of the liquid-vapor mixture back to liquid, and equations for a discrete vapor cavity. Considering the two-phase vapor-liquid mixture as a single-phase mixture with average physical characteristics, the model accords better with the physical process but the algorithm is complex and time-consuming. In order to avoid the unrealistic pressure peaks, Adamkowski & Lewandowski (2012) proposed a new DVCM model, which assumes that the cavity volumes calculated in several computational nodes are shifted to the main cross section without breaking the conservation law. Taking the cavitation transport model into consideration, Sumam et al. (2010) adopted the MacCormack numerical scheme and a new pressure relaxation method to overcome the numerical instability caused in the phase transition process. The results predicted by the two-phase homogeneous equilibrium vaporous cavitation model with frequency dependent friction proposed by Shu (2003) were compared with those of the DVCM model, and the effectiveness of the new model was tested and verified. Great achievements have been obtained for transient cavitating flows, but there are significant differences in pressure amplitudes and phases (Szymkiewicz 2002) between the simulated and experimental results. This is related to the steady friction term or quasi-steady friction term used in these models, which cannot describe the complex characteristics of the shear stress and underestimate the energy dissipation observed. To solve such problems, a series of one-dimensional (1D) unsteady friction models are developed from the point of unsteady wall shear stresses. For a fast closure of the downstream-end valve in a simple reservoir-pipeline-valve (RPV) system, Urbanowicz et al. (2012) came to the conclusion that the combination of the 1D unsteady convolution-based (CB) (Zielke 1968) and DGCM model brought the results much closer to the experimental results. Bergant et al. (2005) introduced the Momentum Correction Factor (MCF) in the momentum equation and verified that the results of the DGCM-CB-MCF model were more consistent with the experiment. Storli & Nielsen (2011) found that the calculated results by using the Brunone unsteady friction model with a constant empirical coefficient k matched only at specific nodes and suggested that a position-dependent and time-dependent k should be employed. As another efficient method, the more accurate two-dimensional (2D) unsteady friction model which takes the velocity distribution of the cross section into consideration attracted more and more attention. Vardy & Hwang (1991) established the quasi 2D five-region model for both laminar and turbulent flows and verified the effectiveness of CB expressions in transient laminar flows. The numerical method is accurate and numerically stable, but it has the disadvantage of low efficiency in solving the sparse matrix in the coupling linear equation. For the same space and temporal discretization, the different numerical schemes proposed by Zhao & Ghidaoui (2003) and Korbar et al. (2014) had the same precision, but with less computation time. To reduce the computational cost, Ohmi et al. (1985) obtained the pressure and average velocity from the 1D continuity equation and the velocity distribution was derived from the 2D momentum equation. A similar procedure was used by Jang et al. (2016). A mixing length algebraic turbulence model was employed in smooth pipes (Silva-Araya 1997) and rough pipes (Silva-Araya & Chaudhry 2001). The two-region turbulence model proposed by Pezzinga (1999) was applied to simulate the transient non-cavitating flows in smooth and rough pipe networks. Zhao & Ghidaoui (2006) developed the quasi 2D model and pointed out that this model was not suitable for transient flow with high initial Reynolds number. Riasi et al. (2013) put forward the quasi 2D model and carried out the dimensionless analysis. The Baldwin-Lomax turbulence model, which was performed in both quasi-steady and frozen forms, was investigated by Wahba (2009) in transient non-cavitating flows. Pezzinga & Santoro (2017) studied two different cavitation models, both in 1D and 2D forms. Moreover, the 2D computational fluid dynamics (CFD) method was adopted by Wang et al. (2016) to analyze the vapor distribution, pressure and velocity fields in pipelines and the results demonstrated that the method could simulate the pressure variations and visualize the associated physical processes. Based on the full cavitation model, the commercial code CFD-ACE+ is used by Pinho et al. (2014) to investigate the corresponding characteristics of the transient cavitating events.

Based on the above facts, further research is carried out in the present work for transient cavitating flows based on the 2D model proposed by Pezzinga & Cannizzaro (2014). The comparisons between the 1D model and 2D model are related to the classical reservoir-pipe-valve system. The effectiveness and superiority of the 2D model is comprehensively tested and verified against a series of experimental data. In order to judge the reliability of the turbulence models and their assumptions in the transient cavitating 2D model, a further sensitivity investigation of different eddy viscosity models (five-region turbulence model and two-region turbulence model) on the accuracy is carried out. Moreover, both the quasi-steady and frozen forms for implementing the five-region turbulence model are examined in order to assess their range of applicability related to the non-dimensional parameter P.

For the 2D model, the pipe is discretized into cylinders with constant area in the radial direction and reaches with constant length in the longitudinal direction. The 1D model is a special case of the 2D model ⁠. The adopted time step is ⁠. For each discrete cylinder, axial velocity u is defined at each radial grid center; shear stress is defined on the inside and outside of the grid surface, as shown in Figure 1.

Indices ‘i’ and ‘j’ refer, respectively, to node numbers indicating longitudinal direction x and radial direction r; index ‘t’ refers to time t; indices ‘pre’ and ‘cor’ refer to the predictor and corrector steps. The values at time t are calculated as the average of the predictor and corrector values.

The numerical simulations are concerned with the fast valve closure in the reservoir-pipe-valve system. At the upstream point, the value of the pressure head is given while the local velocity is obtained by time-marching the momentum equation, and the x-derivative of the head is evaluated forward both in the predictor and the corrector step. At the downstream point, the pressure head during the valve closure is interpolated from a set of tabular experimental results at regular increments of time and evaluated by time-marching the continuity equation, and the x-derivative of the mean velocity is evaluated backward both in the predictor and the corrector step. When the valve is completely closed, the mean velocity is zero and the pressure head is evaluated backward both in the predictor and the corrector step. Moreover, a symmetry boundary condition is applied to the local velocity at the pipeline centerline and a no-slip condition is enforced at the pipe wall.

In order to validate the proposed mathematical model and numerical scheme, the numerical results were compared with experimental tests carried out by Simpson (1986). Layout of the laboratory apparatus is shown in Figure 2. The laboratory apparatus is composed of a smooth copper straight sloping pipe (internal diameter 19.05 mm, length 36 m and thickness 1.6 mm) connecting two pressurized tanks. A specified pressure in each of the tanks is controlled by a computerized pressure control system. The ends of the pipe have an elevation difference of 1.0 m. Three piezoelectric transducers are located at the downstream valve, at the upstream 1/4 point (9.0 m from the upstream end), and at the downstream 1/4 point (9.0 m from the downstream valve), respectively. The wave speed is 1,280 m/s.

Nine different experimental tests were studied and a summary of all corresponding parameters considered is listed in Table 1. All simulated parameters are consistent with the experimental ones.

When several cylinders exist, there are many methods to determine the radial thickness. Bratland (1986) implied that the equal area can give more detail close to the pipe wall than the equal thickness ⁠. Since the radial disturbances are caused primarily by the no-slip condition at the wall and since diffusion processes are much slower than wave processes in low Mach number flows, the equal area has a significant advantage. So all results presented herein are obtained with cylinders of equal area. In addition to the selection of the thicknesses of various cylinders, the total number of them should also be selected. First of all, the grid independence verification of the 1D model is carried out for each test to determine Nx. Eight different numbers of pipe reaches (Nx = 8, 16, 32, 64, 128, 256, 512 and 1024, respectively) were considered and the parameter of the maximum head at the valve was used to assess the performance of the model. The results for R2 are shown in Figure 3(a). The maximum head at the valve increased gradually with the increase of Nx. When Nx > 32, the maximum head at the valve almost remained unchanged, therefore the optimal value of Nx is 32. Then, the grid independence verification of the 2D model is carried out to determine Nx and Nr. Four different numbers of discrete cylinders (Nr = 21, 30, 50 and 100, respectively) at a range of different Nx were considered and the parameter of the maximum head at the valve was used to assess to performance of the model. The results for R2 are shown in Figure 3(b). For each specific Nr value, the maximum head at the valve almost remained unchanged when Nx > 256. For the particular problem considered here, the use of significantly fewer cylinders Nr = 50 can not only save computer resources, but also guarantee the stability of the results. Therefore the value of Nx is 256 and the value of Nr is 50. So for the cavitating flow, the optimum ratio of the number of reaches (Nx) and the number of discrete cylinders (Nr) is found to be about 5.0.

Preliminarily, the mathematical model and numerical scheme were tested for a non-cavitating test (R1), in which the transient pressures stay above the vapor pressure of the liquid. The pressure curves after the rapid closing of the valve at each pipeline monitoring point calculated by the 1D and 2D models are shown in Figure 4. It can be observed from the figures that there exists a small attenuation of the maximum head predicted by the 1D model. The first four peak pressures at the valve are 55.48 m, 55.15 m, 54.82 m and 54.50 m, respectively. And the first four peak pressures at the valve for the 2D model are 55.36 m, 54.33 m, 52.91 m and 51.46 m, respectively. Based on the above-mentioned quantitative analyses, the results predicted by the 2D model are closer to the experimental results because the velocity profile in the cross section is taken into account in the 2D model, which allows the wall shear stress to be more accurately evaluated. It can be noted that, also for upstream 1/4 point and downstream 1/4 point, the improvements made by the 2D model are similar to those previously shown.

The measured head oscillations compared with the results of the 2D model and 1D model are shown in Figure 5 for R2 at all monitoring points. The phenomenon of this test is characterized initially by a Joukowsky pressure rise caused by the rapid closure of the valve, and then the pressure reduced to vapor pressure and the subsequent pressure peak was caused by the collapse of the cavity. It can be noted that the 2D model is much better than the 1D model in reproducing the experimental results at all monitoring points. Compared to the experiment, there is a phase lag phenomenon of the pressure predicted by the 1D model and the discrepancies are magnified with the increase of time. The dissipation caused by friction could be responsible for the shaping and shifting of the pressure trace predicted by the 1D model. On the contrary, the 2D model reproduces the experimental oscillations with very good accuracy throughout the experiment. From the quantitative point of view, both models can predict that the pressure peak caused by the collapse of the first cavity at the valve, ⁠, is greater than the first peak pressure due to Joukowsky pressure rise. is 91.00 m in the experiment, 104.09 m predicted by the 2D model and 106.34 m by the 1D model.

An analogous comparison is shown in Figure 6 for R3. Also, in this case both the two models can reproduce well the above-mentioned phenomenon and predict that the pressure peak caused by the collapse of the first cavity at the valve, ⁠, is greater than the first peak pressure due to Joukowsky pressure rise. is 95.15 m in the experiment, 101.07 m predicted by the 2D model and 102.01 m by the 1D model. For the phase of the pressure trace, there is a good agreement between the calculated pressure curves of the 1D model and the experiment. But the 2D model has a tendency to anticipate the head oscillations in advance with respect to the experimental ones, which had also been observed and attributed to the neglecting of the convective terms and the release of dissolved gas in Pezzinga & Cannizzaro (2014).

Similar comparisons for the rest of the tests were also carried out and analyzed. All the comparisons show that the 2D model has significant improvements in both magnitude and phase shift of pressure head traces compared with the 1D model in most cases. Due to the different evaluation of dissipation, the 1D model allows a good reproduction of the instant of the pressure peaks only for the first and second oscillations, whereas the results of the 2D model are in good agreement with the experimental results in most cases. The oscillation phases computed by the 1D model are more consistent with the experimental ones, only for R3.

A summary of the comparisons between the results of the numerical models and the experimental data is shown in the bar diagrams of Figure 7, where the relative error (in percent) of the 1D model and the 2D model on maximum heads at the valve are reported. It can be observed clearly that the relative errors of the 2D model are generally smaller than those of the 1D model except for R7 and R8. The 2D model is always in better agreement with the experimental results than the 1D model because the 2D model takes the velocity profile into account. However, the reasons for the poor performance of the 2D model for R7 and R8 are worth studying in depth. Various simplifying hypotheses in the 2D model can affect the predicted results. Firstly, the dissipation caused by friction is not the only factor responsible for the shaping and shifting of the pressure trace in the transient flow. The effects of inertia caused by the unsteady velocity profiles, which is not considered, may be far greater than the wall shear stress effect on the pressure trace (Storli & Nielsen 2013). Then, according to Simpson (1986), the vapor fraction at the valve and the time of vapor existence increase with the increase of the initial velocity. The hypothesis of axial-symmetry is likely to lose efficacy in the presence of conspicuous vapor cavities, which would lead to the above-mentioned large relative errors. More detailed experimental results and more advanced models are needed to deepen these aspects.

A comprehensive study of the behavior of axial velocity profiles in the transient non-cavitating flow was numerically investigated by Riasi et al. (2009). A further investigation into the behavior of the velocity profiles in the transient cavitating flow was carried out in the paper. The predicted profiles of velocity divided by the initial mean velocity at the upstream 1/4 point at different instants by the 2D model for R2 (P = 22.99) and R9 (P = 7.94) are shown in Figure 8. When the positive pressure wave passes from the upstream 1/4 point at time t1, there is a uniform shift in the velocity profile in the core region and a strong gradient and a vortex ring (reverse flow) in the near-wall region because of the no-slip condition at the pipe wall, which produces a great wall shear stress. The generated vortex ring is subsequently diffused toward the pipe center and the reverse velocity moves to the core region. When the reflected positive pressure wave passes from the upstream 1/4 point at time t2, the whole cross section has a negative flow, which presents an S shape distribution and another vortex is generated at the pipe wall. The same condition occurs when the negative pressure wave comes from the valve at time t3 and also when the reflected positive wave passes from the upstream 1/4 point at time t4. The process of velocity shift, vortex generation and vortex diffusion will continue to occur as long as there are pressure waves in the pipe. It can be noted that the velocity profiles are different from the corresponding ones in the 1D condition with equal discharge and are very similar to those in the transient non-cavitating flow.

Turbulence models are needed to calculate the shear stress related to energy dissipation in the 2D model. Algebraic turbulence models, especially the two-region and the five-region eddy viscosity models, which relate the eddy viscosity to the mean velocity through some algebraic functions, have been widely used. The five-region model divides the flow cross-sectional region into a viscous sublayer, a buffer zone, and a core region. The buffer zone and core region are further sub-divided to reflect different viscosity characteristics. The two-region model divides the flow cross-sectional region into a viscous sublayer and a turbulent region. Ghidaoui et al. (2002) investigated the two-region and the five-region eddy viscosity models regarding the accuracy of the turbulent non-cavitating transient flows. And further systematic studies to evaluate the effects of these assumptions on cavitating transient flows were carried out in the paper.

Figure 9 depicts the computed head trace results obtained from the two-region and the five-region models and the measured head trace results for R2. Obviously, the pressure peaks as well as their shapes predicted by the two models almost coincide and they are in good agreement with the experimental results. The conclusions of the rest of the tests are consistent with those of R2, which shows that the turbulence models are insensitive to the form of the spatial distribution of the eddy viscosity within the core region. This is because the turbulent stresses are negligible compared to both the inertia and the pressure gradient in the core region. In future investigations, it is suggested that more efforts should be concentrated on the establishment of accurate models for the near-wall region.

The five-region turbulence model is developed from steady flows and assumptions are needed when applying the model to transient flows. For the quasi-steady form, the turbulence model is updated at each time step, providing values for the instantaneous eddy viscosity as a function of the instantaneous velocity profile. For the frozen form, the eddy viscosity values at any spatial location and any time step are kept frozen and equal to their initial steady state values at the same spatial location. The influence of the quasi-steady form on the accuracy of the transient non-cavitating flows was studied by Ghidaoui et al. (2002). In order to evaluate the influence of the quasi-steady and frozen forms regarding transient cavitating flows, a further intensive research of both forms on cavitating flows was carried out in the paper.

Numerical simulations for both forms of five-region turbulence model are performed first for R2 (P = 22.99). The predicted pressure plots at the monitoring points are compared with the corresponding experimental results in Figure 10. The results obtained by both forms are almost on top of each other and in good agreement with the experimental ones. This is because the shear pulse is mainly within the viscous sublayer and has no influence on the turbulence intensity and structure (He & Jackson 2000). That is to say, the flow characteristics are essentially identical to their steady state characteristics after the wave passes the particular location. Therefore, the frozen eddy form and the quasi-steady form are equally acceptable. Greenblatt & Moss (1999) obtained a similar conclusion for the temporally accelerating flow.

Analogous representations are performed in Figure 11 for R9 (P = 7.94). In this case, the results obtained by both forms are close to each other and have a slight difference from the experimental ones. It seems that the two forms work much better for R2 than R9 and this is partly due to the fact that P for the former test is greater than that for the latter test. The conclusion in the cavitating flows is similar to that summarized by Ghidaoui et al. (2002) in non-cavitating flows. This is because when the value of P is of the order of 1, the wave-induced shear pulse diffuses out of the viscous sublayer into the buffer zone and core region. As a result, the velocity profile does not present a logarithmic-like shape and the turbulent production and dissipation are different from their initial steady state characteristics in the region outside the viscous sublayer. Therefore, the structure and intensity of turbulence experience significant changes during the transient event, which makes the frozen and quasi-steady form questionable. Unfortunately, there are no available data that could be used to judge on the better form in such a complex phenomenon. Therefore, there is an urgent need for a systematic experimental study so as to determine the more accurate forms.

Based on the governing equations of the vapor-liquid homogeneous mixture, 1D and 2D models for transient cavitating flow were constructed, respectively. By comparing the experimental pressure curves of the quick closing of the downstream valve in a simple RPV system, the main findings of the paper are summarized below:

• (1)

  For the transient cavitating flow considered here, the optimal ratio of the number of grid reaches (Nx) and the number of discrete cylinders (Nr) is found to be about 5.0 in the 2D model.

• (2)

  The results predicted by the 2D model are in better agreement with the experimental pressure head traces including magnitude and phase than those predicted by the 1D model in most tests. Compared with the experimental results, the pressure curve predicted by the 1D model under most conditions has a phase lag phenomenon and the discrepancies of phase differences are magnified with the increase of time.

• (3)

  The relative errors on maximum heads predicted by the 2D model are generally smaller than those of the 1D model. In cases of conspicuous vapor cavities, the relative errors of the 2D model are relatively larger.

• (4)

  The pressure curves predicted by five-region and two-region turbulence models are similar, which indicates that the transient cavitating flows are insensitive to the magnitude and distribution of the eddy viscosity in the core region. Therefore, future efforts in turbulence modeling of transient flow should be concentrated on the establishment of accurate models for the near-wall region.

• (5)

  If the non-dimensional parameter P is much greater than 1, the turbulence distribution rapidly changes and the velocity profile has enough time to form the fully-developed velocity. So the quasi-steady and frozen forms are both acceptable. If the non-dimensional parameter P is of the order of 1, the quasi-steady and frozen forms become less accurate.

This work was supported by the National Natural Science Foundation of China (Grant Nos 51779257, 51479196, 51179192), and the Program for New Century Excellent Talents in University (NCET) (Grant No. NETC-10-0784).