Application of the Method of Characteristics with Unsteady Friction Model for Numerical Modeling of Transient Flow Induced by Valve Closure

The paper presents a comprehensive methodology for simulating transient flow in pipeline systems induced by valve closure, using the method of characteristics with an unsteady friction model. The research focuses on the instantaneous acceleration-based (IAB) model, a mathematical model employed to describe water hammer behavior in pipeline systems. The methodology includes the development of a mathematical model based on governing fluid dynamics equations, numerical simulation using the proposed model, and validation against experimental data obtained from laboratory-scale pipelines. The study compares the performance of steady and unsteady friction models, revealing the limitations and strengths of each in simulating water hammer events. The paper also discusses the estimation of the damping coefficient (k) using trial-and-error and Reddy's analytical method, and the influence of numerical parameters on the IAB model performance. The numerical results demonstrate good agreement with experimental data, validating the proposed model's accuracy. The methodology presented in this paper can serve as a valuable tool for analyzing and designing pipeline systems subject to water hammer phenomena. It provides insights into transient flow characteristics induced by valve closure and assists in identifying appropriate mitigation measures to prevent damage to the pipeline system.


INTRODUCTION
Water hammer, also known as hydraulic transient, refers to a phenomenon that occurs in closed conduit systems when there is a sudden change in flow rate or pressure.This can lead to significant pressure surges, causing damage to pipes, valves, and other components of the system.It is essential to understand and mitigate the effects of water hammer to ensure the safety and reliability of pipe systems [1], particularly in critical infrastructure such as water distribution networks and power plants.Considering that water hammer is a complex and dangerous phenomenon that requires analysis and damage assessment to find appropriate solutions to minimize it to a certain extent so that pipes can withstand it [2], researchers have studied water hammer for over a century to analyze its behavior and characteristics.Initially, water hammer analysis relied on simple methods such as the use of the wave equation and experimental testing.While these methods have provided great insights into the fundamental physics of water hammer, their capability to predict the behavior of complex systems has been limited [1].Given the limitations of these traditional methods, numerical modeling has emerged as a powerful tool for simulating and analyzing transient events in pipe systems.It provides valuable insights into the underlying physical processes and helps engineers design and optimize pipe systems to reduce the impacts of water hammer.A key aspect of numerical modeling is the selection of an appropriate numerical method and friction models, which are critical for accurately predicting transient behavior in pipe systems [3].One of the main methods used in numerical modeling is the Method of Characteristics (MOC), which is a powerful tool for analyzing water hammer events.This method entails converting the partial differential equations (PDEs) that govern mass and momentum conservation into ordinary differential equations (ODEs) along the characteristic curves, facilitating the analysis of the phenomena [4].Various studies have compared different numerical methods and friction models for water hammer simulations.Kjerrumgaard Jensen et al. [5] investigated water hammer using the MOC and found that unstable friction models, such as Instantaneous Acceleration Based and Convolution Integral models, performed better in capturing the pressure wave compared to the semi-stable model.Similarly, Abdeldayem et al. [3] compared friction models using commercial software and concluded that the unsteady friction model is the most suitable and effective for simulating water hammer, delivering excellent performance in capturing the water hammer wave.Bertaglia et al. [6] conducted a comprehensive comparison of the MOC and the Finite Volume Method using the Quasi-Steady friction model.Their results showed good agreement between all numerical methods and experimental data.On the other hand, Zeidan & Ostfeld [7] compared the MOC with the Wave Characteristics Method (WCM), which is computationally lighter but less accurate.To enhance the accuracy of the WCM, the researchers proposed a new parameter set and calibrated it using available experimental data.However, this method requires initial calibration based on laboratory tests to obtain accurate and suitable results.Duan et al. [8] studied the influence of valve closure time, initial flow rate, and gas content on pressure wave propagation and water hammer effects.They found that increasing valve closure time reduces the maximum pressure wave and delays its occurrence at the pipe ends.Additionally, higher flow rates lead to increased maximum pressure, while increased air content reduces pressure wave oscillation.Kubrak et al. [9] investigated water hammer in a pipeline system comprising steel and high-density polyethylene (HDPE) pipes with different diameters.Their findings revealed that incorporating HDPE pipes into the system can mitigate pressure surges caused by valve closure and reduce the number of wave repetitions during the same time period.These friction models play a crucial role in simulating the pressure wave propagation and energy dissipation within pipe systems.The steady friction model considers only the wall friction losses, while the unsteady friction model takes into account the dynamic changes in fluid flow properties during transient events [10].This research paper concentrates on the development and validation of a numerical model designed to simulate transient flow resulting from valve closure in pipe systems.The performance of steady and unsteady friction models in predicting transient events is investigated, and their accuracy and reliability are compared using experimental data from the literature.The primary goal is to offer a comprehensive understanding of the strengths and limitations of different friction models when simulating water hammer events.Additionally, insights into selecting suitable numerical parameters for precise and efficient numerical simulations are provided.The knowledge gained from this study will contribute to the advancement of more reliable and robust pipe system designs capable of mitigating the adverse effects of water hammer events.

