Simulation of unsteady flow with cavitation in plastic pipes using the discrete bubble cavity and Adamkowski models

The work presents two modified cavitation models for the simulation of transient flow in pressure plastic pipes. The first model is a discrete bubble cavitation model, the prototype of which was presented by Shu, and the second one is the Adamkowski model. In the latter model, the problem encountered in the classical model (DVCM -discrete vapour cavitation model) related to artificial damping of pulsation, which results from approximate timing of cavity opening and collapse. In both models, the corrected efficient calculation of Zielke convolution integral was used to simulate the unsteady wall shear stresses. The numerical results from the two models agree well with the results of measurement in the literature.


Introduction
Transient cavitating flow often occurs in hydraulic piping systems that were not thoroughly analysed at the design stage. This is often the result of abandoning the detailed numerical analysis of possible consequences resulting from the dynamics of flows. Designers focus mostly on stationary flows and the only dynamic equation applied is the simplified Joukowsky formula for the determination of appropriate wall thickness. In fact, not only high pressures but also low pressures resulting from water hammer are a threat. In conditions of vapour pressure, the liquid stream column separation occurs. Two types of transient cavitation can be distinguished in pressure pipelines [1]. The first one is gaseous cavitation, which is a slow-changing phenomenon, the air bubbles emitted from the gas are dissolved in the liquid as a result of the pressure drop (the secretion lasts for a second, and the possible dissolution takes even longer). Gaseous cavitation results in an average reduction of the pressure wave speed, therefore, it helps increase the damping of pressure waves. The above features are conducive to securing the system, hence, this type of cavitation is beneficial even if it does not occur in a steady state in the analysed system. The second type of cavitation, the vaporous cavitation, takes place when the pressure drops to the vapour pressure. This is a fast-changing phenomenon, i.e. the evaporation requires only microseconds to form and to disappear given that there is a sudden increase in pressure. This type of cavitation is responsible for the appearance (in analysed hydraulic systems) of secondary pressure waves, which result from the implosion of vapour areas created in the pipe. In contrast, the interference between the primary (resulting from a fast closure of the valve) and these secondary pressure waves can result in pressure surges exceeding even twice the maximal pressure calculated from the Joukowsky formula. In practice, strengthening and wave suppression depend on the geometry and material properties of the analysed hydraulic systems, initial conditions prevailing in the system before the occurrence of unsteadiness as well as boundary conditions taking place after the occurrence of unsteady flow. When implosion of vapour bubbles occurs at the wall of the pipe, cavitation erosion occurs: when the bubbles collapse (the vapour areas decrease until they completely disappear) the local velocity of the separating phase edges increases rapidly. The collision of the stream surrounding these edges in the last stage of closure is a source of large local pressures, which may be much higher than the pressure calculated from the Joukowsky formula. Due to the above-mentioned impacts of cavitation on hydraulic systems, there is a great need to analyse them at the design stage. In this paper, a one-dimensional flow modelling was undertaken in a simplified hydraulic system consisting of a pressure vessel, a pipeline and a valve. From the literature on unsteady flows occurring in elastic pipes [2][3][4][5][6][7][8], a number of models can be distinguished, including discrete vapour cavity model (DVCM), discrete gas cavity model (DGCM), discrete Adamkowski cavity model (DACM) and discrete bubble cavity model (DBCM). DVCM [3] and DGCM [2,8] have already been modified and used to model transient cavitating flows occurring in plastic pipes [9][10][11][12]. In this paper, the modified DACM and DBCM models will be developed and presented. Unsteady hydraulic resistance (skin friction loss) will be modelled in accordance with the improved solution presented in [13,14]. The calculated results obtained from the modified DBCM and DACM models will be compared with the experimental results found in the literature.

Modelling transient cavitating pipe flow
The unsteady flow of liquids in plastic pipes can be described by the system of equations [15]: Comparing the above system of equations and the system describing the flow in metal conduits [14], the only difference between the systems appears in the equation of continuity (the first of the equations (1)) with an additional convolution integral term describing a retarded strain rate in the plastic pipe wall. The solution of the above system for an idealised case of flow in which the gas areas were not able to form can be found in the literature [15,16]. In this paper, we extend to plastic pipes the following two models: Shu's bubble cavitation model (discrete bubble cavitation model -DBCM) [4], and the Adamkowski model [17,18] (discrete Adamkowski cavitation model -DACM). TDACM is based on the above system of equations (1) assuming only a specific solution defining the average value of flow velocity at the nodes in which cavitation occurs. The DBCM needs a more detailed explanation.

