DOI: 10.5545/sv-jme.2025.1340 337 © The Authors. CC BY-SA 4.0 Int. Licencee: SV-JME Strojniški vestnik - Journal of Mechanical Engineering ▪ VOL 71 ▪ NO 9-10 ▪ Y 2025 Comparison of 1D Euler Equation Based and 3D Navier-Stokes Simulation Methods for Water Hammer Phenomena Nejc V ovk − Jure Ravnik Faculty of Mechanical Engineering, University of Maribor, Slovenia jure.ravnik@um.si Abstract Water hammer phenomena in pipelines can induce significant transient pressure surges, leading to structural failures and operational inefficiencies. This study presents a comparative analyzis of two numerical approaches for simulating water hammer: a one-dimensional (1D) inviscid model with added friction based on the Euler equations and the method of characteristics, and a three-dimensional (3D) viscous model utilizing the Navier-Stokes equations in OpenFOAM. Benchmarking problems are solved first, then both methods are used to study a 3.4 km long DN400 pipeline subject to sudden pump failure by analyzing pressure surges, cavitation, and water column separation. The 1D model effectively predicts transient pressure waves and cavitation conditions with minimal computational cost, while the 3D model provides a detailed representation of multiphase flow dynamics, including cavitation bubble growth and collapse via the volume of fluid method. To mitigate adverse effects, a dynamic combination air valve is introduced, and its effectiveness in reducing pressure surges and cavitation is demonstrated. The results highlight the trade-offs between computational efficiency and accuracy in modelling water hammer events and underscore the importance of protective measures in pipeline systems. Keywords water hammer, cavitation, water column separation, CFD, Euler equation, Navier-Stokes equations, OpenFOAM, method of characteristics Highlights ▪ A 1D inviscid and 3D viscous simulation models were developed for water hammer simulations. ▪ The developed models were compared and used on an example of a 3.4 km long pipeline undergoing sudden pump failure. ▪ Advantages and disadvantages of inviscid versus viscous modelling are discussed. ▪ Results of simulation of cavitation bubble growth on a pipeline with and without a dynamic combination air valve are presented and compared. 1 INTRODUCTION The phenomenon of water hammer in pipelines, particularly during sudden flow blockages, has received significant attention in the field of hydraulic engineering. Water hammer is characterized by transient pressure surges that occur when the flow of fluid is abruptly stopped or altered, often leading to severe mechanical stress on pipeline systems. Traditionally, water hammer analyzis has relied on one- dimensional (1D) inviscid models based on the Euler equations with additional consideration of steady or unsteady friction and the method of characteristics [1]. These models provide a simplified representation of fluid dynamics, allowing for efficient simulations of pressure transients in pipelines and can include cavitation phenomena [2]. More recently, three-dimensional (3D) viscous models based on the Navier-Stokes equations have been employed to provide a more detailed representation of fluid behavior, including the effects of viscosity and turbulence. These models can simulate complex interactions between fluid phases, such as cavitation bubble dynamics and water column separation. The choice between 1D and 3D models often depends on the specific requirements of the analyzis, including computational resources, desired accuracy, and the complexity of the system being studied. While 1D models are computationally efficient and suitable for preliminary assessments, 3D models offer a more comprehensive understanding of fluid dynamics in complex systems. Numerical simulations of water hammer have been widely studied using computational fluid dynamics (CFD). Cao et al. [3] analyzed transient flow in pipelines, emphasizing its importance in urban water systems and hydropower. Khan et al. [4] investigated hydropower penstocks, showcasing CFD’s role in modelling water hammer, cavitation, and column separation. They performed transient CFD simulations for different load rejection conditions using Ansys CFX by modifying the URANS equations. The impact of pipe material properties on water hammer dynamics has been studied extensively. Morvarid et al. [5] analyzed viscoelastic pipe wall effects on pressure fluctuations using the method of characteristics and turbulence modelling. Protective systems like hydropneumatic tanks were investigated by El-Hazek and Halawa [6], showing their effectiveness in damping pressure surges. Air entrainment has been explored as a mitigation strategy. Zhang et al. [7] demonstrated that air pockets can absorb pressure surges in gravitational pipe flows. Additionally, Meng et al. [8] highlighted the influence of flow velocity and pipe wall roughness, finding that higher velocities and roughness exacerbate pressure surges, underscoring their importance in pipeline design. Nikpour et al. [9] emphasized the role of CFD in understanding cavitation and its link to water hammer, crucial for preventing failures in hydraulic systems. They used Ansys Fluent and have shown CFD can be successfully employed in modelling of water hammer phenomena. Ansys Fluent was also used by Han et al. [10] to show that rapid valve closures amplify water hammer pressures, highlighting the need for controlled valve operations. Zhang et al. [11] proposed a dynamic mesh simulation method to analyze transient behavior in pipelines with moving isolation devices, aiding in the design of resilient systems. Aguinaga et al. [12] proposed a mechatronic approach to control water hammer, integrating mechanical, electrical, Process and Thermal Engineering 338 ▪ SV-JME ▪ VOL 71 ▪ NO 9-10 ▪ Y 2025 and hydraulic systems for better transient pressure management. Wu et al. [13] reviewed transient flow percussion theory, emphasizing its role in preventing water hammer in long-distance pipelines. Yang et al. [14] validated 3D CFD simulations as effective tools for analyzing valve-induced water hammer and its impact on pipeline integrity. The interaction between water hammer waves and centrifugal pumps has also been a subject of investigation. Zhang et al. [15] explored the dynamic interactions between valve-closure water hammer waves and pump components, revealing that these interactions can lead to substantial pressure variations and fluid- induced forces on the pump. Malesinska et al. [16] analyzed the effects of sudden cross-section changes on water hammer, showing that abrupt geometry variations significantly influence transient pressure waves. Lupa et al. [17] reviewed water hammer impacts on hydraulic systems, highlighting the importance of empirical validation for simulation reliability. In the present work we focus on the comparison of two methods for assessing flow conditions in a pipeline after a sudden discharge decrease: the 1D inviscid simulation with added friction model and the 3D viscous simulation. In the following subsections we present both methods and then apply them to the analyzis of a 3.4 km long pipeline, which experiences a sudden drop of discharge due to the failure of the electrical grid to deliver power to the pumping station. The 1D solver is based on the method of characteristics and is capable of simulating water hammer in pipelines with various boundary conditions. The 3D model is implemented into OpenFOAM, [18], 3D Navier-Stokes solver, and supports cavitation and multiphase flow. By comparing the results we are able to identify the strengths and weaknesses of both methods. 2 METHODS 2.1 1D Solver When neglecting viscosity the Euler equation describes the momentum balance in the fluid system and at the same time the mass balance is described by the continuity equation. When taking the pipe deformation into account the mass conservation equation reads:       p t c v x  2 0. (1) Here the pressure wave speed is: c D Ee         1 1   , (2) with E the Young’s modulus, e the wall thickness, D the pipe diameter and χ the compressibility. The law of conservation of momentum is established using the force balance and reads as:       v t p x fv D vg 12 2   sgn( )s in , (3) where sgn(v) is the sign of the velocity. The effect of viscosity is modelled with the Fanning friction coefficient f and the wall stress τ w . The wall stress is calculated using the quadratic law of resistance as  w fv  1 2 2 . 2.1.1 Method of Characteristics We solve the coupled Eqs. (1) and (3) using the method of characteristics, [19]. The method of characteristics is a numerical method for solving hyperbolic partial differential equations. It is based on the observation that the solution of a first-order partial differential equation can be represented as a family of curves (characteristics) in the domain of the independent variables. The method of characteristics is well suited for solving these equations as it can capture the shock waves and other discontinuities that arise in the flow. When we combine the two equations with the Lagrange multiplier λ = ±ρ c we obtain a system of two equations named C + and C − : C vv t g c HH t fv D v ij ij ij ij ij ij          :s gn( ) ,, ,, , , 11 11 1 2 1 2    0, (4) C vv t g c HH t fv D v ij ij ij ij ij ij          :s gn( ) ,, ,, , , 11 11 1 2 1 2    0. (5) The equations are discretized in space (∆x) and time (∆t), and the solution is obtained by iterating through the grid points, where the time and position step are connected via ∆x = c∆t. The unknowns in these equations are the flow velocity and piezometric head at location i at time j+1: v i,j+1 and H i,j+1 . Piezometric head is calculated as H = p/ρg+z. Fig. 1 shows the characteristics C + and C − and the intersections where we can calculate the values of the unknown fields. Index i denotes location along the pipe, index j denotes time. Fig. 1. The characteristics C + and C − Boundary conditions can be either known values of head or discharge. For example, if the discharge is known on the left side (i = 0), the boundary condition for head can be calculated from the Eq. (5) and reads: vv j 01 0 , ,   HH c g vv ft D vv jj jj jj 01 10 11 1 2 1 2 ,, ,, ,, sgn( ).           At the other side of the pipe (i = N) the head is known. If there is an open reservoir a there, then the head is equation to the elevation, and the discharge is calculated from the Eq. (4): Hz Nj N    11 1 , , vv g c HH ft D vv Nj Nj Nj Nj Nj Nj     11 1 2 2 ,, ,, ,, () sgn( ).  The friction factor is calculated from the steady state discharge. Assuming known pipe length, the elevation difference between inlet and outlet and the discharge of the friction factor can be calculated from Eq. (4). We developed the 1D model primarily to discover if conditions, which would enable cavitation are present in the pipeline. Thus, when pressure drops below the vapor pressure, we assume that cavitation occurs. If the simulation continues beyond this time instant, the cavitation is not modelled, but rather only the pressure is limited to SV-JME ▪ VOL 71 ▪ NO 9-10 ▪ Y 2025 ▪ 339 Process and Thermal Engineering the vapor pressure. A detailed model of cavitation is implemented in the 3D model. 2.1.2 1D Model Validation To validate the 1D solver, we compared the results to the experimental data from the literature. The researchers [20,21] performed an experiment by creating a pressure wave with a rapid valve closure. A 37.23 m long copper pipe (d = 22.1 mm) connected to two reservoirs was used. The pipe was installed so that it rises by 2.03 m in the flow direction. The fluid flowed through the pipe at a speed of 0.3 m/s. The static pressure in the upstream reservoir was h = 32 m. A ball valve was installed at the end of the pipe, which closed in 0.009 s with the help of a torsion spring. The propagation speed of the pressure waves is given as 1319 m/s. We used a time step of ∆t = 10 −5 s and a spatial step of ∆x = 1.319 cm using 2823 nodes. The friction factor was set to f = 0.009. In Figure 2 we compare the results of the simulations with our method and the experimentally measured values [21] and numerical simulation with steady friction factor [20]. We find a good agreement, especially when the pressure wave arrives first at the measurement point, as the error in pressure surge is less than 2 %. This is the most important part for further calculations, as the highest overpressure and the longest lasting under-pressure are measured at the first pressure wave in the pipeline. We notice that later in the simulation, the simulated wavefronts are sharper in our result compared to the experimental data. This is due to the use of steady state friction factor, which does not account for the transient nature of flow. Our results compare well with the numerical results of Wan et al. [20], who also used a steady state friction factor. The use of an unsteady friction factor would improve the results, but this is not the focus of our study. Fig. 2. Comparison of the temporal development of the static head in the center of the pipeline where time zero corresponds to the moment when the valve starts to close 2.2 3D Navier-Stokes Solver For the 3D pressure surge calculations, we used the open-source software package OpenFOAM v11 [18], which allows simulations of multiphase fluid flows and includes models for cavitation. An analyzis of the numerical results was carried out with the open-source software package ParaView 5.12 [22]. To simulate multiphase flow, we employed the volume of fluid (VOF) method [23–28], which is a numerical technique for capturing the interface between two immiscible fluids. The VOF method models the interface by solving a single set of Navier-Stokes equations for the mixture and introduces a volume fraction field to track the distribution of each fluid within the computational domain. The volume fraction field is a scalar field that represents the fraction of each phase in a given cell. The continuity equation is:       m m t u 0, (6) where u is the flow velocity field and ρ m the mixture density, calculated via the mixing rule from the liquid phase (index l) and gas phase (index v) partial mass densities:    mv l    1, (7) Here α is the gas phase volume fraction. Momentum conservation is described by the Navier-Stokes equations:           m m t pf u uu T , (8) where p is the pressure and f σ = σ κ ∇α the source of momentum due to surface tension between the gas and liquid phase, where σ = 0.07 N/m is the surface tension coefficient for water and water vapor and κ is the interface curvature [18]. Tensor T is the deviatoric part of Cauchy stress tensor, which includes viscosity calculated using the mixing rule, as well as the Reynolds stresses, arising from turbulence, and need to be modelled. To close the system of equations, we use the Menter’s kOmegaSST turbulence model [29]. In the past, it has been discussed, that standard two-equation turbulence models tend to overpredict the eddy viscosity in vapor-liquid mixture zones, suppressing the natural unsteadiness of cavitation [30,31]. This has been solved by Reboud et al. [32], who introduced a correction term for eddy viscosity, improving the modelling of phenomena such as periodic vapor cloud shedding in turbulent cavitation flows [33]. Since our study focuses on pressure-wave propagation and column separation, rather than the detailed structure or dynamics of cavitation, we have not applied the Reboud correction. Moreover, the occurrence of periodic cavitation phenomena, such as cavity shedding, is not expected under the conditions considered in our simulations. The V oF method requires solving an additional equation for the volume fraction field α:        t SS u , (9) where S + and S− are source due to evaporation and condensation. To model turbulence, we chose the Menter kOmegaSST [29] turbulence model. Finally, we solve the energy equation, which includes the effects of phase change in the S T term:           T t E tc E c T k c T mk m vm mk vm m m mp m    11 ,, , uu            () , ,,  m vm m vm T p cc S uu g 11 (10) where E k is the kinetic energy calculated as E k = 1/2|u| 2 and c p,m , c v,m the specific heats. The first term on the right-hand side includes the thermal conductivity of the mixture, k m , that incorporates the molecular thermal conductivity as well as the turbulent thermal conductivity. We model cavitation using the Schnerr-Sauer et. al. [34] model, by modifying the source terms S + and S− in (9) as: SC R pp v b v l       31 2 3 0   max, , (11) Process and Thermal Engineering 340 ▪ SV-JME ▪ VOL 71 ▪ NO 9-10 ▪ Y 2025 SC R pp c b v l       31 2 3 0   max, . (12) Here R b is the cavitation cloud diameter and p v is the vapor pressure. At last, the relation between pressure and density was computed using linear compressibility χ as χ m = αχ v + (1−α)χ l , (13) and dρ m = χ m dp (14) serves as the equation of state. In short, the comparison of the considered phenomena between 1D and 3D simulations goes as follows: • In 1D, no velocity profile develops, since velocity has one compo- nent that points downstream of the pipe. In 3D, we account for the viscosity, which, along with the no-slip boundary condition on the wall, develops a velocity profile. • In 1D, we solve for the wave propagation and do not account for cavitation that might occur as a result of sudden depressuriza- tion. In 3D, we model cavitation through extra terms in the energy equation. • In the 3D simulations we do not account for the deformation of the pipe walls. 2.2.1 3D Model Validation Wang et. al. [35] performed an experiment to investigate the water hammer phenomenon in a pipeline with a sudden valve closure. The pipeline is 5.692 m long, DN 40 and made of plexiglass. The Darcy-Weisbach friction factor of the system ranged between 0.034 and 0.055 (Fanning factor 0.0085 to 0.0138). In their study, multiple scenarios were investigated by varying the static pressures in Tank 1 and Tank 2. This approach allowed them to achieve different initial flow velocities corresponding to different static heads in the pipeline. For the purposes of comparison, we selected the case with an initial velocity of 1.148 m/s and a static head of 1.55 m in Tank 2. After closure, they measured the pressure head versus time. For the remaining details of the experimental apparatus, the reader is directed to the reference [35]. We recreated the experiment numerically using a 2D axisymmetric approach and compared different mesh densities (Fig. 3) and time steps (Fig. 4). We observe good convergence with both mesh density and time step, and the results are in agreement with the experimental data. The error in the prediction of the maximum pressure between simulation and experiment amounts to around 2 % for the coarse and medium meshes, and 1.5 % for the fine mesh. The time step analyzis shows that the solution does not change significantly when decreasing the time step value tenfold. When the time step is decreased by a factor of 100, numerical instability is observed. We attribute this to the fact that for very small time steps approaching machine double- precision limits, the transport phenomena become dominated by the accumulation term. This term can reach disproportionately large values due to the combination of very small time increments and the accumulation of machine precision errors, which are of a similar magnitude to the time step itself. The experimental pressure history, however, exhibits a substantially broader pulse following cavity collapse than predicted in our rigid-wall simulations. The duration of the experimental pulse is governed by the round-trip wave travel time 2L/c eff , where the effective wave speed c eff depends not only on fluid compressibility but also on the compliance of the pipeline and reservoir walls. These elastic effects are absent from our rigid-wall model and therefore shorten the simulated pulse relative to the experiment. In addition, the experimental traces display damped oscillatory tails after the main impulse. The fact that Wang et al. observed the same repeatable waveform across all BV2 valve closure cases strongly supports the interpretation that these features are systematic instrumentation effects rather than random measurement noise. Further broadening of the measured signal could also arise from the influence of elbows and fittings, which extend the effective propagation path and introduce partial reflections and scattering. Taken together, these structural and instrumental effects smear and lengthen the measured time history compared with the rigid-wall model used here. Fig. 3. Pressure head versus time for different mesh densities Fig. 4. Pressure head versus time for different time steps, for mesh with 7500 cells 2.3 Pipeline Length The length of the pipeline is an important parameter in the water hammer analyzis. The longer the pipeline, the longer the time it takes for the pressure wave to travel from the valve to the end of the pipeline. This can result in higher pressure surges and longer duration of underpressure. While for 1D simulations the pipeline length is easily adjusted, for 3D simulations the computational cost increases with the length of the pipeline. Time step analyzis in the previous section showed that good results are achieved when the pressure wave does not travel more than one element within one time step, i.e. time step is limited by the Courant-Friedrichs-Lewy (CFL) condition. This sets the limit for the pipeline length in 3D simulations due to the computational effort required. 3 RESULTS 3.1 The Pipeline To test and compare the 1D and 3D approaches we simulate a pipeline with a length of L = 3408.45 m and a diameter of DN400, SV-JME ▪ VOL 71 ▪ NO 9-10 ▪ Y 2025 ▪ 341 Process and Thermal Engineering which connects a pumping station and a reservoir and assumes an electrical power failure, which stops the pump. The pipeline is made of 17 steel segments; its profile is shown in Fig. 5. The inner diameter of the pipe is 400.1 mm, the outer 406.4 mm, the wall thickness is 6.3 mm. Height difference between the pumping station and the outlet is z 0 = 8.16 m. The Young’s modulus of E = 207·109 Pa is used. In normal operation we consider water (ρ = 999.84 kg/m 3 , χ = 4.54·10 −10 Pa −1 , ν = 1.005·10 −6 m 2 /s) flowing at a rate of Q 0 = 750 m 3 /h with an average velocity of v 0 = 1.66 m/s. The pressure wave speed in these conditions, Eq. (2), is c = 1310 m/s, which gives a characteristic wave travel time of τ = L/c = 2.6 s. The friction factor is calculated for each pipe segment separately; the average is f ave = 0.0067±2.79·10 −5 . Fig. 5. A 17-segment pipeline profile with length of 3408.45 m We assume that electrical power supply fails, which stops the pump. The discharge decreases from Q 0 to Q min = βQ 0 in t stop . The pressure wave travels from the pumping station to the reservoir, where it reflects and travels back. At β = 0 the discharge is Q min = 0 meaning that the power loss completely blocks the flow. This represents the worst-case scenario, as it results in the highest-pressure surges in the pipeline. For 0 < β < 1 the discharge is not completely blocked. We assume that between t 0 and t = t stop linear upstream end discharge variation: Qt QQ t t tt Qt t stop stop stop () () .                       00 0 11 (15) The pump manufacturers estimated the time in which the discharge stops after electrical power failure at t stop is 15 s to 20 s. To estimate the worst-case scenario, we make the following estimate. The pump is rotating at ω 0 = 1488 rpm, has a moment of inertia of I = 3.614 kg/m 2 , its pump efficiency is η = 0.72 and has the electrical power of P el = 132 kW and provides 40.2 m of pressure head. We first estimate the useful work P = ρgΔhQ 0 /η = 14 kW. This gives a normal operation torque of M 0 = P/ω 0 = 732 Nm. The pump stops when the kinetic energy of the rotating parts is converted to the potential energy of the water column. We assume that the average torque is half of the normal operation torque and write a differential equation for the angular acceleration: d dt P I    2 0 . (16) After integration up to time t stop we are able to estimate the time when the pump stops: t I P stop  2 15 0 2  .. s (17) This value serves as the worst-case scenario in the analyzis below. 3.2 1D Simulation Results In this section we present the results obtained using the developed 1D inviscid solver. We focus specifically on the time frame before cavitation occurs as the objective of this work is to identify the pump failure conditions that lead to cavitation. Detailed simulation including cavitation were done with the 3D viscous solver and are presented in the next section. If cavitation does occur in the 1D simulation, we limit the pressure to vapor pressure and let the simulation continue. 3.2.1 Grid Sensitivity Analyzis In Table 1 we show three sets of numerical parameters used in simulations. We compare the results of the simulations with the parameters A, B and C at time t = 3.4τ (Fig. 6). By calculating the relative difference between the head and velocity profiles we obtain the values shown in Table 1. The relative difference norm is calculated as: i i a i b i i a ff f     22 12 /) , / where f is either head or velocity, i is the index of the node and a,b = A, B or a,b = B, C. We observe that the difference in results between numerical parameter sets B and C is very small, which shows that the numerical parameters do not affect the results. We use parameters B in all further simulations. Table 1. Time step, distance between nodes and the number of nodes used for sensitivity analyzis Numerical parameters Time step Δt [s] Δx = cΔt [cm] Number of nodes Norm head Norm velocity A 10 −3 131 2627 B 10 −4 13.1 26047 0.0413 0.0663 C 5·10 −5 6.55 52067 0.0018 0.0029 a) b) Fig. 6. a) Head, and b) velocity profile at t = 3.4 τ Process and Thermal Engineering 342 ▪ SV-JME ▪ VOL 71 ▪ NO 9-10 ▪ Y 2025 a) b) c) d) e) f) g) h) Fig. 7. a), c), e), g) Absolute pressure, and b), d), f), h) flow velocity profiles for three stopping times and four-time instants: a), b) t = 0.5τ = 1.3 s , c), d) t = 1.5τ = 3.9 s ,e), f) t = 2.5τ = 6.5 s and g), h) t = 3.5τ = 9.1 s SV-JME ▪ VOL 71 ▪ NO 9-10 ▪ Y 2025 ▪ 343 Process and Thermal Engineering 3.2.2 Discharge stops completely, β = 0 We simulate the worst-case scenario, when the discharge stops completely, β = 0 and consider three stopping times: t stop = 1.5 s, which represent the worst case scenario, t stop = 15 s, which is the pump manufactures estimate and t stop = 25 s. Figure 7 shows absolute pressure and flow velocity profiles for four-time instances. At half of the characteristic time, the pressure wave has travelled through half of the pipeline. In the case of a fast interruption of the flow, we see that the flow has stopped in the first half of the pipeline and the pressure there has dropped to the vapor pressure. In the case of slow interruption, the flow velocity in the first half of the pipeline only decreases, and the pressure drops, but not yet to the vapor pressure. At t = 1.5τ, the pressure wave has already been reflected from the end of the pipeline and travels back to the pump. In the case of a rapid interruption of the flow, we see that the water in the second half of the pipeline flows towards the pump and the pressure wave consequently also moves towards the pump. When it reaches it at t = 2τ, it will cause a sharp increase in pressure there. If we look at the pressure curve for t stop = 15 s, we can see that the pressure has now also fallen to the vapor pressure in this case. In the case of t stop = 25 s, the pressure is still falling but has not yet reached the vapor pressure. Nevertheless, the vapor pressure also occurs in this case, as can be seen from the pressure curve at t = 2.5τ. At the same time, we see at t = 2.5τ that in the case of t stop = 1.5 s the pressure in the part of the pipeline near the pump has increased significantly, which is a result of the sudden stop of the water flowing towards the pump. This phenomenon is greatly a) b) c) d) e) f) Fig. 8. Time traces of absolute pressure (left) and flow velocity (right) for three flow stopping times: a), b) t stop = 1.5 s, c), d) t stop = 15 s , e), f) t stop = 25 s Process and Thermal Engineering 344 ▪ SV-JME ▪ VOL 71 ▪ NO 9-10 ▪ Y 2025 attenuated in cases where the time for the flow to stop is longer, as it occurs later when the water velocity decreases and then hits the pump at a much lower speed. Figure 8 shows the time histories of absolute pressure and flow velocity at the center of the pipeline and at the beginning of the pipeline near the pump. Green line denotes the values at the pumping station, orange the values at the middle of the pipeline. The results are shown up to 52 s, which is twenty characteristic times. The simulation results are shown for three flow stopping times, for t stop = 1.5 s, t stop = 15 s and t stop = 25 s. We notice considerable differences in the progression. At a very fast stop, t stop = 1.5 s, the pressure wave travels back and forth along the pipeline and at the selected location we observe a clear and distinct velocity fluctuation between positive values (flow from the pump to the end of the pipeline) and negative values (flow back from the end of the pipeline to the pump). The pressure behaves similarly, with the difference that in the part where the pressure drops, it quickly reaches the vapor pressure and the growth of the water vapor column begins. The 1D simulation does not take cavitation into account, so the results in this part show a constant vapor pressure of 2337 Pa. When the flow stop is slower, there is interference between the pressure waves travelling back and forth along the pipeline and between the consequences of the slow decrease in velocity. In both cases, t stop = 15 s and t stop = 25 s the stopping time is longer than the characteristic time τ = 2.6 s. Due to this interaction, we see that up to t stop there is no enormous increase in pressure, but this increase is smaller. Even later, when a more significant increase in pressure occurs, we see that the maximum pressure in the system is much lower than in the case of a rapid flow interruption. a) b) Fig. 9. a) Maximal absolute pressure in the pipeline, and b) the moment when vapor pressure is reached versus flow stopping time t stop for different β In Figure 9 we show the maximum absolute pressure that occurs in the pipeline as a function of the flow stop time tstop. We see that in the case of a very fast flow interruption, when the flow is interrupted before the pressure wave has travelled through the pipeline and back, i.e. when t stop < 2τ = 5.2 s, a very large pressure increase occurs in the system. The increase corresponds to the Joukowsky pressure ∆p = ρcv 0 = 21.7 bar. If the flow is closed slower than 2τ = 5.2 s, the Joukowsky pressure is not reached. At about t stop = 4τ = 10.4 s the maximum pressure drops and reaches the value of normal operation. Interestingly, if the stopping time corresponds to three or four propagation times of the pressure wave through the system 4τ < t stop = 8τ, we see a further, smaller increase in the maximum pressure. This is due to the interference between the propagating wave and the point at which the flow is interrupted. When interpreting this diagram, it must be emphasized that the maximum pressure in the system occurs after the moment when vapor pressure has been reached in the system. Since the numerical model presented in this section limits the pressure to vapor pressure and does not take cavitation into account, the values are too high. A more accurate 3D model, which is presented in the next section, shows that the maximum pressure is lower. At the same time, in Fig. 9, we show the moment when the pressure drops to vapor pressure depending on the time in which the flow is interrupted. We note that regardless of whether the flow is interrupted very quickly (1.5 s) or slowly (45 s), the vapor pressure is always reached. If the discharge reduction is very slow and lasts longer than approximately 45 s, we see that pressure does not decrease to vapor pressure and at the same time it does not increase anywhere in the pipeline. Given that the estimate of the flow interruption time obtained from the pump manufacturer indicates a time of about 20 s, it means that measures are needed in the proposed pipeline to mitigate the consequences of the formation of pressure waves due to the failure of the pump power supply. We see that the case when the flow is completely interrupted (β = 0) is the most demanding case from an engineering point of view, since vapor pressure occurs first in this case and the highest pressure achieved in the system is the highest. Therefore, the results at β = 0 can be considered as the worst possible scenario and if in a real system the pump allows some flow, this alleviates the situation. 3.3 3D Viscous Simulation Results 3.3.1 Reduced Length Pipeline The computational requirements for a 1D simulation of the full length (3.4 km) pipeline for 50 s using a million time steps are about 10 minutes on a single core. The computational requirements for a 3D simulation of the full-length pipeline are much higher, as the number of nodes in the 3D mesh would be at a minimum about 10 6 and the time step is limited by the CFL condition. The computational requirements for a 3D simulation of the full-length pipeline are about 1000 times higher than for a 1D simulation. This means that the computational requirements for a 3D simulation of the full-length pipeline are about 104 hours, which is prohibitive. To reduce the computational cost, we simulate a shorter pipeline (L r = 120 m) with the same diameter (DN400) and the same discharge (Q 0 = 750 m 3 /h) while at the same time at linearly reduced pressure drop. We calculated equivalent flow stopping times for the reduced pipeline so that their ratio to the characteristic time is the same as for the full-length pipeline, i.e. t stop,120 m /(L r /c) = t stop,3400 m /τ. To capture the rapid transient dynamics of column separation and the subsequent water hammer, a time discretization of at least one microsecond was required, resulting in several million time steps for a full-length SV-JME ▪ VOL 71 ▪ NO 9-10 ▪ Y 2025 ▪ 345 Process and Thermal Engineering pipeline simulation. By reducing the numerical pipeline length while preserving equivalent stopping times, the total number of required time steps was decreased, reducing computational time from several weeks to a few days per simulation. The results show that a large cavitation bubble forms behind the pump and due to energy lost for this process the pressure wave reaches smaller absolute pressure values as compared to the inviscid simulation results. Figure 10 shows the pressure curve just downstream of the pump, for different pump stopping times. We find that the magnitude of the pressure surge is largest for stopping time 1.5 s and amounts to 7.4 bar. We tested shorter stopping times as well and found similar maximal pressure surges for them. At a stopping time of 15 s we find that the magnitude of the pressure surge is lower and amounts to 6.3 bar. As the pump stopping time is further extended, the magnitude of the pressure surge decreases. We begin to observe the phenomenon when the pressure wave reflected from the outlet returns to the pump before cavitation occurs. The lower values of the pressure surge can be explained by analyzing the average velocity in the pipeline in Fig. 11. At the moment when the cavitation bubble collapses the flow velocity is lower in the case of longer pump stopping time. A long pump stopping time causes deceleration of the average water velocity in the pipeline, which can be seen in the enlarged image in Fig. 11. The pump stopping time can thus be interpreted as a relaxation rate, which determines the point in time when the cavitation bubble will start to grow. This is evident in Fig. 12, where the cavitation bubble length is shown. The cavitation bubble stops growing when the water velocity in the pipeline is zero. Fig. 10. Pressure time traces for different flow stopping times Fig. 11. Average flow speed for different flow stopping times Figure 13 shows the maximal length of the cavitation bubble and the maximal absolute pressure for different pump stopping times. We see that the maximal length of the cavitation bubble is largest for the shortest pump stopping time and amounts to ≈ 0.55 m. The maximal length of the cavitation bubble decreases with increasing pump stopping time. The maximal absolute pressure is largest for the shortest pump stopping time and amounts to 7.4 bar. The maximal absolute pressure also decreases with increasing pump stopping time. The maximal absolute pressure is lower than in the inviscid simulation, which is due to the energy lost for the cavitation process. Fig. 12. Length of cavitation bubble for different pump stop times Fig. 13. Maximal length of cavitation bubble and maximal absolute pressure for different pump stopping times 3.3.2 Dynamic Combination Air Valve Simulation Both 1D and 3D simulations showed that in the case of sudden loss of electrical power supplying the pump a pressure surge occurs, vapor pressure is reached in the pipeline and cavitation occurs. One possible measure that could be taken to mitigate the consequences of the pressure surge and cavitation, is to install a dynamic combination air valve behind the pump outlet. The dynamic combination air valve opens when the pressure exceeds or is below a certain value and allows the water to flow out or air to be sucked into the pipeline. We simulate the pipeline with the dynamic combination air valve installed and compare the results with the case when the dynamic combination air valve is not installed. We study the use of ARI D-070 dynamic combination air valve [36] for which discharge versus pressure drop curves are available from the manufacturer. The simulation domain is only the first 24 m of the Process and Thermal Engineering 346 ▪ SV-JME ▪ VOL 71 ▪ NO 9-10 ▪ Y 2025 pipeline with the dynamic combination air valve installed, details are shown in Fig. 14. Fig. 14. Geometry of the 3D simulation domain with the dynamic combination air valve installed The boundary conditions used were defined as follows. The discharge at the inlet linearly had the turbulent flow profile shape and the discharge decreased to zero in the chosen flow stopping time tstop. No change in the pressure and the water vapor volume fraction was assumed to occur across the inlet surface. No slip boundary condition was used on the pipe wall. Wall functions were used for turbulent quantities. The flow velocity at the dynamic combination air valve is modelled based on the difference between the pressure in the pipeline and the atmospheric pressure: u K pp pp ms pp y vent atm a atm atm | , / , ,            0 (18) where K = 0.7 was calibrated based on the manufacturer’s pressure drop curve and ρ a = 1.225 kg/m 3 is the air density. The fluid volume fraction depends on the flow direction. When the pressure inside the pipeline exceeds the external pressure, the medium is water; when the external pressure is higher, air enters instead. The negative sign in front of K arises from the chosen coordinate system (see Fig. 14), where the inflow of air is defined in the negative y-direction. At this boundary, a zero pressure gradient is imposed. At the outlet, the pressure is fixed at p = 497,000 Pa, and outflow boundary conditions are applied to the velocity. In Figure 15 we show a comparison of pipelines with and without a dynamic combination air valve. The left panel shows a comparison of the growth of the cavitation bubble. The right panel shows the entry of air through the dynamic combination air valve due to the pressure drop in the pipeline. The results are shown for an equivalent pump stop time t stop = 1.5 s. When the pump stops, we observe the growth of the cavitation bubble in both cases. The ingress of air into the pipeline leads to an increase in the average pressure, so that the cavitation bubble collapses earlier in the case of using a dynamic combination air valve. The results of the pressure surge for the cases with and without a dynamic combination air valve are shown in Fig. 16. The use of a dynamic combination air valve reduces the size of the pressure surge by a third. The reason for the faster collapse of the cavitation bubble and thus the smaller size of the pressure surge can be found in the analyzis of the pressure conditions downstream of the pump before the pressure surge is triggered, in Fig. 17. Here we can see the increase in pressure in the system, which is due to the ingress of air from the environment. The pressure rises above the saturation pressure, which prevents further growth of the cavitation bubble and causes its premature collapse. Fig. 15. Comparison of water vapor content in the pipeline with and without the dynamic combination air valve Fig. 16. Pressure surge in the pipeline with and without the dynamic combination air valve Fig. 17. Pressure at the start of the pipeline, directly after the pump 4 CONCLUSIONS We have developed and compared a 1D inviscid and a 3D viscous numerical simulation tool to model the pressure surge and pressure drop in a pipeline subjected to a sudden suspension of flow. The main advantage of the inviscid simulation is that it can be performed with minimal computational resources for pipeline systems at an engineering level. It can correctly predict the pressure surge and the presence of the conditions that would lead to pressure dropping to vapor pressure. However, it does not model cavitation dynamics. On the other hand, viscous 3D simulation is severely limited by SV-JME ▪ VOL 71 ▪ NO 9-10 ▪ Y 2025 ▪ 347 Process and Thermal Engineering computational resources, making it impossible to simulate a pipeline over its entire length. We have shown that with the necessary steps it is possible to obtain good results for an equivalent short pipeline. The main advantage of viscous simulations is the fact that it is possible to model water phase change and the influence of these changes on the flow dynamics in a detailed 3D model. This advantage becomes particularly clear when the case of simulation of the pipeline with a dynamic combination air valve. Such a simulation is not possible with an inviscid simulation but gives an important insight into the efficiency of such a valve as protection against pressure surges or cavitation. Our simulations have shown that such a measure significantly alters the flow dynamics in the pipeline. The valve releases air into the system when the pressure drops, which prevents the formation of a vacuum and consequently prevents the pipe from flattening. At the same time, the valve allows air and water to be released from the pipeline at the moment when the pressure wave returns and the pressure rises. This reduces the maximum pressure reached in the system and protects the pipeline from rupture. References [1] Kjerrumgaard Jensen, R., Kær Larsen, J., Lindgren Lassen, K., Mandø, M., and Andreasen, A. Implementation and validation of a free open source 1D water hammer code. Fluids 3 64 (2018) DOI:10.3390/fluids3030064. [2] Sadafi, M., Riasi, A., Nourbakhsh, S. A. Cavitating flow during water hammer using a generalized interface vaporous cavitation model. J Fluids Struct 34 190-201 (2012) DOI:10.1016/j.jfluidstructs.2012.05.014. [3] Cao, Y., Zhou, L., Ou, C., Fang, H., Liu, D. 3D CFD simulation and analyzis of transient flow in a water pipeline. J Water Supply Res Technol-Aqua 71 751-767 (2022) DOI:10.2166/aqua.2022.023. [4] Khan, F., Kumar, A., Abbas, A., Gupta, D.P., Binjolkar, P. Computational fluid dynamics based transient investigation of the penstock in a hydropower plant. IOP Conf Ser-Earth Environ Sci 1411 012059 (2024) DOI:10.1088/1755- 1315/1411/1/012059. [5] Morvarid, M., Rezghi, A., Riasi, A., Haghighi Yazdi, M. 3D numerical simulation of laminar water hammer considering pipe wall viscoelasticity and the arbitrary Lagrangian-Eulerian method. World J Eng 15 298-305 (2018) DOI:10.1108/WJE- 08-2017-0236. [6] El-Hazek, A., Halawa, M. Optimum hydro pneumatic tank sizing to protect transmission pipelines supply system against water hammer. Eng Res J (Shoubra) 53 212-221 (2024) DOI:10.21608/erjsh.2023.241739.1228. [7] Zhang, B., Wan, W., Shi, M. Experimental and numerical simulation of water hammer in gravitational pipe flow with continuous air entrainment. Water 10 928 (2018) DOI:10.3390/w10070928. [8] Meng, Q.Z., Yang, C., Liu, L.S., Wang, Z. Effect of wall roughness and flow velocity on water hammer. Adv Mater Res 1014 185-191 (2014) DOI:10.4028/www. scientific.net/AMR.1014.185. [9] Nikpour, M.R., Nazemi, A.H., Dalir, A.H., Shoja, F., Varjavand, P. Experimental and numerical simulation of water hammer. Arab J Sci Eng 39 2669-2375 (2014) DOI:10.1007/s13369-013-0942-1. [10] Han, Y., Shi, W., Xu, H., Wang, J., Zhou, L. Effects of closing times and laws on water hammer in a ball valve pipeline. Water 14 1497 (2022) DOI:10.3390/ w14091497. [11] Zhang, K., Zeng, W., Simpson, A.R., Zhang, S., Wang, C. Water hammer simulation method in pressurized pipeline with a moving isolation device. Water 13 1794 (2021) DOI:10.3390/w13131794. [12] Aguinaga, A., Cando, E., Orquera, E. The Mechatronic approach in the mathematical modelling and simulation to control the water hammer in hydraulic facilities. Proc 9 th World Congr Mech Chem Mater Eng 101 (2023) DOI:10.11159/ icmie23.101. [13] Wu, L., You, R., Wang, W., Zhang, F. Status of research on water hammer effect in long distance pressure pipelines. J Eng Res Rep 23 204-211 (2022) DOI:10.9734/jerr/2022/v23i12778. [14] Yang, S., Wu, D., Lai, Z., Du, T. Three-dimensional computational fluid dynamics simulation of valve-induced water hammer. Proc Inst Mech Eng C-J Mech Eng Sci 231 2263-2274 (2017) DOI:10.1177/0954406216631780. [15] Zhang, W., Yang, S., Wu, D., Xu, Z. Dynamic interaction between valve-closure water hammer wave and centrifugal pump. Proc Inst Mech Eng C-J Mech Eng Sci 235 6767-6781 (2021) DOI:10.1177/09544062211000768. [16] Malesinska, A., Kubrak, M., Rogulski, M., Puntorieri, P., Fiamma, V., Barbaro G. Water hammer simulation in a steel pipeline system with a sudden cross section change. J Fluids Eng 143 091204 (2021) DOI:10.1115/1.4050728. [17] Lupa, S.-I., Gagnon, M., Muntean, S., Abdul-Nour, G. The impact of water hammer on hydraulic power units. Energies 15 1526 (2022) DOI:10.3390/en15041526. [18] Weller, H.G., Tabor, G., Jasak, H., Fureby, C. A tensorial approach to computational continuum mechanics using object-oriented techniques. Comp Phys 12 620-631 (1998) DOI:10.1063/1.168744. [19] Douglas, J.Jr., Russell, T.F. Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures. SIAM J Numer Anal 19 871-885 (1982) DOI:10.1137/0719063. [20] Wan, W. Huang, W. Water hammer simulation of a series pipe system using the MacCormack time marching scheme. Acta Mech 229 3143-3160 (2018) DOI:10.1007/s00707-018-2179-2. [21] Bergant, A., Simpson, R.A., Vìtkovsky, J. Developments in unsteady pipe flow friction modelling. J Hydraul Res 39 (2001) 249-257 DOI:10.1080/00221680109499828. [22] Ayachit, U., Geveci, B., Avila, L. The ParaView guide: Updated for ParaView version 4.3, Technical report (2015). [23] Cazzoli, G., Falfari, S., Bianchi, G.M., Forte, C., Catellani, C. Assessment of the cavitation models implemented in OpenFOAM® under DI-like conditions. Energy Proc 101 638-645 (2016) DOI:10.1016/j.egypro.2016.11.081. [24] Koch, M. Lechner, C., Reuter, F., Köhler, K., Mettin, R., Lauterborn, W. Numerical modeling of laser generated cavitation bubbles with the finite volume and volume of fluid method, using OpenFOAM. Comput Fluids 126 71-90 (2016) DOI:10.1016/j.compfluid.2015.11.008. [25] Ngoc, T.V., Luan, M.N., Hieu, N.K. Numerical modelling of freestream cavitating flow through ship propeller using OpenFOAM. IOP Conf Ser Mater Sci Eng 1109 (2021) 012048 DOI:10.1088/1757-899X/1109/1/012048. [26] Pendar, M.-R., Roohi, E. Investigation of cavitation around 3D hemispherical head-form body and conical cavitators using different turbulence and cavitation models. Ocean Eng 112 287-306 (2016) DOI:10.1016/j.oceaneng.2015.12.010. [27] Zhang, H., Zhang, L. Numerical simulation of cavitating turbulent flow in a high head Francis turbine at part load operation with OpenFOAM. Procedia Eng 31 156-165 (2012) DOI:10.1016/j.proeng.2012.01.1006. [28] Savio, A., Cianferra, M., Armenio, V. Analysis of performance of cavitation models with analytically calculated coefficients. Energies 14 6425 (2021) DOI:10.3390/ en14196425. [29] Menter, F.R. Elements of industrial heat transfer predictions. 16 th Braz Congr Mech Eng (COBEM) Uberlandia (2001). [30] Sorgüven, E., Schnerr, G.H. Modified k - ω model for simulation of cavitating flows. Proc Appl Math Mech 2 386-387 (2003) DOI:10.1002/pamm.200310177. [31] Kevorkijan, L., Pezdevsek, M., Bilus, I., Hren, G. Cavitation erosion modelling - comparison of different driving pressure approaches. Int J Simul Model 22 233- 244 (2023) DOI:10.2507/IJSIMM22-2-639 233. [32] Reboud, J.-L., Stutz, B., Coutier, O. Two phase flow structure of cavitation: Experiment and modeling of unsteady effects. 3 rd Int Symp Cavit, Grenoble 28 1-8 (1998). [33] Peters, A., El Moctar, O. Numerical assessment of cavitation-induced erosion using a multi-scale Euler-Lagrange method. J Fluid Mech 894 A19 (2020) DOI:10.1017/jfm.2020.273. [34] Schnerr, G.H., Sauer, J., Physical and numerical modeling of unsteady cavitation dynamics. Proc 4 th Int Conf Multiphase Flow, New Orleans (2001). [35] Wang, H., Zhou, L., Liu, D. Experimental investigation of column separation using rapid closure of an upstream valve. Appl Sci 13 12874 (2023) DOI:10.3390/ app132312874. [36] Aquestia.com. https://www.arivalves.com/products/water-supply/item/d-070- dynamic-combination-air-valve, accessed 2024-09-01. Acknowledgements  The authors would like to acknowledge the support of their respective institutions and the company Rudis, l.l.c. Trbovlje. Received: 2025-04-08, revised: 2025-07-04, 2025-09-13, accepted: 2025- 09-19 as Original Scientific Paper 1.01. Conflict of interest The authors declare that there is no conflict of interest. Process and Thermal Engineering 348 ▪ SV-JME ▪ VOL 71 ▪ NO 9-10 ▪ Y 2025 Data availability The data that support the findings of this study are available from the corresponding author upon reasonable request. Author contribution Nejc Vovk: Investigation, Writing - review & editing; Jure Ravnik: Methodology, Software, Supervision, Writing - original draft, Writing - review & editing. Primerjava simulacij vodnega udara na podlagi 1D Eulerjeve enačbe in 3D Navier-Stokesove enačbe Povzetek   Poja v  v o dnega  udara  v  ce v o v o dih  lahk o   po vzr o či  po ras t  tlak a,  k ar  vodi do strukturnih okvar cevovodov. Ta študija predstavlja primerjalno analizo dv eh  numeričnih  pris t o po v  za  simulacijo   v o dnega  udara:  eno dimenzionalni  (1D)  ne visk ozni  mo del  z  do danim  trenjem,  ki  t emelji  na  Eulerje vih  enačb ah  in metodi karakteristik, ter tridimenzionalni (3D) viskozni model, ki uporablja N a vier -S t ok eso v e  enačbe  v  OpenF O AM  simulacijsk em  o k o lju.  N ajprej  je  prikazana validacija pristopov, nato obe metodi uporabimo za simulacijo 3,4  km  do lgega  ce v o v o da  DN400,  ki  je  izpo s ta vljen  nenadni  o kv ari  čr palk e,  kjer  anal iziramo   tlačni  udar  s  pre trganjem  v o dnega  s t o plca.  1D  mo del  učink o vit o   napo v eduje  preho dne  tlačne  v alo v e  in  k a vitacijsk e  po go je  z  minimalnimi  računskimi  s tr o ški,  medt em  k o   3D  mo del  zago ta vlja  po dr o bno   študijo  dinamik e  v ečfaznega  t o k a,  vključno   z  ras tjo   in  k o lapsom  k a vitacijskih  mehurje v  po   me t o di  k o nčnih  v o lumno v .  Za  ublažit e v  neželenih  ef ekt o v  je  predlagan  k o mbinirani  zračni  v entil,  za  k at erega  smo  do k azali  učink o vit o s t  pri  zmanjše v anju  tlačnega  udara  in  k a vitacije.  R ezultati  po udarjajo   k o m pr o mise  med  računsk o   učink o vit o s tjo   in  natančno s tjo   pri  mo deliranju  po ja v o v  v o dnega  udara in po udarjajo  po men zaščitnih ukrepo v v ce v o v o dnih sis t emih.  Ključne besede vodni udar, kavitacija, pretrganje vodnega stolpca, CFD,  Eulerje v a  enačba,  N a vier -S t o k eso v e  enačbe,  OpenF O AM,  me t o da  karakteristik.