Numerical modeling and simulation of water hammer phenomena using the MacCormack method

Water hammer refers to pressure ﬂ uctuations caused by changes in ﬂ uid speed in hydraulic systems. This research proposes a model using the MacCormack method to study and control water hammer, speci ﬁ cally focusing on managing valve closure to reduce shock wave pressure during transient ﬂ ow. Numerical simulations compare linear and quadratic closure laws for both rapid and gradual valve closure, with the goal of identifying the optimal laws. The ﬁ ndings show that with faster closure times, the overpressure at the valve section decreases linearly in the second quarter of the return period (t4). The application of the quadratic law results in reduced pressures at the valve compared to the linear law, with a maximum pressure difference of 0.05 bar between the highest and lowest values. When searching for the optimal valve closure law to mitigate overpressure, it is found that the exponent ‘ m ’ falls within the range of 1.117 – 1.599. As the slow closing time increases gradually, the range of variation for ‘ m ’ decreases. Furthermore, in the case of underpressure prevention, the exponent ‘ m ’ ranges from 0.95 to 1.419, and the range of ‘ m ’ remains relatively constant with the increase in closing time.


INTRODUCTION
Water hammer refers to the occurrence of a pressure wave traveling at the speed of sound in water supply systems.It arises from abrupt changes in water flow, which can be caused by various factors such as the starting or stopping of a pump, pipe ruptures, opening or closing of valves, and so on (Khalighi et al. 2015;Mungkasi et al. 2015;Jensen et al. 2018;Firkowski et al. 2019).In addition, water hammer can also result from the failure of hydraulic devices like valves and pumps, variations in water demand, human errors, and similar factors (Bergant et al. 2010;Malekpour et al. 2015).According to Twyman (2018aTwyman ( , 2018b)), the recommended approach for protection involves gradually closing valves.However, in the case of automatic valves, a malfunction can result in rapid closure.Han et al. (2022) demonstrated that the closing law employed significantly impacts the pressure value, with a higher water hammer pressure observed when the valve closes rapidly at the beginning of the closing process.Various numerical methods have been proposed to solve the transient flow equation.These include the finite volume method (Christen et al. 2010;Berger et al. 2011), the finite element method (Chen & Lee 2003), and the finite difference method (Liang et al. 2006;Beguería et al. 2009).
Several studies have investigated the influence of different types of valves and network configurations on transient pressures.It was found that extreme pressure values were primarily related to the type of valve, while the network configuration had minimal impact on attenuating transient pressure (Herasymov et al. 2019).Various closure techniques were examined by different researchers.Herasymov et al. (2019) recommended linear closure with a breaking point to reduce surges, while Wan et al. (2019aWan et al. ( , 2019b) ) emphasized the importance of staggered valve closing times in mitigating water hammer.Zhao et al. (2020) demonstrated that flow rate and valve closing speed significantly influenced water hammer pressure, with higher flow rates and faster closing speeds leading to greater surge pressure.Arefi et al. (2021) focused on water transport pipelines characterized by large diameter and high flow rates.They analyzed the effects of factors such as valve type, abrupt valve closing time, friction coefficient, wave velocity, and viscoelasticity.The study found that exclusive use of air valves was not effective, but proper installation of water flow control devices and hydropneumatic tanks along the pipeline helped prevent negative pressure.Jinhao et al. (2021) investigated six key parameters influencing transient flow and found that flow rate and Young's modulus had a positive correlation with maximum pressure, while pipe diameter, valve closing time, and wall thickness showed a negative correlation.Flow rate, diameter, and valve closing time were identified as key factors impacting water hammer simulation.Yong et al. (2022) discovered that extending the valve closing time effectively reduced maximum water hammer pressure.They also observed that pressure fluctuations were influenced by the valve closing speed, with faster closing speeds initially resulting in higher water hammer pressure.Zheng et al. (2022) identified the distinct characteristics of the optimal nonlinear closure law of the valve during fast and slow closure, providing guidance for real-time valve control and design aimed at attenuating surge waves.
Optimizing valve closure strategies and coordinating hydraulic system performance have been the commonly used methods to control the water hammer phenomenon (Bazargan-Lari et al. 2013;Wan & Li 2016;Zhou et al. 2017;Wan & Zhang 2018;Wan et al. 2019aWan et al. , 2019b)).Researchers have focused on developing optimal valve closure laws (OVCLs) to mitigate maximum water hammer pressure, as this pressure depends primarily on the valve closure characteristics.Implementing an OVCL is a useful approach to control extreme pressures, improve operational reliability, and extend the lifespan of water distribution systems.
In their study, Toumi & Sekiou (2024) employed the characteristic method of characteristics (MOC) and a mixed scheme for simulating transient flow.Their objective was to identify the OVCL for mitigating overpressure.The results revealed that for each slow closing time, there exists a unique convex law that effectively reduces maximum overpressure.Furthermore, the evolution of the optimal pressure at the valve is governed by two models: exponential and linear.
The MacCormack method (MacCormack 1982) is an implicit approach proposed within the finite difference framework to solve the complete unsteady Navier-Stokes equations.This method offers several notable advantages, such as unconditional stability and the elimination of the need to solve equations involving block tridiagonal matrices.Several researchers have investigated the application of the MacCormack finite difference scheme in one dimension and have demonstrated its feasibility and robustness as a computational solution for solving the water hammer equation (Tseng & Chu 2000;Wang et al. 2000;Vincent et al. 2001).It is possible to explicitly calculate numerical solutions for the discretized governing equations using these methods, and thus they take up less computational resources than methods involving implicit solutions, which has not yet been highlighted by the current literature (Magdalena & Pebriansyah 2022).Furthermore, the choice of the number of grid points may be important, as well as the specification of the boundary conditions Breuss (2004).
The main objective of this article is to explore the application of the MacCormack method for the numerical simulation of one-dimensional transient flow in an inclined straight pipeline.Specifically, the study focuses on simulating the flow resulting from the sudden closure of a valve located at the end of the pipeline.The research aims to investigate the impact of both rapid and slow valve closures with the goal of minimizing the occurrence of overpressure and underpressure peaks.In addition, the study seeks to identify the OVCLs that effectively protect against both overpressure and underpressure conditions.
The structure of this paper is as follows: The initial section provides a fundamental description of the water hammer phenomenon and presents its governing equation.Subsequently, the 1D MacCormack method is employed to discretize the finite computational domain along with boundary conditions.The second section focuses on exploring the optimal closing law among various commonly used closing laws, quadratic and linear.The third section centers on the search for the ideal closing law and the depiction of optimal flow variations.Finally, the results obtained for both fast and slow closure of the valve are compared and discussed.