Discrete bubble cavity model (DBCM)
The equation of motion of two-phase vaporous cavitating flow in the horizontal pipeline is: The equation that describes the mixture density is: Writing the continuity equation separately for the vaporous and liquid phase whereaverage velocity of vapour flow [m/s] andaverage velocity of liquid flow [m/s]. By adding the above equations to each other, one gets: Next assuming the non-slip flow, in which the vapour phase has the same velocity as the liquid one = = : After differentiation: The total derivative of equation (3) is: Using equations (11) and (12) in equation (9.1) yields: The first term on the left-hand side of the above equation, in square brackets, is a formula for the pressure wave speed. Due to the fact that in this paper the numerical solution will be based on a constant rectangular grid, the value of the pressure wave propagation velocity will be determined for the flow occurring in the initial time. Subsequently, we deal with a single-phase liquid flow, which determines Therefore, the equation of continuity can be finally described as: In non-slip flows, the share of the dispersed phase is of the statistical nature, i.e. the volume and mass concentration are equal to the relevant proportions of dynamic nature, i.e. transport concentration and the degree of dryness [20]; then the following relationship applies: = , where is superficial velocity of the liquid phase.

Discrete Adamkowski Cavity Model (DACM)
The main assumption of the Adamkowski model [17,18] is that a segment with few vapour zones may be replaced with only one vaporous zone, which is equivalent with regard to the mass and energy conservation laws. When the pressure falls to a vapour pressure in a given cross-section, the conservation law requires the equality between the cavity volume in the corresponding pipe segment and the sum of the cavity volumes in the original pipe segment. This law also requires the wave period along the pipe length to be preserved. In the intermediate cross-sections " " of the pipeline, two streams of the flowing liquid with velocities + and − are replaced with one continuous fluid stream, represented with the velocity . The energy conservation law requires that the sum of the kinetic and potential energies of elementary liquid masses separated by a vaporous zone has to be equal to the total energy following the reconnection of the two separate liquid streams, which is why: In the case of the horizontal pipeline, the difference ( + − − ) is equal to 0. The + and − are altitudes of the centres of cross-sectional surfaces limiting the cavity within the computational cross-section "j" from the right and left sides, respectively. As in the water hammer case, the mean flow direction changes numerous times until final flow suppression occurs, the Adamkowski equation for mean velocity is expressed by: This equation (18) is used when the pressure falls to the vapour pressure = , i.e. when the column separation of liquid stream takes place.

Wall shear stress model
In DACM the wall shear stress is modelled as for a single-phase flow. The equation used in this work which represents wall shear stress as a sum of quasi-steady component and an unsteady component is: and its effective numerical solution, employed in this work is [14]: , , are constant during simulation and are calculated prior to the simulation with the coefficient describing the effective weighting function; Δ̂is a dimensionless time; . ( ) and . ( ) are respectively the classic inefficient weighting function (laminar or turbulent depending on the type of flow before the transients) and efficient weighting function (are the approximation of the classic weighting function). In this paper, the three-term efficient weighting function, presented in a former study [13] was used. These coefficients as the simulation was in the turbulent regime were properly scaled using the universal weighting function procedure. The second analysed DBCM is a two-phase flow model, which required the following modification in the effective numerical solution: The is the mixture viscosity calculated in present work from a simple Dukler's solution [4]:

Numerical scheme
With the formulation of cavitation and shear stress model above, the method of characteristics (MOC) based on a rectangular grid is adopted in this study to solve the continuity and motion equations (e.g. equations (10) and (15) when the first cavitation model is used). The principle of this MOC is to convert the original set of partial differential equations (PDE) into an ordinary differential equation (ODE) by considering the wave evolution process along the characteristic lines ( ± 0 ), at which the solution can be integrated from known information at last-time step. The obtained ODE is then solved by the central finite difference (FD) scheme to determine the numerical solutions at the current time step along the pipeline. To ensure the stability and convergence, a Courant number for the simulation process is less than or equal to 1 for the reservoir-pipeline-valve system used in this study. More details about the mathematical formulation and application of MOC can be found in a former work [16]. The final computer programmes were written in MatLab, which is a multi-paradigm numerical computing environment and proprietary programming language developed by MathWorks (plain MatLab, no extra toolboxes were needed). The identification of the calculation algorithm forced to filter the creep compliance function according to the procedure presented in former studies [15,16].