2.
RESEARCH METHODOLOGY This paper starts by providing an overview of the development of the numerical model, including the MOC, friction models, and governing equations.The paper proceeds to discuss the simulation results, conducting a comparative analysis of the performance between steady and unsteady friction models in simulating water hammer waves.Additionally, it explores the estimation of the damping coefficient through trial-and-error and analytical methods.To compare the accuracy of different friction models, the Root Mean Squared Error (RMSE) has been employed as a statistical analysis tool.The calculation of RMSE requires determining the square root of the average of the squared discrepancies is in which   and   denote the predicted and observed values, respectively, and  is the total number of observations.A lower RMSE value signifies a better fit of the model to the experimental data, thus indicating a precise prediction.The implementation of the RMSE in this study assures a reliable and comprehensive metric for assessing the precision of the friction models employed in water hammer analysis.

3.
NUMERICAL MODEL DEVELOPMENT For water hammer modeling, numerical model development refers to the process of creating a computational framework for simulating and analyzing transient events in pipe systems.This involves selecting an appropriate numerical method, incorporating friction models, formulating governing equations, equations interaction in boundary conditions, and acquire experimental data [11].The development of an accurate and efficient numerical model is crucial for understanding and mitigating the impacts of water hammer on pipe systems [12].This research paper aims to investigate the numerical modeling of transient flow induced by valve closure in pipe systems using both steady and unsteady friction models through the application of the MOC.

a. Numerical Method
For this study, we utilized the MOC to solve the governing equations of unsteady pipe flow.The MOC is a popular numerical technique used in water hammer, or hydraulic transient analysis to solve PDEs that arise in the study of unsteady fluid flow in closed conduits by transforming it into a set of ODEs along characteristic lines, or curves, in the space-time domain [13 & 14].These characteristic lines represent the path that information (pressure and flow rate) propagates through the system.By solving the ODEs along these characteristic lines, the method provides an accurate and efficient solution to the unsteady flow problem.A crucial aspect of the MOC is the selection of the Courant number   , which is the ratio of the fluid's wave speed to the numerical speed (∆x/∆t) [4].This dimensionless number affects the stability and accuracy of the numerical solution, the Courant number should generally be less than or equal to 1 [15].In this study, we set the Courant number to 1, which is a common value used in water hammer analysis.This choice strikes a balance between numerical stability, accuracy, and computational efficiency.The appropriate choice of the Courant number is essential for obtaining a reliable and accurate numerical solution, ultimately contributing to the success of the model in representing the physical system. Here Despite these different approaches, determining the appropriate coefficient for specific situations remains a challenge in the application of the IAB model.In this study, the trial-and-error technique will be employed to determine the value of the IAB damping coefficient and compare the results with the analytical method proposed by Reddy et al. [24].