Mathematical model
The fundamental hyperbolic equations that accurately describe transient flow phenomena in a hydraulic pipeline are derived consecutively from the continuity and dynamic equations (Streeter & Wylie 1983;Toumi & Remini 2022).
where P is the pressure (Pa); ρ is the density of the liquid, for water ρ ¼ 10 3 kg/m 3 ; u is the velocity in the x direction (m/s); g is the acceleration due to gravity (m/s 2 ); α is the pipe inclination angle; j is the unit pressure drop ; and a is the wave velocity (m/s).
For the case of a buried pipe, wave celerity is done by the following expression (Ait Sliman 2018): where with ρ being the density of the liquid, for water ρ ¼ 10 3 kg/m 3 ; E m is the Young's modulus of metal (Pa); v m is Poisson's coefficient of the metal; e m is the armor thickness (m); K is the modulus of elasticity of water (Pa); r is the pipe radius (m); P is the hydrostatic pressure (Pa), and P a is the backfill soil load (Pa).
The backfill soil load is calculated by the formula given by Marston written as follows: where γ is the soil density (KN/m 3 ), B d is the trench width (m), H is the height of backfill soil above pipe (m), and C d is the load coefficient.
μ 0 is the vertical coefficient friction between the backfill and the trench walls and K is Rankine's coefficient (Moser 2001).
j is expressed by the following relation: where d is the internal diameter of the pipe (m), λ is the coefficient of friction, S is the pipe cross section, and Q is the flow rate.

