Modeling Transient Pipe Flow in Plastic Pipes with Modified Discrete Bubble Cavitation Model

: Most of today’s water supply systems are based on plastic pipes. They are characterized by the retarded strain (RS) that takes place in the walls of these pipes. The occurrence of RS increases energy losses and leads to a different form of the basic equations describing the transient pipe flow. In this paper, the RS is calculated with the use of convolution integral of the local derivative of pressure and creep function that describes the viscoelastic behavior of the pipe-wall material. The main equations of a discrete bubble cavity model (DBCM) are based on a momentum equation of two-phase vaporous cavitating flow and continuity equations written initially separately for the gas and liquid phase. In transient flows, another important source of pressure damping is skin friction. Accordingly, the wall shear stress model also required necessary modifications. The final partial derivative set of equations was solved with the use of the method of characteristics (MOC), which transforms the original set of partial differential equations (PDE) into a set of ordinary differential equations (ODE). The developed numerical solutions along with the appropriate boundary conditions formed a basis to write a computer program that was used in comparison analysis. The comparisons between computed and measured results showed that the novel modified DBCM predicts pressure and velocity waveforms including cavitation and retarded strain effects with an acceptable accuracy. It was noticed that the influence of unsteady friction on damping of pressure waves was much smaller than the influence of retarded strain.


Gaseous and Vaporous Cavitation
Cavitation is one of the natural phenomena whose thorough understanding should be a scientific priority. Among others, it takes place when gas is released from the liquid. It occurs in hydraulic systems (water supply, hydropower, heating, cooling, etc.) in which the flow (forced by the pressure gradient) takes place through pressurized pipes. There are two types of cavitation: gaseous and vaporous [1][2][3].
More dangerous is vaporous cavitation, which occurs when the pressure drops to the saturated vapor pressure. This type of cavitation is rapidly changing, as it only takes place during the duration of the reduced pressure. In the literature, there is a group of mathematical models based on this type of cavitation, the so-called discrete vapor cavity models (DVCM) [1][2][3]. In the event of a water hammer, the reduction of pressure to the saturated vapor pressure takes a relatively short time, which is followed by an implosion of the resulting vapor regions. The implosion is accompanied by large local increases in the velocity of the liquid, because the cavitation space must rapidly be filled with liquid at the time of pressure increase (above the vapor pressure). The impact of the liquid against the walls of the pipe (as well as the walls of other elements of the systems: valves, turbines, pumps, flow meters, etc.) results in cavitation erosion in the long term. Irreversible losses appear in the material of the walls of the pipes and other elements of the system. Sections in which such erosion takes place are systematically weakened in terms of strength, and it is in these places that leaks or, in extreme cases, complete damage of the structure can occur. Cavitation also leads to a reduction of the efficiency of hydraulic systems, contributing to the deterioration of the operation of energy-saving systems in hydraulic drives [4].
The second type of cavitation, namely gaseous cavitation, is a slowly changing phenomenon occurring in systems with unsteady flows (dynamic, rapid changes of velocity and pressure) or large pressure drops along the length of the system. Each liquid dissolves a certain amount of air (possibly a different gas). In water systems (water supply networks), the average amount of dissolved air is about 2%. In oils, on the other hand, the amount of dissolved air can reach up to about 10%. Hence, the influence of this type of cavitation is much more noticeable in oil-hydraulic systems than in water supply systems. Interestingly, during water hammer, such cavitation areas, due to the large time necessary for desorption and absorption, are beneficial. Their presence causes a faster damping of dynamic waveforms, as the "air bags" emitted by their action resemble local air-liquid shock absorbers. The influence of this type of cavitation is still poorly understood both experimentally and theoretically. There is a group of models called discrete gas cavity models (DGCM) [1,2], which take into account the influence of free gas in a simplified way.
The type of pipe material also significantly affects the intensity and timing of transient phenomena [5]. The flows in metal pipes with vapor-gas cavitation areas are well recognized and described. However, if we look at plastic pipes, which are now starting to displace metal pipes (especially in water supply systems), the researchers have mostly used the two basic cavitation models, i.e., the DVCM and the DGCM. Apart from these two models, alternative models have been developed, including a revised version of the DVCM model proposed by Adamkowski [6,7] as well as a model based on two-phase flow equations that can be called a discrete bubble cavity model (DBCM), which was developed by Shu [8]. Shu's model does not generate the unrealistic pressure spikes due to flow discontinuity at each computational section [8] that have been found in DVCM simulations. In DGCM, it is difficult to assign the physical amount of free air at computational sections along the pipeline. The model that is based on two-phase pipe flow equations is in principle more realistic than the model that is based on single-phase pipe flow equations with cavities lumped at computational sections. However, the aforementioned discrete Adamkowski cavity (DACM) and DBCM models have not been previously used to model transient cavitating flow in plastic pipes. The main objective of this paper is to present a novel DBCM that will enable the simulation of transient cavitating flows in plastic pipes.