c. Governing Equations
To obtain the governing equation of the MOC with unsteady friction models, we start by examining the following set of one-dimensional PDEs that describe unsteady flow in a pipe.Our objective is to determine a linear combination of these equations: Mass Continuity equation: Momentum equation: Here, H is the piezometric head.
The Linear Combination of Equations: Where,  is an unknown multiplier [4].By using Eq.7 and Eq.8 in Eq.9, will get: To convert the governing PDEs into ODEs we will use the material derivative concept: Here,  is the fluid characteristics (, ) By using material derivative concept in the Eq.10, we will get: Where,  is the fluid flow rate and  is the pipe area.By the following steps, we will find the slope of characteristic lines, which represents the rate at which information propagates through the fluid along the characteristic lines: Depending on what the ∅ value is it, the slope of characteristic lines will be as shown in table below.
Table 1: Possible cases for slopes of the characteristic lines By determining the slopes of the characteristic lines, we can introduce two distinct sets of lines, C+ and C-, along which the ODEs can be solved as shown in Figure 1.By substituting  (+) value in Eq.12, we get: Here we will assume that last term of this equation will not change by time because there is no explicit solution to this term except through approximation.By integrating of this equation, we get: Using this approach, we can derive the equation for the positive characteristic line  + , which will be utilized to determine the discharge at the P node in time  + ∆ as shown below: Since the fluid properties are known only at the nodes within the characteristic grid in the spacetime domain, we must employ interpolation techniques to determine the characteristics at points where the characteristic lines intersect the grid as follow: Here,  is a fluid properties (, ).
• When ∅ < 0,  (+) = For this case, the positive characteristic line equation  + that will be used to determine the discharge at the P node in time  + ∆ is shown below: Similarly, the negative characteristic line equation  − , which we will be used to determine the discharge at the P node in time  + ∆ is:

Boundary Conditions
To simulate transient events in pipe systems, this requires the interaction of the governing equations with boundary conditions in the appropriate manner.A detailed explanation of the interaction of equations with these boundary conditions for case of water hammer induced by downstream valve closure is as follow: i.Upstream Boundary Condition (Reservoir): At the upstream boundary, the reservoir's piezometric head (H) is assumed to be constant as it typically has a large volume compared to the flow in the pipe system (Figure 6).At the downstream boundary, the valve controls the flow, influencing the discharge (Q) and piezometric head (H) in the pipe system (Figure 7).Then Eq.24 will be: • When ∅ ≥ 0 Form Eq.14 head will be: By substitute   value in Eq.25, we get: Neglecting the negative root of   , we get: To determine the fluid properties at the downstream boundary nodes, we will use the following equations Eq.13, Eq.26, and Eq.27, respectively.
• When ∅ < 0 Form Eq.19, head will be: By substitute   value in Eq.25, we get: It is important to mention here that there are no negative characteristic lines present in the downstream boundary nodes, and ∅ value at the downstream boundary nodes will be found by the following approach: iii.Interior Nodes: Fig. 8 Interior nodes boundary conditions.
At the interior nodes (Figure 8), both the C+ and C-characteristic lines intersect.For this, the governing equations must be simultaneously solved for these lines to determine the fluid properties at these nodes, accounting for any changes in pipe geometry or fluid properties.The resulting solution provides the updated discharge (Q) and piezometric head (H) values at the interior nodes for the next time step.There are four cases to determine fluid properties at the interior boundary nodes, the cases are: • Case 1: When ∅ ≥ 0 at C+ & C-lines (Figure 9).By adding Eq.14 to Eq.16, we obtain the following equation.
In addition, by using Eq.16, Eq.13 and Eq.15, we will find the value of   ,  + and  − , respectively:  By adding Eq.14 to Eq.22, we obtain the following equation.
Eq.19, Eq.18 and Eq.15 have been used to find the value of   ,  + and  − , respectively: These experiments offer a diverse range of scenarios that allow us to assess the performance of the numerical models for transient flow using Steady and Unsteady Friction Models by the Characteristics Method.In the following sections, we will present the results of our numerical simulations and compare them with the experimental data to evaluate the accuracy and reliability of our models under different operating conditions.