Transformation of the transient flow equations
The transformation of equations for transient flow in pressure pipes, written as a function of velocity and pressure (see Equation ( 1)), to equations written as a function of piezometric head and flow rate without neglecting the terms u(∂u/∂x) and u(∂P/∂x) in front of the terms ∂u/∂t and ∂P/∂t, and the pipe may have an inclination with respect to the horizontal plane, i.e., sin(α) ≠ 0, must use the following mathematical expressions and introduce them into Equation (1).
where H and Z are the manometric head and the geodesic level, respectively.
Multiplying the first equation by a/g and the second by S, we get These two equations represent the variations in head and flow rate with respect to time in the direction of flow.This transformation was performed because, in practical scenarios, the available data are given in terms of heights and flows, and all intermediate parameters are calculated using consistent units.

Resolution method
The MacCormack method is a widely used discretization scheme for the numerical solution of hyperbolic partial differential equations.This second-order finite difference method was introduced by Robert W. MacCormack in 1969. 1 The MacCormack method is elegant and easy to understand and program.
The MacCormack method consists of using two steps: in the prediction step it uses a finite difference scheme decentered upstream, and in the corrector step it uses a scheme in finite differences decentered downstream.These two steps can be summarized mathematically as follows: • Predictor step: • Corrector step: Applying these two steps to Equation (8) leads to the following expressions: Predictor step: This step allows the prediction of the values of H and Q at time (t þ Δt).
Let us multiply both sides by Δt and rearrange: Corrector step: This step allows the correction of the values of H and Q at time (t þ Δt) obtained in the prediction step.
Let us put where Here, Δt is the time step, Δx is the space step, S is the cross section of the pipe, d is the pipe diameter, λ is the coefficient of friction, i is the space index, and k is the time index.
The application of the MacCormack method diagram to the interior nodes will not cause any problem, but at the boundary limits (on the left and on the right), it results in the following expressions: At node i ¼ 1, the problem only arises in the correction step and the expression giving the flow rate is written as follows: Corrector step: At node i ¼ n, the problem arises only in the prediction step and the expression giving the piezometric head is written as follows: Predictor step: These two expressions clearly highlight the appearance of two fictitious nodes (i ¼ 0) and (i ¼ n þ 1) with coordinates (Q 0 , H 0 ) and (Q nþ1 , H nþ1 ), respectively.For this and to calculate the values of the flow rate and the head at the nodes (i ¼ 1) and (i ¼ n), the MOC is used to overcome the fictitious nodes.Figure 1 elucidates the combination of the application of the MacCormack method and the MOC.

Problem statement
To illustrate the mathematical model derived, the RPV (Reservoir-Pipeline-Valve) system is adopted (Figure 2).This system comprises a large reservoir positioned at the upstream section of the inclined pipeline, and a valve situated at the downstream end that discharges into the atmosphere.The reservoir (R), equipped with an overflow at an elevation of 200 m, supplies a flow rate of 350 l/s through a gravity pipe (P) spanning a length of 1,000 m and inclined at an angle of 10°.At the opposite end of the pipe, at an elevation of 50 m, there is a motorized valve with controlled closure.The transient flow is generated by the fast and slow closing of the valve.The schematic representation of the problem is depicted in Figure 3.

Boundary conditions and stability criterion
When applying the MacCormack method in finite differences, they appear, at the limits (on the left and on the right), indefinite variables in the domain of interest, this pane highlighted how to get rid of this in definite.
H(2) and Q(2) are quantities calculated at time (t), so where H r is the piezometric head at the reservoir and CN is the negative characteristic of the MOC.
Here,  where the parameter R represents the flow resistance effect, and the elements of this parameter are as follows: a is the wave velocity (m/s) and S is the pipe cross section (m 2 ).
2.4.2.2.At the valve.In the case of instantaneous closing (valve closing time equals zero) Q(n) ¼ 0, the piezometric head at node number N is written as follows: where the positive characteristic, CP, is written as follows: In the case where closing is controlled, the flow at the valve follows, over time, a function depending on the chosen closing law, and the flow, in this case, is written as a function of Q 0 , Δt, Ct, and m by the following expression (Provenzano et al. 2011;Subani & Norsarahaida 2015): where m is the exponent determining the valve closing law, Δt is the time step, k is the time index, Ct is the total valve closing time, and Q k and Q kþ1 are the flow rates at times t and t þ Δt, respectively.
Figure 4 shows the variation of flow over the time for different curves of the possible closing valve laws.
The exponent m, of the expression of Q at node n, generates the infinity of functions of Q, and consequently the search for the optimal value of this exponent turns out to be a major priority for an optimal management.In this case, the piezometric head at the valve is calculated by the following relationship: with The stability criterion of the MacCormack method must ensure the stability of the finite difference method and the MOC, which allows us to write the following expression: where a is the wave velocity (m/s), Δt is the time step, Δx is the space step, N is the number of discretization nodes, and L is the length of the pipe (Goldberg & Benjamin 1983;Shampine 2005;John 2018).The flowchart in Figure 5 gives the schematic representation of the solution to the transient flow problem in pressurized pipelines.
The execution of this algorithm requires knowledge of the characteristics of the materials used, the reservoir, and the volume flow transported by the pipe and the liquid.All the latter characteristics are represented in Table 1.