Simulation examples
The obtained mathematical formulas described in detail in the previous section were used to write computer programmes for modelling unsteady runs in a hydraulic system built from a reservoir-pipevalve. The number of experimental studies on unsteady flows in this type of systems during which cavitation areas were created is limited, with only Güney [9,10] and Soares [12] being known to the authors. In this paper, the studies obtained by Güney will be used for comparison. The test stand was composed of a 43.1 meters long (R=0.0208 m and e= 0.0042 m) LDPE pipe. Two cases of turbulent flow were analysed, in which the main difference was a significant difference in the temperature of the flowing liquidin Case 1 the temperature was 1 = 13.8℃ and in Case 2 -2 = 38.5℃. The temperature significantly affects the mechanical properties of the pipe material: a drastic change in the Values of the necessary input parameters for computer programmes are summarised in Table 1, while coefficients describing the creep function for the two analysed temperatures are given in Table 2. The method of characteristics was used with a constant number of reaches N=64. Such a selected number of reaches the following time steps ∆ 1 = 0.0022 [ ] and ∆ 2 = 0.0031 [ ]. The results of the simulation tests obtained in relation to the experimental data are presented in Figure 1.   From the presented results, it can be seen that in the analysed Case 1 (Figures 1a and 1b) both flow models produce accurate results. The detailed remarks regarding Case 1 are as follows: a) the visible peak of pressure at the top of the first amplitude ( Figure 1a) as well as the subsequent drop followed by a gentle experimentally observed increase in pressure was modelled with both DBCM and DACM; b) the slight peak observed in experimental studies at the beginning of the second amplitude ( = 0.96 ) was modelled correctly by both models (Figure 1b); c) in the unsteady flow without the cavitational liquid stream column ( > 1.65 ), it can be seen that the phase shift of the simulated amplitudes increases over time. The result of the above behaviour should be sought in the change of the pressure wave speed during the water hammer or measurement and approximation errors at the stage of defining the experimental creep function by Güney. Detailed comments on analysed Case 2 were similar, although it was noted that: a) in the first experimental amplitude, as in Case 1 after the original peak, a gentle pressure increase. Which, remembering given by Güney, the pressure drops on the length of the pipe (Case I: ℎ = 2.88 and Case II: ℎ = 2.93 ) is unusual. For similar initial parameters, one expects a similar course of pressure changes, as modelled by the analysed numerical models ( Figure 1c); b) both numerical models (DBCM and DACM) starting from the third amplitude overestimate the minimum pressures occurring between successive amplitudes (see the results for = 3.1 ); c) in the analysed case, it was harder to find simulation differences between the analysed numerical models (Figure 1d).

Conclusions
The paper presents the modelling of unsteady flows of liquids in pressure polymer pipes. During pressure drop to the saturated vapour pressure, there were areas of cavitation that were modelled with two modified models, namely: DBCMbased on two-phase flow equations and DACMan improved version of the classical DVCM model. The usefulness of the presented models for flow analysis in plastic pipes was examined for the first time. From the research, it becomes clear that the presented models provided high-accuracy simulation of the unsteady flows with the cavitational liquid column separation, and that both similarly modelled cavitation areas. It was hard to see simulation differences between the results obtained from the DBCM and DACM models.
Subsequently, a precise quantitative analysis of the results obtained will be carried out in relation to other experimental results known in the literature in which cavitation areas appeared in plastic pipes. Applied experimental creep functions (developed by Güney, and filtered using method described in a former paper [16]) showed high efficiency in modelling the analysed flows. This further indicates that the exact "maps" of the creep function in a wide range of times and temperatures will allow in the future a complete exclusion of the calibration method (estimation of creep compliance coefficients) commonly used today.