f. Models Programming
In contemporary hydraulic engineering, programming has taken a pivotal role in modeling complex physical processes, including phenomena like the water hammer.With the integration of numerical methods and advancements in computing infrastructure and programming languages, intricate and highly accurate simulations have become more streamlined and readily accessible.In the past, the colossal complexities and computational demands associated with numerical methods rendered them beyond the realm of practical use.Researchers primarily depended on simpler analytical techniques or empirical formulas, which, while easier to handle, often lacked precision.However, the advent of modern computational technologies has placed numerical methods within reach, paving the way for detailed and complex simulations of phenomena such as the water hammer with heightened accuracy.The algorithm employed to develop the Python code for simulating the water hammer in this research, is clearly detailed in the figures below (Figures 13 and 14).The procedure initiates with the designation of inputs, encompassing elements like discharge, pressure head, wave speed, and the like.Subsequently, an array is established for fluid properties variables, providing a storage mechanism for the characteristics at each individual node.The initial step involves identifying the fluid properties at the nodes in the grid at the initial time.This is followed by deducing the fluid properties for the ensuing temporal step for all nodes, depending on the boundary conditions and the ∅ signal.After the calculation of fluid properties, the data are preserved in an Excel file.The process then transitions to the subsequent temporal step, instigating a repetition of the aforementioned procedure.