RESULTS AND DISCUSSION
To begin transient flow calculations, the resolution method has to use the results obtained under steady-state conditions (Figure 6).The results are shown in Table 2.

Case of fast closure
Figures 7 and 8 present the evolution of pressure for rapid and slow closure, respectively.It is evident that during rapid closure, the pressure increases linearly in the first quarter of the shock wave return period, followed by a similarly linear descent in the second quarter of the return period.The quadratic law results in lower pressures at the valve compared to the linear law.The linear law exhibits a high pressure gradient in the first quarter of the return period and a lower gradient in the second quarter.Both valve closure laws show a difference of no more than 0.05 bars between the maximum and minimum overpressure values recorded.

Case of slow closure
This section is focused on exploring different valve closing times (t) ranging from 0.51 t 4 to 10 t 4 s, with a step size of Δt ¼ 0.001 s, to determine the optimal valve closing law that minimizes overpressure and underpressure.To illustrate this process, calculations are performed using two commonly used closure laws: the linear law and the quadratic law.Figures 7(a) and 7((b)) show the time evolution of valve pressure for a closing time equal to three times the return period (3 t 4 ).From these two graphs we can deduce the maximum overpressure and the maximum underpressure.It is clear that for the slow closing time tested (3 t 4 ), the pressure generated under the quadratic law is higher than  that generated under the linear law.Moreover, the minimum pressure is observed in the case of quadratic valve closure.
Extending the previous result, Figures 8 and 9 show the evolution of maximums overpressure and depression, respectively, obtained for different slow closing times ranging from 0.51 t 4 to 10 t 4 with a time step equal to 0.01 t 4 , (i.e., 950 valve closing times).
It appears that as the slow closing time increases, overpressure at the valve decreases and tends asymptotically along the x-axis toward equilibrium pressure.The overpressures generated by the quadratic law are significantly    higher than those of the linear law. Figure 9, which shows the variation in maximum depression, clearly demonstrates the decrease in maximum depression as the valve's slow closing time increases.The depression curve of the linear law appears sinusoidal, unlike the quadratic law, which has a smoother appearance.
3.2.1.1.Search for optimal laws to prevent overpressure.The search for the optimal law that generates the minimum overpressure, for any slow valve closing time, involves testing all possible closing laws by varying the exponent m of Equation ( 20), and recording the value of m for the law that minimizes the maximum overpressure.
To do this, we examined slow valve closing times ranging from 0.51 t 4 to 10 t 4 , using a time step of 0.05 t 4 (191 closing times were tested).Figure 10 illustrates the variation of the exponent m of the optimal law as a function of slow closing times.The curve describing the evolution of m shows a periodic pattern, with each slow closing time exhibiting a unique optimal law that produces the minimum overpressure.For low closing times, the exponent m is between 1.117 and 1.599.As the closing time increases, the range of variation of m decreases.
This curve of values of m of the optimal laws is put under the envelope of two curves, one upper with equation m sup ¼ 1.8055Ct À0,101 (R 2 ¼ 0.9913) and the other lower with equation m inf ¼ À9E À 05Ct 2 þ 0:0017Ct þ 1:2699 (R 2 ¼ 0.9075).The gap between these two curves begins large and narrows over time.3.2.1.2.Search for optimal laws to prevent underpressure.The search for the optimal law to generate the minimum underpressure, for any slow valve closing time, consists also in testing all possible closing laws by varying the exponent m and recording the m that minimizes the highest underpressure peaks.To do this, we examine the slow closing times ranging from 0.51t 4 to 10t 4 , taking a time step of 0.05t 4 (191 closing times tested).
As shown in Figure 12, the optimal values of the exponent m in the case of prevention against depression range from 0.913 to 1.405, and they are located between two envelope curves.The equations of these envelope curves are successively m sup ¼ 1.5288Ct À0.051 and m inf ¼ À3E-05Ct 2 À 0.0004Ct þ 0.9557.In contrast to the results obtained for overpressure protection, the range of variation of the exponent m appears to be only slightly influenced by the increase in the valve closing time and the most of optimal values are between 0.91 and 1.0.For each slow closing time (t .0.5t 4 ), two optimal laws are derived.The first law generates the minimum overpressure and controls the flow evolution during closure, while the second law generates the minimum depression and also governs the flow evolution during closure.
Figure 14(a) displays the variation of flow during valve closure for closing times ranging from 1t 4 to 9t 4 , based on the optimal laws that generate the minimum overpressure.Figure 14(b), however, shows the flow evolution during closure using the optimal laws that generate the minimum depression.