Transient Cavitating Flow in Plastic Pipes
Among the phenomena accompanying transient flows, the most important are (1) unsteady friction (UF, other name: skin friction), (2) cavitation (CAV), (3) viscoelastic property of pipe deformation (flow in plastic conduits) associated with the retarded strain (RS), and (4) mutual fluid-structure interaction (FSI) of the flowing liquid with the vibrations of the pipe walls. In this work, we will focus in detail on the first three phenomena. They will be implemented in the revised DBCM model. The continuity equation with the retarded strain term was originally proposed by Rieutord and Blanchard [9]. Cavitation was modeled by Güney [10] using the column-separation modeling assumption proposed by Swaffield [11] and Safwat [12]. Another scientist examining the effect of cavitation occurring during transient flow in PE and PVC pipes was Mitosek [13,14], who showed that an increased pressure reduction is accompanied with gas desorption (reduced pressure oscillations with the increase time period of their existence). Hadj-Taïeb with Taïeb [15] proposed initially a numerical model based on the conservative finite difference method to solve the nonlinear system of hyperbolic partial differential equations describing the transient flow in which the degasification takes place (according to Henry's law). Their study showed that the degasification area is strongly connected with the wall elasticity. The same two authors [16] proposed an alternative modified mathematical model that includes retarded strain and cavitation, which was solved with the second-order finite difference scheme. The mixture density in this model was expressed by means of a nonlinear expression of the liquid volume fraction. Borga et al. [17,18] conducted several transient tests with localized gas cavities in around 200 m long HDPE pipe and concluded that the presence of the leak (or air valves) in cavitating flow induces a greater damping and dispersion of transient pressure waves. Soares et al. [19,20] continued the research of Borga (which was done under the supervision of H. Ramos) and compared the effect of used cavitation models DVCM and DGCM for the prediction of transient flows with cavitation in HDPE pipes. The results indicated that the assumption of the ideal gas law (DGCM) is more appropriate than a simple adoption of vapor pressure when the pressure reaches vapor pressure (DVCM) and induces more attenuation and dispersion of transient pressures. For flows with cavitation, a new set of pipe-wall viscoelastic parameters was determined (calibration technique). The unsteady friction losses, pipe-wall viscoelasticity, and wave speed variation due to the formation of localized gas cavities were described only by the creep function. Such an approach lumped all these important phenomena in the coefficients of the creep function. Keramat et al. [21] utilized DVCM and modeled RS uses a modified Kelvin-Voigt model to study the unsteady flow with cavitation in plastic pipes. His model did not include at that time the unsteady friction effects. The main conclusion from the presented simulations (compared to simulated results with Covas [22] and Soares data [19]) is that viscoelastic pipes strongly diminish the dangers of column separation: "First, cavity opening and collapse occur only one or two times instead of tens of times (as inelastic pipes)". Two years later, Keramat and Tijsseling [23] were first to present a numerical model that included all four important phenomena that take place in transient pipe flows: UF, CAV, RS, and FSI. Unfortunately, to date, there are no experimental results that are conducted to check this interesting model in the full extent. In 2018, Urbanowicz and Firkowski developed the foundations of the model presented in this paper [24]. A year later, Urbanowicz et al. [25] compared the DACM and DBCM models, which had not been used before, for the analysis of cavitation flows in plastic pipes, indicating that they both model the phenomena in a similar way, despite the fact that they are characterized by a significantly different mathematical notation.

Recent Progress in Cavitation Modeling in Metal Pipes
Liu et al. [26] analyzed cavitation that can take place in long-distance transport pipelines. The water hammer due to the collapse of air cavities in the pipeline was discussed when the pump unit is shut down due to an incident. The theoretical and numerical analysis pointed out that it is very important to prevent the occurrence of large water hammer loads due to the collapse of flow interruption in such a system. Santoro et al. [27] tested the DVCM model, writing the continuity equation in terms of mass balance instead of volume balance. Such an assumption allowed calculations with appropriate computational fine grids. Additionally, the flow field was assumed to be two-dimensional (2D axial-symmetric flow), in order to evaluate unsteady friction without the need of parameters calibration. This research pointed out that one-dimensional (1D) models are weakly sensitive to grid size, whereas 2D model results are practically grid-independent, and in the opinion of the authors, the 2D model performs better than the 1D ones. Shankar et al. [28] studied the optimal operation of centrifugal pumps to avoid the major harmful issues as cavitation and water hammering. These authors built a system with a cascade parallel pumping setup. The extensive experimental study reveals that the preferable operating region enhances reliability as well as reduces the occurrence of faults. This paper can serve as a reference to VFD pumping systems and paves the way for sensor-less control. Zhao et al. [29] built an experimental test stand to realize a water hammer event with multipoint collapse. The influencing factors and laws of the cavity length and water hammer pressure have been summarized using the experimental data. They also reveal that the initial flow rate and valve-closing speed greatly affect the water hammer pressure rise and cavity length. In their next work [30], the authors presented a new water hammer velocity formula, a new cavity model, and introduced a floating grid method. An in-house program written in C++ confirmed that the simulation results of the new model matched the measured values.
Sun et al. [31] proposed a quasi-two-dimensional transient model coupled with DVCM which, according to the authors' analysis, can provide a better fit than classic 1D solutions. Warda et al. [32] performed three-dimensional computer fluid dynamics (CFD) simulations based on the finite volume numerical approach. The cavitation was modeled with the use of two models: the Volume of Fluid (VOF) and Schnerr-Sauer. They concluded that the 3D model that was adopted is "deemed physically superior to the existing 1D models as it removes the restriction of the 1D models that vapor cavities, when formed, fill the whole cross-section of the pipe without radial variation". Sanín-Villa et al. [33] considered the influence of the convective terms in the momentum and continuity equations (which standardly are neglected). The cavitation problem has been evaluated by use of the DVCM model. In conclusion, they stated that the influence of the convective term is small compared with a simple model where those terms are neglected. Tang et al. [34] used Fluent software to investigate the cavitation flow in the pipeline. A density-pressure model has been implemented into the continuity equation by using the further development of a user-defined function, which gives the possibility of studying the effects of the variable wave speed on the transient cavitation flow. The weakly compressible fluid RANS model (CFD) results agree well with the measured results. Saidani et al. [35] analyzed the temperature effect (in a range from 4 to 95 °C) on unsteady flow with cavitation. These authors simulated single-phase and two-phase transient flows in a hydraulic copper pipe system. The DVCM and DGCM models were used. From the performed simulations, it was evident that the water hammer is considerably sensitive to the temperature, and its proper value needs to be considered at the design stage of hydraulic systems. Yang et al. [36] used a uniform cavitation distribution model in which the critical flow velocity gradients are calculated both in front and at the back of the section and are the sufficient condition to define water column separation. Dynamic meshes were applied for tracking the change of vaporous cavitation. However, multidimensional models are computationally expensive.

Mathematical Model Derivation
The original bubble cavitation flow model was introduced by Shu more than fifteen years ago [8]. Its applicability was limited to systems made of metal pipes. Today, when plastic pipes are replacing traditional pipes (this trend is especially visible in water supply systems), there is a strong need to modify this interesting numerical model. In the model discussed and presented below, a single-phase flow is treated as its special solution.

Momentum and Continuity Equations
The derivation starts from the equation of momentum of two-phase vapor-liquid flow in a freely oriented conduit: in which being the mixture density is calculated: Please note that this starting momentum equation, as well as the set of the continuity equations (Equation (4)), is identical to the one discussed in the IAHR Synthetic Report [37]. After the differentiation and ordering, one gets the following form: For the sake of simplicity, let us assume that the pipe is horizontal, i.e., θ = 0. Thus, the last term on the left-hand side of the above equation (Equation (3)) is zero. Now, let us write the continuity equations written separately for the gas and liquid phase, respectively: where α-volumetric fraction of liquid.
Adding the continuity equations (Equation (4)) together for the respective phases gives: Next, assuming that a homogeneous bubbly flow takes place, then the dispersed vapor phase does have the same velocity as the surrounding continuous liquid phase = = : By making differentiation and ordering in Equation (6), the following result is obtained: Dividing Equation (7) by , the first useful form of this equation is derived: However, when one multiplies Equation (7) by and then divides by A, we get the second useful form of Equation (7): Using the above Equation (9) and the fact that the analyzed system is horizontal, the momentum Equation (3) can be reduced to a simpler form: Please note that Equation (10) reduces to the equation of a single-phase flow (continuous liquid phase) when no cavitation occurs (mean values of pressure in an analyzed cross-section are larger than the vapor pressure).
The next step is to derive the continuity equation. From the works [22,37,38], it follows that in bubble flow and plastic pipes, the fluid elasticity (separately defined for liquid and vapor) and the pipe deformation can be defined as follows: where Ξ = -enhanced pipe restraint factor. The first two equations represent liquid and vapor elasticity, respectively, whereas the third one defines pipe deformation (the right-hand one). The derivation of the first two ones is straightforward. This is not the case for the third one-its derivation is presented in Appendix A. The total derivative of Equation (2) mixture density is: When Equations (11) and (12) are used in the continuity equation, Equation (8), one gets: In addition, a constant pressure wave speed is assumed. In the proposed model, the value of the speed will be assumed for the steady flow occurring before the water hammer event. Then, there is only pure liquid phase ( = 1). The last term of the square bracket in Equation (13) vanishes, and the formula under square bracket reduces to: which is the pressure wave speed of the pure liquid phase. The above wave speed equation includes elastic effects of the fluid ( ) and of the pipe wall ( ). Enhanced pipe restraint factor Ξ is calculated in a different way in thin (( / ) < 25) and thick-walled pipelines [2]. Equation (14) governs the final form of the continuity equation for unsteady flows in plastic pipes: In non-slip flow conditions, the proportion of the dispersed phase is of a statistical nature; i.e., the volumetric concentration and the mass are equal to the corresponding dynamic shares-the transport concentration and the degree of dryness [39,40]. Then, the following relationship applies in non-slip flows: = / , and the final set of fundamental equations (appropriately momentum and continuity) is as follows: The term / indicates the difference between the velocities of the liquid and vapor phase.

Wall Shear Stress and Retarded Strain
The wall shear stress in transient pipe flow can be calculated with the help of convolutional theory. Zielke [41] for laminar flow and later Vardy-Brown [42] for turbulent flow presented an initial version of this equation. For homogeneous bubble flow, the mixture density should be taken into account: The function ( − ) [-] is a so-called weighting function. In our work, this function is identical to the functions used in other cavitational and single-phase flow models.
For detailed form and more information about weighting functions, please refer to paper [43]. The numerical modified effective solution of the above convolution integral, which is used in this work, is based on the improved solution for single-phase flows [44,45]: The and coefficients are determined from the analytical formulas presented in a recent paper [47], and Δ̂= ∆ is dimensionless time. In the case of turbulent flow, these coefficients must be scaled in accordance with the guidelines presented in paper [48]. The next step is to evaluate the partial derivative of retarded strain using a convolution integral. According to [49,50], the retarded strain can be written in a simpler form than the original one [38]: where (  . The above equation (Equation (21)) may be written in a simpler form: where: The use of convolution integrals (Equations (17) and (20)) in the basic system of Equation (16) results in:

Numerical Solution for Inner Nodes
The convective terms are less important and are omitted from the main analyzed set of equations. This procedure will ensure that the use of interpolation, which greatly affects the numerical solution [2], is excluded. The simplified equations of continuity and momentum are identified as and , respectively: Combining linearly = + , these equations by the unknown multiplier gives: By examination of the above Equation (27) with the definition of total derivatives, it can be noted that with: it becomes the ordinary differential equation: The solution of Equation (28) yields two particular values of , By inserting the above values into Equation (28) the particular relation between and is given as follows: This leads to a set of two equations on positive and negative characteristic lines: From Equation (2), it follows that: Considering the above Equation (34), the third term on the left-hand side in Equations (32) and (33) can take the form: Shu [8] writes the partial derivative of the density of the mixture in a logarithmic form: Note that under the logarithm, we have the division of by , but the could be exchanged to any other constant, and the above equation would be satisfied anyway.
From both characteristics at the inner node, the following explicit system of equations is obtained: The system can be further rewritten after introducing: where: By transforming the system of Equation (39) in a way that the parameters searched for a given inner node D of the characteristics grid ( Figure 1) remain on the left-hand side, one obtains: where as well are time-dependent functions that are iteratively calculated using the values known from the previous numerical time step: From the above equations, the final solutions for the inner node D (Figure 1) of the grid of characteristics are obtained. The solution for the mean velocity at cross-section is: and the solution for the pressure is: The analysis of the above formulas shows that in order for > (note that then = ), the condition that ≥ must be met. When Having the instantaneous value of the mixture density , the vapor density , and the liquid density , the instantaneous value of the liquid phase concentration should be determined from the formula (Equation (2)):

Boundary conditions
The next step is to solve the boundary conditions. According to Figure 1, the instantaneous closing valve of an RVP system is at the left-hand side of the system (x = 0). The valve boundary condition is derived from the negative characteristic: Please note that the value of is based only on known values from the previous time steps. The velocity at the valve section for time > 0 has zero value, i.e., = 0 (closed valve). The above Equation (49) takes the form: which finally reduces to: When the pressure at this boundary is higher than the vapor pressure , then the natural logarithm is equal to 0 as = ; this means that there is no cavitation when < 0, and the pressure at the valve section can be calculated from the following equation: Otherwise, when ≥ 0 and = , then the bubble mixture density and volumetric fraction of the liquid, respectively, should be calculated as follows: Next, the dynamic viscosity of the homogeneous bubble mixture using Dukler's formula [46] should be calculated: At the opposite end (x = L) of the RPV system, Figure 1 is the reservoir. At the crosssection connecting the pipe with the reservoir, the pressure is assumed to be of constant value, i.e., = during the complete transient event associated with the analyzed water hammer phenomenon. As the pressure does not pulsate at this cross-section, the retarded strain is neglected here. The final equation for the velocity pulsation at this section in which the pressure is always higher than the vapor pressure > is:

Experimental Verification of New Model
In order to demonstrate the effectiveness of the newly presented model, in this section, the results of simulation tests will be compared with the experimental results presented by Güney [10]. The Güney experimental test stand located at the INSA research center in Lyon (France) was a simple system consisting of three main components: reservoir-pipe-valve ( Figure 2). In the analyzed RPV system, in steady flow, water flowed directly into the atmosphere. The pipe had a total length of L = 43.1 m and an internal diameter D = 0.0416 m (the wall thickness of the pipe was e = 0.0042 m). The test pipe was made of low-density polyethylene (LDPE). The experimental tests of the water hammer forced by the sudden (momentary) closure of the valve (shutting off the flow) have been carried out for five different temperatures of the flowing liquid (water). In Table 1, the parameters required to simulate the analyzed unsteady flows with cavitation are tabulated. It can be seen that although the initial flow velocity was similar, due to the change in viscosity, the value of the L=43.1 m Reynolds number increased with the temperature increase (almost twice as high for 05: ≈ 82000 than for 01: ≈ 45500). After the temperature change, not only do the parameters related to the flowing liquid change (Table 1) but also the values of the parameters representing the mechanical properties of the pipe; thus, it is necessary to compare their values ( creep compliances and retardation time coefficients values are presented in Table 2). Güney used the time-temperature superposition principle (also known as time-temperature reducibility) to derive his creep compliance functions for different temperatures. During initial simulations, complete Güney creep compliance functions were used (three exponential terms) that can be found in the works [10,52]. As the initially obtained simulation results indicated that this creep function is a source of simulation error, we had a detailed look at the original coefficients. We noticed that the corresponding creep compliance values of small retardation times ( < 1.5 • 10 [ ]-original and coefficients) are out of the frequency range of the used dynamic viscoelastometer RHEOVIBRON. Filtering out this coefficient for small retardation times (rejecting from the analysis original and coefficients without changing all other experimentally defined creep coefficients) helped to receive corrected comparisons results. The creep functions for LDPE have different characteristics (Figure 3) than those for the typical currently used plastic material, namely HDPE. The LDPE material has higher values of creep compliance than the HDPE material. Additionally, we may see (Figure 3) that an increase in temperature increases the creep compliance values. The HDPE traces which are presented for comparison in Figure 3 were obtained experimentally by Covas et al. [22]. The pressure wave speeds were estimated based on the empirically observed duration of the first pressure amplitudes. Their values summarized in Table 1 enabled the determination of (see Table 2) from the transformed formula of the pressure wave speed: where = 0.97; Ξ = = 9.61.
The method of characteristics was used with a constant number of reaches N = 64. The selected number of reaches meets the computational compliance criteria discussed in paper [53], i.e., N > 10. Extra simulation studies performed during the preparation of this paper whose purpose was to investigate the impact of the number of reaches showed that there are no significant differences between the results of N = 16, 32, and the selected 64.   The qualitative analysis of the obtained results ( Figure 4) indicates the following: -In systems based on plastic pipes, as the temperature of the flowing medium increases, the maximum value of the pressure at the first amplitudes decreases (assuming a similar value of the initial velocity in the steady flow just before the quick valve closure). The above is due to the decrease in pressure wave speed with increasing temperature; -The decrease of pressure wave speed caused by the increase in temperature is also responsible for the change in the frequency of the water hammer itself. Based on the research carried out, one may notice that for a higher value of the pressure wave speed, the number of pressure amplitudes appearing in the same time interval increases (Figure 4a-five amplitudes), while for a small value, this amount is smaller (Figure 4d and Figure 4e-four amplitudes); - The omission of unsteady hydraulic resistances negatively affects the modeled waveforms, which are overestimated starting from the second amplitudes (Figure 4a-e).
Large discrepancies are visible in the modeled pressures at the peaks and at the valleys of these pressure amplitudes; - The modified proposed numerical solution modeled the first peak of pressure visible at the beginning of all tops of first amplitudes (Figure 5a) as well the small peak at the beginning of second amplitudes (Figure 5b). This proves the physics of the analyzed phenomena; -In 03, 04, and 05, where only single column separation takes place (after the first amplitude), it can be seen that the phase shift of the simulated pressure increased over time. This behavior can be explained by the fact that in a real situation, the pressure wave speed does not remain constant during the entire transient event but rather slightly changes. The change of the pressure wave speed is not included in the current version of the mathematical model, as it would force the use of interpolation, which would introduce additional numerical damping [2]; -Although qualitative studies indicate the advantage of the model taking into account unsteady friction, it is necessary to carry out quantitative studies to confirm the above hypothesis. Such research will be carried out in the next Chapter 4; - The largest model discrepancies occurred in the runs carried out for 02, 03, 04, and 05 at the top of the first amplitude. In the experimental studies in the final phase of the pressure increase (first amplitudes), no increase in pressure was observed as in 01 (Figure 4a), while from what we can see in Figure 4b to Figure 4e, such an increase was modeled by the numerical model. This increase was not influenced by the way of taking into account skin friction (quasi-steady or unsteady resistances); hence, the applied creep functions were responsible for them.

Quantitative Analysis of Results
In this section, quantitative research is performed whose role is to define and determine important criteria parameters of the analyzed flow. It is difficult to find any favorable quantitative method in the literature on the subject of transient pipe flows. Here, we present a new methodology that results in two criteria parameters. The role of the final

Pressure [Pa]
qualitative parameters is to determine the compliance of the simulated histories with respect to the experimental ones in a simplified mathematical way.
A MATLAB subprogram was written to search for maxima (peaks) and minimum values of pressure histories and their occurrence times (calculated from the beginning of the analyzed transient state). In a demonstrative way, Figure 6 illustrates the working idea of this proposed "collecting" subprogram. As can be seen, the pressure drops to the saturated vapor pressure were not taken into account, as the final results would be false. Additionally, when determining the time compliance, the time of the first pressure peak at the first amplitude was omitted, as it would also cause the final result to be distorted.
The pressure compliance parameter determining the compliance of the maximum and minimum simulated pressures is calculated by the following formula: where , -simulated maximal and minimal pressures and , -experimentally predicted maximal and minimal pressures The time compliance parameter that determines the time fit of subsequent simulated amplitudes was calculated using the following formula: where , -times of occurrence of maximal and minimal simulated pressures and ,experimentally observed times of maximal and minimal pressures (note that times , and , representing the maximum at first amplitude are not taken into account, because in some cases, their fit could distort the whole analysis). The degree of simulation compatibility increases with decreasing values of the above coefficients. In Table 3  The results of quantitative research indicate the following: -Unsteady friction losses contribute to a significant reduction of pressure compliance errors (parameter -values of simulated pressures). When only quasi-steady resistances were taken into account, the average error from all the tests carried out was about 18%, while when the model of UF losses was taken into account, then the average error was about 6.5%; i.e., almost three times smaller; - The unsteady friction also influences the second analyzed parameter, i.e., the phase compatibility . However, in this case, there was a slight increase in the value of the time fit error of the modeled waveforms. When the quasi-steady nature of the resistance was taken into account, the average error was 0.71%, while when the unsteady nature of the resistance was taken into account, the average error was 0.95%. The difference is very small and can be neglected; however, it is recommended to use models of unsteady friction when the experimentally obtained creep functions are used during modeling; - The deterioration of the quantitative parameter , which was noted in the previous paragraph, after taking into account the unsteady hydraulic resistances, prompted us to analyze another quantitative parameter, which is the duration of the cavitation phenomenon at the analyzed cross-section (cross-section at the valve). From the data presented in Table 3 (∆ ), it can be seen that the duration of the cavitation phenomenon at the cross-section at the valve decreased with increasing temperature (decreasing the speed of pressure wave propagation). It can also be seen that the numerical model which takes into account unsteady friction predicts the duration of cavitation areas slightly longer than it was in the experiments. The quasi-steady resistance model overestimated the duration of cavitation quite significantly.

Conclusions
The paper presents a modified unsteady discrete bubble cavity model DBCM. The new model is developed in a very simple form, which makes it easy for implementation in commercial programs for unsteady pipe flow analysis. The new model takes into account three very important phenomena: unsteady wall shear stress, vaporous cavitation, and pipe-wall retarded strain. The latter mentioned phenomenon occurs in plastic pipes.
The conducted comparative studies have shown that with the help of the presented model, it is possible to simulate pressure and velocity waveforms in which vapor areas appear as a result of the cavitation phenomenon in plastic pipes. It was noticed that the influence of unsteady friction was much smaller than the influence of retarded strain. An innovative method of calculating the convolutional integral describing the retarded strain was applied by using the analogy to the convolutional integral defining the wall shear stress (Schohl's method). Taking into account the commonly known boundary conditions related to the method of characteristics enables the use of the novel model in complex networks: water supply, oil hydraulics, heating, etc.
The modified solution is an alternative to the two commonly used transient cavitating pipe flow models, namely the DGCM (discrete gas cavity model) and the DVCM (discrete vapor cavity model). In our future work, we are planning to execute broad comparisons of the presented new model with the existing ones.
The use of the experimentally determined creep functions (obtained by Güney) showed that such functions can be an alternative to the calibration methods commonly developed today. This indicates that the creation of the so-called "maps" of creep function curves for various polymeric materials currently used in the world for pressure pipes will significantly help designers at the design stage and will enable the study of the most dangerous unsteady cases "a priori". Presentation of such "maps" obtained for different temperatures should be a priority of the current scientific research. Data Availability Statement: All codes generated during the study and experimental data are available from the corresponding author by request.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
These changes in plastic pipes differ from those reported in elastic pipes. The relation between the cross-section derivative and strain derivative is: The circumferential strain can be decomposed into an instantaneous elastic strain and a retarded strain : Using the above decomposition (Equation (A3)) in Equation (A2) gives: Typically, the instantaneous strain , which is assumed to be linear-elastic, can be related to the hoop stress as follows: where -elastic component of the hoop stress. The hoop stress is also related to the fluid pressure and the ratio between the pipe inner diameter and wall thickness: Since in most practical applications /2 ≪ /Ξ, then: This equation has been used in the manuscript during the derivation of the continuity equation-see Equation (11).