RESULTS AND DISCUSSIONS a. Comparison of Steady and Unsteady Friction Models
In this section, we validate the numerical models by comparing and discussing the results of the steady and unsteady friction models with the available experimental data in the literature to highlight the strengths and limitations of each in simulating water hammer events.For the steady friction model, the results, as shown in Figure 15, indicate that the model effectively captured the first two peaks of the water hammer event when compared to all experimental data.However, the model showed less accuracy in describing the actual amount of damping that occurs during the water hammer propagation, which can be attributed to the inclusion of only wall friction losses in the energy loss calculation [  The steady friction model successfully detects the first pressure peak in the transient event in Figure 15 can be attributed to the transient event's insensitivity to dissipation at its onset [19], while the detection of the second peak as in Figure 15 (a to d) can be associated with a complex processes occurring in the transient event, such as cavitation, which involves the boiling and formation of vapor bubbles within the liquid due to pressure reduction to the vapor pressure level or below.The explosion of these bubbles could produce pressure wave propagation similar to the first pressure wave, as is the case in Figure 15, or even higher, as in Figure 16, generating what is known as a two-phase water hammer.Therefore, we observe the second pressure peak to be close to the first pressure peak and then stabilize in the remaining pressure peaks of the transient event due to the pressure not dropping to the required level [21].
Fig. 16: The effect of cavitation on pressure waves heads.
In 15 (e -f), the presence of a second peak can potentially be attributed to the specific pipe length utilized in the experimental setup.This observation is supported by the understanding that shorter pipes tend to result in shorter travel distances for pressure waves, resulting in reduced frictional losses along the pipeline.
In contrast, the results of the unsteady friction model showed a more accurate representation of damping and a more comprehensive and precise depiction of the water hammer phenomenon throughout the duration of the event.The model succeeded in capturing the actual damping behavior of pressure waves accurately, which can be attributed to the inclusion of unsteady friction terms in the governing equations to explain and describe the losses resulting from dynamic changes in fluid flow characteristics during the water hammer.This occurred due to the fact that rapid changes in transient pressures and flow rates led to increased turbulent shear stresses between the flowing liquid layers [1] and the resistance of the liquid mass to these changes [19].
The inability of the unsteady friction model to detect the second or third pressure peak in some cases is attributed to the occurrence of cavitation phenomenon after the first wave peak, which led to the generation of a wave higher than or equal to the first wave peak itself, as previously discussed.This is due to the fact that the term cavitation phenomenon was not included in the governing equations, resulting in neglecting its calculations during the simulation.
During the propagation of pressure waves in the unsteady friction IAB model, certain instances reveal a temporal discrepancy in the timing of pressure waves compared to experimental data.This discrepancy tends to magnify over time, leading to decreased accuracy in simulating the frequency and speed of water hammer waves.Consequently, this temporal discrepancy stands as a significant limitation of the 1-k IAB unsteady friction model [27].
Upon comparing the model results with the experimental data at the middle section of the pipe, a consistent pattern was observed, aligning with the previous comparison at the valve, as illustrated 17.The successfully the initial pressure but in accurately the and of subsequent pressure peaks.In contrast, the unsteady model demonstrated the ability to effectively simulate the authentic damping behavior of the pressure wave peaks, providing a more reliable representation of the transient flow conditions.17: Comparison between the experimental data and numerical results at the middle section of the pipe.Although the primary focus in the engineering design of piping systems is generally to prioritize the endurance of the first pressure peak of the water hammer, which is provided by the steady friction model, there are many cases that require adequate attention to secondary pressure peaks as well [28].These include multiple operations such as opening a valve after closing it or starting a pump after turning it off, etc. [4].Furthermore, it is important to gain a better understanding of the mechanisms of dynamic changes that occur during transient events to understand and avoid more complex phenomena, such as the two-phase water hammer, which may generate higher pressure waves than normal conditions [29].After providing a visual comparison of the friction models through figures, the subsequent table offers a detailed statistical assessment of these models.Table 3 provides a comparative analysis between the friction models, utilizing the RMSE as a robust statistical measure.This method allows for an accurate quantitative assessment of the models, bringing attention to the nuanced differences in their performance that might not be immediately apparent in the graphical representations.Therefore, it offers a deeper understanding of the precision of the two friction models in simulating the water hammer phenomenon.Choosing the damping coefficient value has a significant impact on the accuracy of simulating pressure wave damping and propagation, which is one of the main weaknesses of the IAB friction model.In this study, the trial-and-error method is used to estimate the damping coefficient for a specific pipeline system.The process begins with an initial value for k and adjusted it repeatedly until the simulated pressure wave closely matched the experimental data.Although this method is simple and easy to implement, it is considered one of the most important techniques to obtain the optimal k value.To overcome the limitations of the trial-and-error technique, which requires prior experimental data, the analytical method proposed by Reddy et al [24] is applied to estimate the damping coefficient.This method is an improved version of Vardy & Brown's method [23], which improved specifically for use in smooth pipes, and when the valve closed instantaneously at downstream of pipe.This method involves calculating the k coefficient based on Vardy's shear decay coefficient (C * ), which in turn depends on the Reynolds number.The experimental data from (Bergant et al. 2001) research paper will not be considered in the context of the Reddy method, as this technique is proposed for analyzing turbulent flows with a Reynolds number exceeding 4,000.The given formula was used to calculate k coefficient for pipeline system of the USC experimental data in the literature and compared the results with those obtained using the trial-and-error technique.Upon comparison, we found that Reddy et al. [24] method effectively described the pressure wave damping, closely resembling the trial-and-error technique in three of the four cases, as depicted in Figures 18 (a, c and d).However, this approach yielded less accurate results in the fourth case, shown in Figure 18 (b).This discrepancy can be attributed to the more significant impact of the cavitation phenomenon in Experiment 2, which may cause an increase in the generated pressure wave magnitude compared to when it is absent.In this case, the second pressure peak surpasses the first peak, which is a dominant outcome of this phenomenon.Consequently, a larger damping coefficient was employed using the trial-and-error technique when attempting to simulate this case compared a scenario without the phenomenon, in order to compensate for the lack of inclusion of the cavitation term in the governing equations.

CONCLUSION
This research paper aimed to develop and validate numerical models for simulating transient flow in pipe systems using steady and unsteady friction models based on the MOC.The study focused on analyzing and comparing the performance of different friction models in predicting transient events caused by valve closure.The results demonstrated that the steady friction model effectively captured the first pressure peaks of the water hammer event in most cases.However, it showed limitations in accurately describing the damping and propagation of subsequent pressure peaks.On the other hand, the unsteady friction model, specifically the 1-k IAB model, provided a more accurate representation of damping and a comprehensive depiction of water hammer events throughout the duration of the transient flow.Additionally, the unsteady friction model was found to be capable of simulating water hammer events at different locations in the pipe, including the middle of the pipe.This further demonstrates the versatility and applicability of this model in capturing transient flow events across various pipe systems.A significant challenge in the application of the IAB model was the determination of the appropriate damping coefficient value.In this study, both the trial-and-error technique and Reddy et al. [24] analytical method were employed to estimate the damping coefficient.The comparison of both methods showed that Reddy et method described the pressure damping in most cases, with some discrepancies arising due to the impact of cavitation phenomenon in specific scenarios.In cases where experimental data is not available, the Reddy et al.'s method can be relied upon to find the damping coefficient in smooth pipes with a sudden valve closure at the downstream end in the case of turbulent flow.This offers a practical and efficient solution for estimating the damping coefficient in the absence of experimental data.In conclusion, this research paper highlights the importance of selecting appropriate friction models and numerical parameters for accurately simulating transient events in pipe systems.The unsteady friction model offers a more reliable and accurate representation of water hammer events compared to the steady friction model, particularly in scenarios where secondary pressure peaks need to be considered.Further research can focus on developing the IAB model by employing a two-coefficient IAB model to eliminate the variation in the propagation speed of pressure waves compared to experimental data.By addressing these limitations and enhancing the accuracy of the numerical models, future studies can contribute to a more comprehensive understanding of transient flow events and their impact on pipe systems.

Fig. 1 Fig. 2
Fig. 1 Characteristic lines in the x-t plan.These characteristic lines represent the paths along which fluid properties change consistently over time, enabling efficient and accurate analysis of the transient behavior in pipe systems.The characteristic lines, C+ and C-, are derived by

Fig. 3
Fig. 3 Negative characteristic line when ∅ ≥ 0 in the x-t plan.

Fig. 4
Fig. 4 Positive characteristic line when ∅ < 0 in the x-t plan.

Fig. 5
Fig. 5 Negative characteristic line when ∅ < 0 in the x-t plan.

Fig. 6
Fig. 6 Upstream boundary conditions.The governing equations for the C-lines are then used to calculate the discharge (Q) at the upstream boundary nodes by depending on the following states:

Fig. 7
Fig. 7 Downstream boundary conditions.The valve's behavior can be represented by a state of orifice states (Small Orifice), which relates the discharge (Q) to the valve opening and other fluid properties.The governing equations for the C+ lines are used in conjunction with the orifice equation to determine the fluid properties at the downstream boundary nodes by the following steps: Dividing the semi-open orifice discharge equation by the fully open orifice discharge equation.    = (  )  √2  (  )  √2  … … … 24

Fig. 13
Fig. 13 Algorithm of steady friction model

Fig. 15 :
Fig. 15: Comparison between the experimental data and numerical results at valve location downstream of the pipeline.
Fig.17: Comparison between the experimental data and numerical results at the middle section of the pipe.Although the primary focus in the engineering design of piping systems is generally to prioritize the endurance of the first pressure peak of the water hammer, which is provided by the steady friction model, there are many cases that require adequate attention to secondary pressure peaks as well[28].These include multiple operations such as opening a valve after closing it or starting a pump after turning it off, etc. [4].Furthermore, it is important to gain a better understanding of the mechanisms of dynamic changes that occur during transient events to understand and avoid more complex phenomena, such as the two-phase water hammer, which may generate higher pressure waves than normal conditions [29].After providing a visual comparison of the friction models through figures, the subsequent table offers a detailed statistical assessment of these models.Table3provides a comparative analysis between the friction models, utilizing the RMSE as a robust statistical measure.This method allows for an accurate quantitative assessment of the models, bringing attention to the nuanced differences in Fig. 18 Comparison of Damping Coefficient (k) Estimation Methods.
,   is the Courant number,  is the wave speed in the pipe, ∆ is a spatial step size, ∆ is a time step size.Brown [23], and subsequently refined by Reddy et al. [24] specifically for downstream valve instantaneous closures in smooth pipes under turbulent transient flow, as shown in Equations 5 and 6.
* ), initially introduced by Vardy & Vo= Initial velocity in the pipe; Hs = Static head at the tank; f values calculated using Colebrook-White Equation in all test except test no.5, f=64/Re because Re<2000.
19], which may not account for the dynamic changes in fluid flow properties during transient events.

Table 3 :
Statistical comparison of friction models . b.