CONCLUSION
The results obtained from studying the linear and quadratic laws demonstrate that as the rapid closure valve increases, the overpressure at the valve decreases linearly during the second quarter of t4, following an initial increase in the first period.In the case of slow closure, it is observed that as the closing time increases, both overpressure and depression decrease, approaching a steady-state pressure.The study also reveals that different laws are optimal for mitigating overpressure and depression during slow closure, and as the duration of slow closure increases, both overpressure and depression decrease.The exponent 'm' of the optimal law, which is crucial for effective protection against overpressure and depression, falls within a defined range determined by two boundary curves.This parameter is dependent on the duration of slow closure, which governs the optimal flow variation.When searching for the OVCL to mitigate overpressure, it is found that the exponent 'm' falls within the range of 1.117-1.599.As the slow closing time increases gradually, the range of variation for 'm' decreases.Furthermore, in the case of underpressure prevention, the exponent 'm' ranges from 0.95 to 1.419, and the range of 'm' remains relatively constant with the increase in closing time.Finally, the findings confirm the accuracy of the MacCormack method in numerically simulating water hammer phenomena.

Figure 1 |
Figure 1 | Cross section of the gravity pipe taken as an example.

Figure 2 |
Figure 2 | Meshing in time and space of MacCormack and MOC methods.

Figure 3 |
Figure 3 | Schematic of the RPV system.

Figure 4 |
Figure 4 | Variation of flow versus time for different closing valve laws.

Figure 5 |
Figure 5 | Flowchart for calculating hydraulic parameters in transient flow.

Figure 6 |
Figure 6 | Simulation result of pressure versus time for linear and quadratic closing laws.

Figure 7 |Figure 8 |
Figure 7 | (a) Simulation result of pressure at the valve versus time for quadratic linear law.(b) Simulation result of pressure at the valve versus time for quadratic closing law.

Figure 9 |
Figure 9 | Simulation result of depression at the valve versus time for linear and quadratic closing laws.

Figure 10 |
Figure 10 | Simulation result of optimal value of the exponent m versus valve closing time in the case of overpressure protection.

Figure 11 |
Figure 11 | (a) Variation of overpressure versus closing time for optimal m.(b) Variation of corresponding depression.

Figure 13 |
Figure 13 | (a) Variation of depression versus valve closing time for optimal m.(b) Variation of corresponding overpressure.

Figure 14 |
Figure 14 | (a) Variation of flow during valve closure based on the optimal laws that generate the minimum overpressure.(b) Variation of flow during valve closure based on the optimal laws that generate minimum depression.

Figure 13
Figure 13(a) shows the variation, under optimum depression protection conditions, of depression for different valve slow closing times.It can be seen that as the slow closing time increases, minimum depression decreases and tends to steady-state pressure for times greater than 10t 4 .The corresponding overpressures are shown in Figure 13(b).For each slow closing time (t .0.5t 4 ), two optimal laws are derived.The first law generates the minimum overpressure and controls the flow evolution during closure, while the second law generates the minimum depression and also governs the flow evolution during closure.Figure14(a) displays the variation of flow during valve closure for closing times ranging from 1t 4 to 9t 4 , based on the optimal laws that generate the minimum overpressure.Figure14(b), however, shows the flow evolution during closure using the optimal laws that generate the minimum depression.

Table 1 |
Information regarding the proposed problem and used method

Table 2 |
Simulation result in steady state