Insights into CFD Modelling of Water Hammer

: A problem with 1-D water hammer modelling is in the application of accurate unsteady friction. Moreover, investigating the time response of ﬂuid dynamics and unsteady turbulence structures during the water hammer is not possible with a 1-D model. This review article provides a summary of 1-D modelling using the recent ﬁnite volume approach and the discussion extends to a quasi-2-D model and historical developments as well as recent advancements in 3-D CFD simulations of water hammer. The eddy viscosity model is excellent in capturing pressure proﬁles but it is computationally intensive and requires more computational time. This article reviews 3-D CFD simulations with sliding mesh, an immersed solid approach, and dynamic mesh approaches for modelling valve closures. Despite prediction accuracy, a huge computational time and high computer resources are required to execute 3-D ﬂow simulations with advanced valve modelling techniques. Experimental validation shows that a 3-D CFD simulation with a ﬂow rate reduction curve as a boundary condition predicted accurate pressure variation results. Finally, a brief overview of the transient ﬂow turbulence structures for a rapidly accelerated and decelerated pipe ﬂow using DNS (Direct numerical simulation) data sets is presented. Overall, this paper summarises past developments and future scope in the ﬁeld of water hammer modelling using CFD.


Introduction
The water hammer phenomenon is a hydraulic shock condition that occurs in the form of a pressure wave that travels through the pipeline and causes significant damage unless the pressure surge is contained.It mostly occurs in piping systems where there is a sudden change in flow rate caused by operations such as rapid valve opening or closure and pump failure.The pipeline and pipe devices will be damaged as a result of the rate of change of momentum in these operations.Therefore, it is crucial to precisely estimate water hammer pressures using numerical techniques [1][2][3][4].The available numerical methods to study transient pipe flow are graphical methods, analytical methods (available only for very simple boundary conditions), finite difference methods (FDM), finite volume methods (FVM), and so on.The method of characteristics is the most widely used one-dimensional approach.The 1-D finite volume method represents the conservative forms of partial differential equations in the algebraic form [5][6][7][8][9].The FVM is generated by taking into account a small volume that surrounds each nodal point in a geometric mesh and then using the divergence theorem to relate the volume integral to the surface integrals.The surface integrals act as fluxes at the boundaries of each finite volume.
The one-dimensional model underestimates the energy dissipation and the actual decay of the pressure is not captured correctly [10][11][12][13][14].The non-availability of accurate unsteady friction models is a challenge to 1-D water hammer prediction.However, several studies have been employed to develop 1-D pipe unsteady wall-shear stress models.The convolution-based and instantaneous acceleration-based [15] 1-D unsteady friction models were tested only for a limited range of unsteady flow conditions and their performance in general transient flow conditions is unknown.Due to limitations in acquiring precise pressure damping in 1-D simulations, considerable attention has been drawn to using 2-D simulations with accurate turbulence models.Even though these turbulence models were first derived for steady-state flows, it is generally accepted that they may also be valid under quasi-steady conditions.However, the fact is that different turbulent models are based on different presumptions, and the solution is not uniform [1].
Advanced 3-D modelling is essential to understand unsteady friction in transient pipe flows.Computational fluid dynamics (CFD) have been implemented recently for the same because of the rapid advancement of computers.They have been proven to provide precise information about how the flow field changes and have been effectively used to solve unsteady flow problems inside fluid machines [1,16].
This paper gives a general overview of various approaches, especially CFD, used for water hammer modelling.It summarizes the advantages and limitations of various turbulence closure models for water hammer modelling.This article concludes with future research directions.

The Finite Volume Approach
The Finite Volume technique is utilized in solving hyperbolic equations using a Godunov-type algorithm [6][7][8][17][18][19].The first-order Godunov approach assumes the average value of U n i over each cell i at a time level n to be known and the average value of U n+1 i over each cell i at the next time level n is to be found.However, the distribution of U over the cell is unknown.The reconstruction step consists of guessing the distribution of U within the cell i at the time level n.Therefore, the value in each cell is presumed to be constant in this approach.But in the second-order Godunov scheme, a linear value, rather than a constant value, is presumed in each cell.The reconstruction of cell values was utilized in the second-order Godunov scheme.A straight line with a given slope is used to represent the reconstruction in the (Monotonic Upstream-centered Scheme for Conservations Laws (MUSCL)-Hancock scheme).The solution is implemented to calculate the flux through the cell interface and there is a need for a proper limiter function to maintain the cell value within the limit.

Water Hammer Equation with Unsteady Friction
The standard water hammer equation that accounts for friction [20] is as follows: A ρ ∂p ∂x where p is the pressure; Q is the discharge; A is the area of pipe; a is the wave speed; ρ is the fluid density; t is the time; D is the diameter of pipe; J s is the head loss due to the quasi-steady friction, determined by using the Hagen Poiseuille law and Colebrook-White formula [21]; and J u is the head loss attributed to unsteady friction.Zelki [22] introduced a convolution integral approach for one-dimensional models by establishing an exact solution of laminar unsteady friction.The above techniques make use of past local accelerations and weighting functions.The obtained outcomes are tedious and need a lot of storage space on a computer.Hence, the simplified versions of Zielke's approach with less accuracy were proposed by [23][24][25].Vardy and Brown [26][27][28] and Pal et al. [29] extended the convolution integral approach to turbulent flow for smooth as well as rough pipes.This approach yields satisfactory outcomes with less numerical accuracy as the convolution integral is approximated by a small number of weighted coefficients [30,31].Carstens et al. [32] introduced instantaneous acceleration-based methods based on the assumptions that unsteady friction damping results from transient local and convective accelerations.Several alternative formulations [33][34][35] have been proposed since then.One and two-coefficient model models provide satisfactory results for unsteady friction terms.An expression for unsteady friction term for a one-coefficient model can be expressed as follows: where g is the acceleration due to gravity.An expression for unsteady friction term for a two-coefficient model [35,36] can be written as follows: where K ut is the decay coefficients belonging to local accelerations and K ux is the decay coefficient belonging to convective acceleration.

Numerical Solution with Second-Order Finite Volume Method
The first-order fixed-grid method of characteristics approach is widely used in the solution of water hammer equations, which typically involve unsteady friction factors [34,37].The space-line interpolation fixed-grid Method of Characteristics (MOC) approach [38] can provide the variables along C + and C − characteristic lines when the effects of the courant number is considered.P n+1 i and Q n+1 i are calculated by bringing together Equations ( 6) and ( 7) A robust second-order finite volume Godunov-type scheme is utilized for the solution of the water hammer equation with an unsteady friction term.The Riemann problem for a general hyperbolic system is represented as follows: where ; and The pipeline is segmented into N sections using the fixed-grid length ∆x as displayed in Figure 1.The integration of Equation ( 8) for the ith control volume between i − 1/2 and i + 1/2 control interfaces yields where i−1/2 udx; the superscripts n and n + 1 represent the t and t + 1 time levels, respectively; f i+1/2 and f i−1/2 represent the mass and momentum fluxes at the control interfaces i − cell interface [38][39][40][41].The exact solution of the variables at the interface i + 1/2 is obtained after applying Rankine-Hugoniot conditions The following relations are obtained across each wave speed   , after applying the Rankine-Hugoniot condition:  The unsteady friction term is added to the source term ().Source terms () are incorporated into the solution by time splitting utilising a second-order Runge-Kutta discretization to the following explicit methods.
First step: Thus, the mass and momentum fluxes calculated at i + 1/2 for all internal nodes and t t n , t n+1 are as follows: where U n L is the average value of u to the left of interface i + 1/2 at n and U n R is the average value of u to the right of interface i + 1/2 at n.A polynomial reconstruction approach is used to estimate U n L and U n R whose order defines the accuracy of the numerical scheme.In Godunov's first-order accuracy approach, the condition U n L = U n i and U n R = U n i+1 is used.In the second-order Godunov scheme, the MUSCL-Hancock scheme is introduced [38,39,41] The MINMOD limiter is recommended to stop false oscillations, the features of which are found in Toro [41]. Figure 2 shows a flow chart of the second-order Godunov method for water hammer modelling with the MINMOD limiter; however, it can be replaced with any other limiter.

Quasi-2-D Eddy Viscosity Model
During unsteady flow conditions, wall friction affects both the magnitude and shape of pressure waves in the pipes.The fluid flow parameters, like the velocity profile, wall shear stress, and energy dissipation rate, are affected by pipe roughness during transient flow.The accurate prediction of transient waves is required for the design of water pipelines, natural gas pipelines, and pressurized sewage pipeline systems.Pezzinga [42] was the first to present a quasi-2-D turbulence model on hydraulically rough pipes.However, the quasi-2-D model [14,15,42] is computationally intensive and has only been The unsteady friction term is added to the source term S(u).Source terms S(u) are incorporated into the solution by time splitting utilising a second-order Runge-Kutta discretization to the following explicit methods.
First step: Water 2023, 15, 3988 6 of 25 Second step: where S U , and has been linked to previous variables U n i rather than U n+1 i .Third step: where S U , and i is computed by the previous variables U n i .The Courant-Friedrichs-Lewy criterion C r = a∆t ∆x ≤ 1 is already achieved by considering that ∆t is small as compared with ∆x.Although the second-order Runge-kutta method exhibits the stability criterion as well: Hence, the permissible time step for the source term (∆t max,s ) can be provided by the following: Further, the maximum permissible time step with the convective term and the source term is provided by the following: ∆t max,s =min(∆t max,CFL, ∆t max,s ) Since there is unsteady friction in the source term S(u), the water hammer simulations may have a little stability, when ∆t = ∆x a .Hence the proposed stability criterion is C r = a∆t/∆x ≤ 0.9.

Quasi-2-D Eddy Viscosity Model
During unsteady flow conditions, wall friction affects both the magnitude and shape of pressure waves in the pipes.The fluid flow parameters, like the velocity profile, wall shear stress, and energy dissipation rate, are affected by pipe roughness during transient flow.The accurate prediction of transient waves is required for the design of water pipelines, natural gas pipelines, and pressurized sewage pipeline systems.Pezzinga [42] was the first to present a quasi-2-D turbulence model on hydraulically rough pipes.However, the quasi-2-D model [14,15,42] is computationally intensive and has only been used for simple pipe systems and boundary conditions.The model was developed by Silva-Araya and Chaudhry [14] to model unsteady friction during a slow valve closure.The model computes velocity iteratively and the convergence of iterative velocity profiles is subjected to satisfying the continuity equation.The actual friction factor product of the energy dissipation factor and steady-state friction factor at any instant was used in the method of characteristics equations to compute flow variables in the time and space domain.Hanmaiahgari and Maji [43] implemented the model developed by Silva-Araya and Chaudhry [14] to simulate sudden valve closure experiments [21] conducted at different Reynolds numbers.In this model [42,43], a quasi-steady approximation is proposed to compute the frictional losses in unsteady pipe flow.Though this method accurately predicts maximum or minimum pressure without a two-phase flow, it cannot predict the time history of the pressure profile.The quasi-steady approximation underestimates pressure damping as compared with the real system.Hence, the eddy viscosity model is used in the simulation of transient pipe flows [43].The general governing equation for the conservation of mass and momentum [31] is expressed in the following form: A ρ ∂p ∂x where the energy dissipation factor The total shear stress is given as follows: where µ du dr is the shear stress due to viscosity, V is the average velocity, and −ρu v is the Reynolds shear stress.The Reynolds shear stress is estimated from Prandtl's mixing length concept as follows: The Reynolds shear stress is also estimated using the eddy viscosity principle proposed by Boussinesq Both the eddy viscosity and mixing length models match well.The inner region uses the mixing length model and the outer region uses the eddy viscosity model.The distribution of the eddy viscosity is shown in Figure 3.The computation of the mixing length distribution is necessary for the computation of the eddy viscosity.The equation showing the relation between eddy viscosity and mixing length [44] is given as follows: The mixing length equation [44,45] is given as follows: Water 2023, 15, 3988 8 of 25 where λ + 1 is mixing length damping factor for flows with zero pressure gradient [45] λ The mixing length equation [44,45] is given as follows: where  1 + is mixing length damping factor for flows with zero pressure gradient [45]  The more details about the quasi 2D eddy viscosity model is found in Silva-Araya and Chaudhry [14].Hanmaiahgari and Maji [43] evaluated the performance of the eddy viscosity model during sudden valve closure and validated it with the experimental measurements of Adamkowski and Lewandowski [21].They compared the pressure profiles at the downstream end of pipelines for the initial flow velocity  0 = 0.63 m/s and compared these with the experimental data, as shown in Figure 4.The dimensionless pressure and dimensionless time were represented by ( −  0 )/( 0 ) and /, respectively.The more details about the quasi 2D eddy viscosity model is found in Silva-Araya and Chaudhry [14].Hanmaiahgari and Maji [43] evaluated the performance of the eddy viscosity model during sudden valve closure and validated it with the experimental measurements of Adamkowski and Lewandowski [21].They compared the pressure profiles at the downstream end of pipelines for the initial flow velocity V 0 = 0.63 m/s and compared these with the experimental data, as shown in Figure 4.The dimensionless pressure and dimensionless time were represented by (H − H 0 )g/(aV 0 ) and ta/L, respectively.The eddy viscosity model is advantageous as it computes 2D velocity profiles during unsteady flow.The limitation of the eddy viscosity model is that they are computationally intensive and need more computational time compared with convolution-based and instantenous acceleration-based methods [43].They can be used in extended simulations for accurately locating leakages in water distribution systems.

Governing Equations for 3-D CFD Model
The mass and momentum equations for the 3-D CFD water hammer model are the Reynolds-Averaged Navier-Stokes equations represented in tensor notation as follows.
∂ρ ∂t where ρ is the fluid density, u i and u j are the velocity components, t is the time, p is the pressure, and µ is the fluid dynamic viscosity.For turbulent flows, the Reynolds stress tensor ρu i u j is modelled using a turbulence model, but it disappears for laminar flows.It is not required to solve the energy equation during the water hammer calculations since the variation in temperature is insignificant [34,46,47].The compressibility effect is taken into consideration during the simulation, and the relationship is expressed as where K f is the fluid bulk modulus, ρ is water density, and p is absolute pressure.However, this relationship ignores the influence of pipe elasticity, which results in overestimating the wave speed and generating a higher pressure rise than the experimental data.To overcome this issue, the bulk modulus is replaced by [48] where K f is the adjusted bulk modulus; D is the pipe inner diameter; E is the Young modulus of elastic pipe; and e is the pipe wall thickness.The above correction now considers both the fluid's compressibility and the pipe's elasticity.Hence the wave velocity decreases as

Mesh Analysis
Results from 3-D CFD simulations are less accurate and less efficient when the mesh quality is poor.Several mesh-independent analyses with different accuracy standards were conducted in previous studies to determine the most suitable mesh size.The viscous sublayer must also be carefully selected to examine the impact of wall shear stress on other variables.Most of the studies considered the three mesh parameters, as suggested by Martins et al. [49], were λ x = ∆x D measured in the axial direction, λ θ = ∆r θ πD measured in the tangential direction, and λ r = FLT/δ, the ratio of first layer thickness (FLT) to the thickness of the turbulent flow viscous sublayer determined in the radial direction, beginning near the wall and advancing in a core direction.The distribution of grid nodes is concentrated in the proximity to the pipe wall, as shown in Figure 5.
is concentrated in the proximity to the pipe wall, as shown in Figure 5.
Cao et al. [50] conducted a greater number of meshing techniques to analyse the effects of grid density on 3-D CFD simulation results.It was found that both the numerical accuracy and computing cost rise with higher grid densities.The 3-D simulations can achieve high precision and cost-effective computation at some suggested value of grid densities.Higher grid density does not improve accuracy when it exceeds the recommended values; instead, it increases calculations and leads to numerical instability.

Turbulence Modelling
Solving the RANS equation directly can lead to a closure problem because the Reynolds stresses are non-linear and unknown [51].Turbulence models are employed to provide the mathematical closure to equations of motion.They are defined as semiempirical formulations of Reynolds stresses in terms of average velocity gradients.Ghidaoui et al. [2] presented an interesting nondimensional parameter  defined as the ratio of the time scale of the delay in turbulence production to the time scale of the pressure wave.It is given as follows: Cao et al. [50] conducted a greater number of meshing techniques to analyse the effects of grid density on 3-D CFD simulation results.It was found that both the numerical accuracy and computing cost rise with higher grid densities.The 3-D simulations can achieve high precision and cost-effective computation at some suggested value of grid densities.Higher grid density does not improve accuracy when it exceeds the recommended values; instead, it increases calculations and leads to numerical instability.

Turbulence Modelling
Solving the RANS equation directly can lead to a closure problem because the Reynolds stresses are non-linear and unknown [51].Turbulence models are employed to provide the mathematical closure to equations of motion.They are defined as semi-empirical formulations of Reynolds stresses in terms of average velocity gradients.Ghidaoui et al. [2] presented an interesting nondimensional parameter P defined as the ratio of the time scale of the delay in turbulence production to the time scale of the pressure wave.It is given as follows: where D is the inside diameter of the pipe, L is the pipe length, u * is friction velocity dependent on the Darcy-Weisbach friction factor and the mean steady-state velocity, and a is the speed of the wave.Based on the magnitude of the parameter P, three regimes were observed by Ghidaoui et al. [2].When P 1, the steady-state flow's structure and intensity of turbulence are identical to those of transient flow.This indicates that in these types of flows, the steady-state turbulence model remains valid.When P ∼ 1, variations in the structure and intensity of turbulence during water hammer can be seen.There is sufficient time for the turbulence to propagate towards the centre of the pipe, thereby the viability of the steady-state turbulence models should be examined.When P 1, the consequences of the water hammer are negligible.
Of all the turbulence models that the ANSYS FLUENT 2022 R2 (22.2) offers, the k-ε model utilizes less computational resources in modelling flows and is most preferred in engineering applications [13,47,52].The mathematical formulation of the (Re-Normalization Group) RNG k-ε has a similarity with the standard k-ε but with some modifications [53,54].The RNG-based k-ε turbulence model is generated from the instantaneous Navier-Stokes equations utilising a mathematical method called "renormalization group" (RNG) methods.The eddy viscosity approximation is used in the RNG k-ε turbulence model to calculate the Reynolds stress tensor in the momentum equation.Saemi et al. [55] carried out 2-D axisymmetric calculations based on unsteady Reynolds-averaged Navier-Stokes (URANS) equations and a high Reynolds RNG k-ε model.They found that the high Reynolds RNG k-ε model cannot reproduce more accurate wall shear stress distributions and other characteristics of turbulent water hammers close to the wall.
The realizable k-ε model fulfils certain mathematical limitations on the Reynolds stresses that are consistent with turbulent flow physics.The Boussinesq hypothesis is used to establish the relationship between the Reynolds stresses and the mean velocity gradients in the model.The Boussinesq hypothesis is preferred over the Reynolds stress model (RSM) because of the low computational cost involved in the calculation of turbulent viscosity.Previous works on separated flows and flows with complex secondary flow characteristics have demonstrated that the realizable model yields good results compared with all other k-ε models.
where k is the turbulent kinetic energy, υ T is the kinetic eddy viscosity, and S ij is mean velocity gradient, defined as The extra transport equations (k and ε) specific for the realizable k-ε model are expressed as follows: where G k is the turbulence kinetic energy generation due to the mean strain rate; Y M is the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate; G b is the generation of turbulence kinetic energy due to buoyancy; σ k and σ ε are the turbulent Prandtl numbers for k and ε, respectively; µ T = ρC µ k 2 /ε is the turbulent viscosity; C 1 = max 0.43, η η+5 ; η = Sk ε ; S = 2S ij S ij ; and C 2 , C 1ε , and C 3ε are constants.Yang et al. [1], Martins et al. [13,56], Wu et al. [57], and Ferreira et al. [58] employed the realizable k-ε model in the 3-D CFD model for the simulation of water hammer problems.The realizable k-ε model has the drawback of producing non-physical turbulent viscosities in case of multiple reference frames and rotating sliding meshes where the computational domain contains rotating and stationary fluid zones The other drawback of the k-ε model is that it needs to be modified to be appropriate for unsteady turbulent flows with adverse and favourable pressure gradients [59].
A more accurate prediction of the transient pressure during valve closure can be possible with the shear stress transport k-ω model [50,60].The refined mesh and the contribution of the boundary layer both have a considerable impact on the correct estimation of shear stress, which is easily incorporated into the SST k-ω model.The Boussinesq approximation is used in k-ω model to derive the unknown Reynolds stress tensor as where ω is the specific rate of dissipation, υ T is the kinematic eddy viscosity, and δ ij is the Kronecker delta.
The extra transport equations (k and ω) specific for the SST k-ω model are expressed as follows: where

Various 3-D CFD Approaches in Water Hammer Modelling
CFD solvers are set up with suitable boundary conditions which define the physics of fluid flow.At the inlet boundary, a constant pressure is set and at the outlet boundary, a time domain mass outflow is given, as shown in Figure 6.Three different cases of valve closure: instantaneous, linear, and flow rate reduction curves obtained from 3D valve modelling are also presented in Figure 7.
The flow velocity is reduced linearly from a steady state to zero at the outlet to represent the valve closure in the simulation of laminar water hammer flow by Karney et al. [63].There was a good agreement in pressure variations between the experimental data and CFD results as compared with MOC models.Martins et al. [47] simulated water hammer flow using a 3-D URANS and k-ε turbulence model.The study utilised a three-dimensional effective mesh created by Martins et al. [49] and the discharge reduction curve as the outlet boundary condition as defined by Wylie and Streeter [64].The discharge reduction curve only enhanced the slope of the pressure change and no variation was found between the results of 3-D water hammer predictions and 2-D axisymmetric computations utilising such a simple boundary condition.Saemi et al. [46] investigated 2-D axisymmetric simulations for water hammer flow in FLUENT where a discharge reduction curve was used as an outlet boundary condition.The study showed that using the discharge curve generated from 3-D simulations increases the accuracy of 2-D axisymmetric calculations resulting in identical prediction in 2-D and 3-D.In the present work, the results of the 3-D CFD simulation using the flow rate curve as an outlet boundary condition are validated with the experimental data of Bergant et al. [34].The pipe system was made of 37.2 m long copper pipe, with 0.022 m diameter and 1.6 mm wall thickness.The initial steady-state flow rate was taken to be 0.076 × 10 −5 m 3 /s and the upstream reservoir head H 0 as 32 m.A downstream ball valve was closed rapidly in 0.009 s. Figure 8a,b shows the excellent match of the 3-D CFD pressure history at the valve and midpoint of the pipe with its experimental counterpart.
The 3-D valve modelling method is incapable of visualizing the 3-D flow structures generated by the valve closure inside the domain.Advanced valve modelling techniques are necessary to obtain accurate results.Hence, several studies focused on modelling valve closure based on the dynamic mesh technique.Dynamic mesh is often used in simulations of rotating equipment, oscillating structures, fluid-structure interactions, and free-surface flows with changing fluid domains due to motion on domain boundaries.ANSYS FLUENT 2022 R2 (22.2) let users specify domain motion with boundary profiles or user-defined functions.The geometry consists of two domains: Domain 1 represents the portion of pipe upstream and downstream to the valve and Domain 2 represents the valve's space in the pipe.The lower region of the valve is represented by the upper region of domain 2. As the valve moves into the pipe represented by the downward moving arrow in Figure 9, the volume of domain 2 starts reducing.The valve (solid) part is seen between domains 1 and 2 as it moves inside the pipe.Interfaces connecting the domains are updated at every time step.
ANSYS FLUENT 2022 R2 (22.2) is used in this study since it allows both smoothing and re-meshing at every time step.Neyestanaki et al. [65] employed the spring/Laplacian-based smoothing method.This method brings each mesh vertex closer to the geometric centre of its surrounding vertex.Using smoothing for larger displacement causes cell deformation to become greatly skewed.As a result, re-meshing with triangular or tetrahedral mesh is employed when the mesh is more sensitive to mesh changes [65].
Saemi et al. [62,66,67] employed a dynamic mesh approach for modelling axial gate valve closure.The static pressure at the outlet and the total pressure at the inlet are used in the simulation.The limitation of the dynamic mesh technique is that re-meshing is time-consuming and can result in divergence near the end of the closure of the valve when re-meshing is carried out in a smaller region [61].In summary, dynamic mesh is a powerful technique that enables fluid dynamics in complicated systems to be simulated with realistic and greater precision.ANSYS FLUENT 2022 R2 (22.2) is used in this study since it allows both smoothing and re-meshing at every time step.Neyestanaki et al. [65] employed the spring/Laplacianbased smoothing method.This method brings each mesh vertex closer to the geometric centre of its surrounding vertex.Using smoothing for larger displacement causes cell deformation to become greatly skewed.As a result, re-meshing with triangular or tetrahedral mesh is employed when the mesh is more sensitive to mesh changes [65].
Saemi et al. [62,66,67] employed a dynamic mesh approach for modelling axial gate valve closure.The static pressure at the outlet and the total pressure at the inlet are used in the simulation.The limitation of the dynamic mesh technique is that re-meshing is timeconsuming and can result in divergence near the end of the closure of the valve when remeshing is carried out in a smaller region [61].In summary, dynamic mesh is a powerful technique that enables fluid dynamics in complicated systems to be simulated with realistic and greater precision.The sliding mesh approach provides accurate results when transient interaction between fixed and rotating components is desired in the flow simulations.More than one cell zones are utilised in the sliding mesh approach; the mesh is created locally in each zone.Each zone (stationary or moving) is surrounded by one or more "interface zones" that connect to the neighbouring cell zone.The two cell zones move relative to one another along the mesh interface during the calculation.The sliding mesh technique is utilised for modelling the dynamic course of the ball valve closing as shown in Figure 10.The valve body rotates around its centre during the closing operation and the two pairs of spherical surfaces act as interface surfaces.The area of contact between the two interfaces is reduced as the valve body starts rotating, and after a certain angle, the interface surfaces are completely separated.Due to the loss of contact between the interfaces at their intersection, no fluid can move between these two regions, which denotes complete valve closure.The sliding mesh approach provides accurate results when transient interaction between fixed and rotating components is desired in the flow simulations.More than one cell zones are utilised in the sliding mesh approach; the mesh is created locally in each zone.Each zone (stationary or moving) is surrounded by one or more "interface zones" that connect to the neighbouring cell zone.The two cell zones move relative to one another along the mesh interface during the calculation.The sliding mesh technique is utilised for modelling the dynamic course of the ball valve closing as shown in Figure 10.The valve body rotates around its centre during the closing operation and the two pairs of spherical surfaces act as interface surfaces.The area of contact between the two interfaces is reduced as the valve body starts rotating, and after a certain angle, the interface surfaces are completely separated.Due to the loss of contact between the interfaces at their intersection, no fluid can move between these two regions, which denotes complete valve closure.In this approach, the domain used to model the valve body is known as an immersed solid domain.The solid domain is immersed into the fluid region represented by the pipe indicating the valve closure process shown clearly in Figure 11b.In the overlapping of solid and fluid regions, the immersed solid technique introduces a source term in the momentum equation to force the fluid velocity to match the velocity of the immersed solid.Kalantar et al. [61] employed the immersed solid method for modelling the valve closure to study the pressure-time method during the water hammer phenomenon.The study showed that there could be a chance for the leaking of water through the immersed solid as the valve starts closing.Neyestanaki et al. [65] presented and compared various methods for axial valve closure in a 3-D water hammer simulation.The result showed that the immersed solid approach underestimates differential pressure rise and overestimates wall shear stress close the end of valve closing due to a lag in discharge reduction.In the case of axial valve motion, the immersed solid approach is easier and preferable to the dynamic mesh approach because the need for mesh deformation and re-meshing is eliminated in the method [65,69].Yang et al. [1] used a 3-D computational fluid dynamics model with the sliding mesh method for modelling water hammer and compared those results with 1-D MOC calculations based on unsteady friction models.The results also demonstrated that sliding mesh is a reliable technique for simulating the rapid closing of rotating ball valves.Saemi et al. [46] investigated 3-D numerical simulation for laminar and turbulent water hammer flow which utilized a sliding mesh approach for closing the valve.In this study, the axisymmetric model was also used in which the flow rate reduction was used as an outlet boundary condition.However, compared with the axisymmetric simulations, 3-D CFD simulations using sliding mesh provided accurate pressure variation results; however, they required vast computational resources.The sliding mesh technique was utilized to realize the rotation process of the ball valve to analyse the impact of various closing laws and closing times on the water hammer [68].Cao et al. [50] developed a 3-D CFD model to study the 3-D dynamic features of the ball valve closure.The ball valve closure was accomplished by utilising the sliding mesh technique.It was observed that during the closing operation, mean velocity variations and differential pressure at the valve underwent three different stages.The axial motion of a thin gate valve can be modelled as a mesh movement without any deformation using the sliding mesh technique [65].This approach predicts more accurate outcomes that are closer to the experimental measurements.The sliding mesh is the most preferred approach for valve closure modelling in the literature as they are less expensive and present more accurate results.
In this approach, the domain used to model the valve body is known as an immersed solid domain.The solid domain is immersed into the fluid region represented by the pipe indicating the valve closure process shown clearly in Figure 11b.In the overlapping of solid and fluid regions, the immersed solid technique introduces a source term in the momentum equation to force the fluid velocity to match the velocity of the immersed solid.Kalantar et al. [61] employed the immersed solid method for modelling the valve closure to study the pressure-time method during the water hammer phenomenon.The study showed that there could be a chance for the leaking of water through the immersed solid as the valve starts closing.Neyestanaki et al. [65] presented and compared various methods for axial valve closure in a 3-D water hammer simulation.The result showed that the immersed solid approach underestimates differential pressure rise and overestimates wall shear stress close to the end of valve closing due to a lag in discharge reduction.In the case of axial valve motion, the immersed solid approach is easier and preferable to the dynamic mesh approach because the need for mesh deformation and re-meshing is eliminated in the method [65,69].

Turbulence Structures in Accelerating and Decelerating Flows
The accelerating turbulent pipe flow occurs when there is a sudden valve opening, a step increase in flow rate, a pump starting, and load acceptance in turbines.A turbulent flow in a pipe suddenly accelerates subjected to a total acceleration greater than zero, as in Equation (47) where  is the average velocity of the flow.
The decelerating turbulent pipe flow occurs when there is a sudden valve closing, a step decrease in flow rate, a pump tripping, and load rejection in turbines.A steady flow in a pipe is suddenly decelerated if total acceleration is less than zero.

Turbulence Structures in Accelerating and Decelerating Flows
The accelerating turbulent pipe flow occurs when there is a sudden valve opening, a step increase in flow rate, a pump starting, and load acceptance in turbines.A turbulent flow in a pipe suddenly accelerates subjected to a total acceleration greater than zero, as in Equation (47) where V is the average velocity of the flow.
The decelerating turbulent pipe flow occurs when there is a sudden valve closing, a step decrease in flow rate, a pump tripping, and load rejection in turbines.A steady flow in a pipe is suddenly decelerated if total acceleration is less than zero.
Figures 12 and 13 show the velocity profile in the case of accelerating and decelerating flows, respectively.Most fluid conveyance methods through pipes are not steady from the industrial perspective.In addition, control valves or pumps continuously adjust the flow through pipes to satisfy certain demands.Most engineering designs exclude transient effects in unsteady flows because of a lack of understanding about such flows.Hence, to estimate the frictional drag in accelerating, as well as decelerating, pipe flow, not many valid unstable friction models are available.Uncertainty regarding the transitional mechanisms and durations of the various transitional stages may result from the chosen step-down size or the masking impact caused by the flow structures at the start of the transient.The study of turbulence structures in accelerated and decelerated flows is important to an in-depth understanding of unsteady friction.
The comparison between sudden accelerated flows and sudden decelerated flows indicates that the velocity profile is more stable in an accelerating flow than in that of a decelerating flow [70,71].Therefore, there may be a significant variation in how accelerating and decelerating flows behave.Vardy and Brown [71] investigated a specific type of unsteady flow where a steady flow through the pipe is abruptly accelerated.The general flow pattern and the presence of the radially travelling viscosity front, in particular, are the properties of suddenly accelerating flows.The accelerated flow is defined by three different stages.There is a lag in the turbulence response at stage 1, then a strong turbulence response at stage 2, and in stage 3 there exists a persistent lag in the response that resembles a succession of quasi-steady states.Turbulence's axial components show a sudden response, while the corresponding changes in radial and tangential components occur much later.The effects of varying the pipe diameter, the initial Reynolds number, and the acceleration rate on the flow characteristics were investigated by He et al. [70].Their study showed that the most essential pattern of behaviour during stage 1 remained the same across all the cases, despite major variations existing in stages 2 and 3.The turbulence models utilized in the analysis of accelerated and decelerated flows by Vardy and Brown [71] predicted similar qualitative behaviour, but there were major variations in when and how the delayed turbulence response to the acceleration triggers the radially moving front.The effect of a water-hammer wavefront is identical to a sudden shift in axial velocity at a given axial location, and then a slow shift in the radial velocity distribution [71].The two common characteristics regarding the prescribed turbulent viscosity distributions are that they approximate the radial distributions of the effective viscosity and with regard to the selected distributions as frozen.They concluded that only the first of the three flow stages can be approximated by assuming unaffected effective viscosity from the previous steady flow.This stage's predictions are slightly insensitive to the turbulence model; therefore, they can be used to support simplified models.
The behaviour of the mean wall shear stress and its sporadic nature directly affect the energy cost of transportation networks.Using the direct numerical simulation data sets of turbulent flow in a pipe at various Reynolds numbers, Guerrero et al. [72] investigated the highly intermittent nature of wall shear stress and its relation with other flow parameters.The analysis provided a thorough understanding of the flow dynamics responsible for backflow events, their connection to the larger scales of motion in the outer region, and their impact on turbulence production, turbulent inertia, and viscous dissipation.They demonstrated that the presence of two large-scale anticlockwise rotating motions in the buffer zone plays a crucial role in the prevalence of backflow events.Extremely positive streamwise wall shear stress regions are produced by an intense quasi-streamwise vortex which carries streamwise momentum from overlap and outer regions towards the wall.

Accelerating Flows
Guerrero et al. [73] investigated the characteristics of accelerating turbulent pipe flows after a quick increase in flow rate from a more detailed perspective.A high-spatiotemporal DNS data set was employed for this purpose.Maruyama et al. [74] described a delayed response in turbulence in the seminal investigation of a step-wise accelerating turbulent flow through pipes linked to two lags one with new turbulence and the second with the wall-normal diffusion of vorticity caused by the acceleration exerted on an already turbulent base flow.Further, He and Jackson [75] observed a third lag due to the redistribution of TKE among the three orthogonal velocity components.Greenblatt and Moss [76], using singlecomponent LDV, collected time series data sets of quickly accelerating pipe flow.The result revealed that displacement thickness, as well as momentum thickness, showed coherence with the progression of time.Jung and Hung [77] used large eddy simulation data sets to replicate the investigations of He and Jackson [75].These simulation results also supported the three-phase lags in the turbulence response.Using DNS, Guerrero et al. [73] examined the phases that occur in turbulent flow as the flow rate is rapidly increased, as measured by the temporal evolution of dynamic variables in the mean momentum balance.Inertial, pre-transition, transition, and core relaxation are the four distinct phases that characterise the transient behaviour of a turbulent flow after a rapid temporal acceleration [76,78].According to the results, an equilibrium in the skin friction coefficient is not always a reliable predictor of when the flow is to reach its final steady state.By analysing the transient nature of the largest uniform momentum zone of the flow, and the core region during the last transitional phase, this study explains the delayed regeneration of turbulence in the flow's wake zone utilizing the methods proposed by Guala et al. [79].These findings pave the way for the creation of more accurate unsteady friction models and provide a deeper insight into the dynamics of unsteady turbulent pipe flow [27,71,80].The temporal evolution of turbulent eddies, flow dynamics, momentum distribution, and Reynolds and viscous shear stresses is analyzed to determine the relationship between flow dynamics and frictional drag, particularly at inertial and core-relaxation periods [73].Figure 14a shows a linear increment in the Bulk Reynolds number (Re b = 2U b R/ν) on a time domain scaled by the initial viscous units (t +0 = tu 2 τ,0 /ν) and Figure 14b depicts the time response of the skin friction coefficient (C f = 2u 2 τ /U 2 b ) [73].
Water 2023, 15, x FOR PEER REVIEW 20 of 25 core relaxation are the four distinct phases that characterise the transient behaviour of a turbulent flow after a rapid temporal acceleration [76,78].According to the results, an equilibrium in the skin friction coefficient is not always a reliable predictor of when the flow is to reach its final steady state.By analysing the transient nature of the largest uniform momentum zone of the flow, and the core region during the last transitional phase, this study explains the delayed regeneration of turbulence in the flow's wake zone utilizing the methods proposed by Guala et al. [79].These findings pave the way for the creation of more accurate unsteady friction models and provide a deeper insight into the dynamics of unsteady turbulent pipe flow [27,71,80].The temporal evolution of turbulent eddies, flow dynamics, momentum distribution, and Reynolds and viscous shear stresses is analyzed to determine the relationship between flow dynamics and frictional drag, particularly at inertial and core-relaxation periods [73].Figure 14a  Over the past ten years, there has been an emphasis on the temporal response of frictional drag in accelerating turbulent pipe and channel flows [81].Several numerical analyses have been conducted in this area [78,82].Using DNS data sets, He and Seddighi [78] investigated the time dependency of the mean skin friction coefficient   as well as other mean flow statistics to identify the phases encountered by a turbulent channel flow after a quick increase in discharge.According to the findings,   () in a suddenly accelerated turbulent channel flow resembles the bypass transition of a turbulent boundary layer.Also, a suddenly accelerated channel flow goes through two transient phases: pretransition and transition.During the pretransition phase, a perturbation boundary layer forms as a result of a plug-like inflow into the turbulent base flow, whose turbulence remains fixed at first.Hence the skin friction coefficient   declines, reaching its lowest value at a point where the pre-transitional phase ends and the transitional phase Over the past ten years, there has been an emphasis on the temporal response of frictional drag in accelerating turbulent pipe and channel flows [81].Several numerical analyses have been conducted in this area [78,82].Using DNS data sets, He and Seddighi [78] investigated the time dependency of the mean skin friction coefficient C f as well as other mean flow statistics to identify the phases encountered by a turbulent channel flow after a quick increase in discharge.According to the findings, C f (t) in a suddenly accelerated turbulent channel flow resembles the bypass transition of a turbulent boundary layer.Also, a suddenly accelerated channel flow goes through two transient phases: pretransition and transition.During the pretransition phase, a perturbation boundary layer forms as a result of a plug-like inflow into the turbulent base flow, whose turbulence remains fixed at first.Hence the skin friction coefficient C f declines, reaching its lowest value at a point where the pre-transitional phase ends and the transitional phase starts.The skin friction coefficient returns to its original value during the transition phase as a result of the formation of 'new' turbulent areas that expand, merge, and spread across the whole flow domain.The skin friction coefficient C f reaches an equilibrium as the flow is fully developed.Guerrero et al. [73] utilized DNS data sets to analyze the transient flow dynamics of a series of suddenly accelerated turbulent pipe flows between two steady Reynolds numbers (Re).The temporal variation in the mean flow dynamics showed coherence and demonstrated that a suddenly accelerated internal flow goes through four transient phases, namely inertial (phase I), defined by viscous forces that grow rapidly and a turbulent behaviour freezes; pre-transition (phase II), a quick decay of viscous forces and a slow turbulence response near the wall; transition (phase III), viscous and turbulent forces rising proportionally in the inner region; and core relaxation (phase IV), turbulence moving slowly from the wall to the wake zone.

Decelerating Flows
The decelerating flows show 'frozen' turbulence behaviour similar to accelerating flows in the early flow excursion [74,81].Maruyama et al. [74] carried out a similar experimental investigation that shows that the total kinetic energy decays at a roughly linear rate following the initial "frozen turbulence" period.By analyzing the time series of the average frictional velocity u τ (t) and the average centerline velocity U c , a decelerating flow is followed by two transient phases that have rapid and slow time responses [83].However, it does not explain the mechanisms underlying these two distinct relaxation times.More recently, the total kinetic energy production budget has been divided into three parts and their behaviour is analyzed for accelerating and decelerating flows by Joel and Sundstrom [84].They concluded that the turbulence establishment in those two unsteady flows varies significantly.Only a limited number of studies are available on the time variation in the mean skin friction coefficient and the near-wall dynamics related to the time-varying decelerating internal flow, and the majority of them do not examine frictional drag throughout the entire transient phase [81,85].There have been very few attempts to analyze the transient frictional drag in decelerated internal flows using DNS [82,83].
Guerrero et al. [86] investigated the turbulence flow dynamics and time response experienced by sudden decelerating turbulent flow through pipes using direct numerical simulation data sets of linearly decelerated turbulent pipe flows.The turbulence response of a fully developed turbulent flow exhibits a phase lag when the flow is subjected to linear acceleration.It is concluded from the study that (i) there is a phase lag in the turbulence response of a quickly decelerating flow, as is shown in studies of accelerating flows as well; (ii) the presence of large-scale layers of negative azimuthal vorticity near the wall is indicated by the fact that the viscous stress can become negative there; (iii) when a flow excursion occurs, the viscous force's sign switches in the viscous sublayer from negative to positive, and this change propagates in the direction normal to the wall by viscous diffusion; and (iv) the shape of the average velocity profile is changed within the near-wall zone and its linear behaviour is recovered during this time.However, the overlap and wake regions of the flow do not change, which indicates that the outer flow turbulence freezes.(v) Meanwhile, the near-wall region displays a stationary mean velocity profile.Figure 15a shows a ramp down in flow rate was enforced at t +1 = 0 to rapidly reduce the Reynolds number and Figure 15b depicts that the time response of the skin friction coefficient indicates that flow undergoes several transient behaviours before becoming stationary [86].

Future Research Directions
Parallel computing techniques are receiving attention as they reduce computational time and costs.Parallel computing has been increasingly accessible, affordable, and widely applied as the number of available supercomputers has grown.Graphical processing units (GPUs) offer appealing features for large-scale computing [87].Parallel computing with the GPUs together can improve the computational resources required for 3-D CFD studies.
Advancing our understanding of fluid-structure interaction may be possible with the use of 2D and 3D finite volume modelling that can capture turbulence-driven vortices in transient flow It can help in analyzing problems like pipe incrustation, and wall erosion enabling preventive and rehabilitation actions.Hence, the developments in numerically precise models can help in cost savings and prevent issues in pipe design, management, and operational domains.

Conclusions
Designing, operating, and assessing the risks associated with a pipeline system require transient analysis.This article provides historical developments of water hammer modelling, focussing particularly on the recent CFD simulation in improving the transient flow numerical prediction.This study briefly reviewed the application of Godunov's second-order finite volume method (FVM) technique in solving water hammer equations with unsteady friction and the discussion is extended to the application of the quasi-2-D eddy viscosity turbulence closure model.This paper reviewed the pros and cons of the sliding mesh, immersed solid approach, and the dynamic mesh approaches for modelling valve closures.Based on findings from previous literature, this study concludes that the sliding mesh approach is less expensive than other methods and is utilised mostly for valve closure modelling.However, the valve modelling techniques require vast computational resources.This study shows the 3-D CFD simulation with a flow rate reduction curve as a boundary condition accurately predicts the pressure variation close to the experimental results.Finally, a brief overview of the transient flow turbulence

Future Research Directions
Parallel computing techniques are receiving attention as they reduce computational time and costs.Parallel computing has been increasingly accessible, affordable, and widely applied as the number of available supercomputers has grown.Graphical processing units (GPUs) offer appealing features for large-scale computing [87].Parallel computing with the GPUs together can improve the computational resources required for 3-D CFD studies.
Advancing our understanding of fluid-structure interaction may be possible with the use of 2D and 3D finite volume modelling that can capture turbulence-driven vortices in transient flow It can help in analyzing problems like pipe incrustation, and wall erosion enabling preventive and rehabilitation actions.Hence, the developments in numerically precise models can help in cost savings and prevent issues in pipe design, management, and operational domains.

Conclusions
Designing, operating, and assessing the risks associated with a pipeline system require transient analysis.This article provides historical developments of water hammer modelling, focussing particularly on the recent CFD simulation in improving the transient flow numerical prediction.This study briefly reviewed the application of Godunov's second-order finite volume method (FVM) technique in solving water hammer equations with unsteady friction and the discussion is extended to the application of the quasi-2-D eddy viscosity turbulence closure model.This paper reviewed the pros and cons of the sliding mesh, immersed solid approach, and the dynamic mesh approaches for modelling valve closures.Based on findings from previous literature, this study concludes that the sliding mesh approach is less expensive than other methods and is utilised mostly for valve closure modelling.However, the valve modelling techniques require vast computational resources.This study shows the 3-D CFD simulation with a flow rate reduction curve as a boundary condition accurately predicts the pressure variation close to the experimental results.Finally, a brief overview of the transient flow turbulence structures for a rapidly accelerated and decelerated pipe flow using DNS (direct numerical simulation) data sets is /2 −   ) − ( +1/2 −   ) (13) For the upstream boundary,  1/2 () = ( 1/2 ,  1/2 ),  1/2 =  1/2  1/2 ().Similarly, for the downstream boundary,  +1/2 () = ( +1/2 ,  +1/2 ),  +1/2 =  +1/2  +1/2 () attained by combining the Riemann invariant of the equation with a pressure-flow boundary condition.Here,  −1 +1 =  0 +1 =  1/2 and  +1 +1 =  +2 +1 =  +1/2 , where  1/2 and  +1/2 are computed through the combination of the Riemann invariant with a pressure-flow boundary relation at time .

Water 2023 , 25 Figure 4 .
Figure 4. Measured and computed pressure profile for sudden valve closure with initial velocity  0 = 0.63 m/s [43].The eddy viscosity model is advantageous as it computes 2D velocity profiles during unsteady flow.The limitation of the eddy viscosity model is that they are computationally intensive and need more computational time compared with convolution-based and

Figure 4 .
Figure 4. Measured and computed pressure profile for sudden valve closure with initial velocity V 0 = 0.63 m/s [43].

Water 2023 ,
15, x FOR PEER REVIEW 13 of 25 2-D axisymmetric calculations resulting in identical prediction in 2-D and 3-D.In the present work, the results of the 3-D CFD simulation using the flow rate curve as an outlet boundary condition are validated with the experimental data of Bergant et al. [34].The pipe system was made of 37.2 m long copper pipe, with 0.022 m diameter and 1.6 mm wall thickness.The initial steady-state flow rate was taken to be 0.076 × 10 −5 m 3 /s and the upstream reservoir head  0 as 32 m.A downstream ball valve was closed rapidly in 0.009 s. Figure 8a,b shows the excellent match of the 3-D CFD pressure history at the valve and midpoint of the pipe with its experimental counterpart.

Figure 8 .
Figure 8.Comparison of pressure head (h) variation with time ( t) between CFD simulations and experimental data [34]: (a) at valve (b) at the midpoint of the copper pipeline.Water 2023, 15, x FOR PEER REVIEW 15 of 25

Figure 9 .
Figure 9. Part of the pipe geometry shows the gate valve moving inside the pipe [65].

Figure 9 .
Figure 9. Part of the pipe geometry shows the gate valve moving inside the pipe [65].

Figure 11 .
Figure 11.(a) fluid domain and immersed solid domain at  1 , (b) overlapping of the fluid domain and immersed solid domain at the time  2 for modelling valve closure [65].

Figure 11 .
Figure 11.(a) fluid domain and immersed solid domain at t 1 , (b) overlapping of the fluid domain and immersed solid domain at the time t 2 for modelling valve closure [65].

Figure 13 .
Figure 13.Typical velocity profile in a decelerating pipe flow ( 1 >  2 >  3 ).Most fluid conveyance methods through pipes are not steady from the industrial perspective.In addition, control valves or pumps continuously adjust the flow through pipes to satisfy certain demands.Most engineering designs exclude transient effects in unsteady flows because of a lack of understanding about such flows.Hence, to estimate the frictional drag in accelerating, as well as decelerating, pipe flow, not many valid unstable friction models are available.Uncertainty regarding the transitional mechanisms and durations of the various transitional stages may result from the chosen step-down size or the masking impact caused by the flow structures at the start of the transient.The study of turbulence structures in accelerated and decelerated flows is important to an in-

Figure 13 .
Figure 13.Typical velocity profile in a decelerating pipe flow ( 1 >  2 >  3 ).Most fluid conveyance methods through pipes are not steady from the industrial perspective.In addition, control valves or pumps continuously adjust the flow through pipes to satisfy certain demands.Most engineering designs exclude transient effects in unsteady flows because of a lack of understanding about such flows.Hence, to estimate the frictional drag in accelerating, as well as decelerating, pipe flow, not many valid unstable friction models are available.Uncertainty regarding the transitional mechanisms and durations of the various transitional stages may result from the chosen step-down size or the masking impact caused by the flow structures at the start of the transient.The study of turbulence structures in accelerated and decelerated flows is important to an in-
k is the turbulence kinetic energy generation due to the mean strain rate; Y k and Y ω are the dissipation of k and ω due to turbulence; G ω is the generation of ω; Γ k and Γ ω are the effective diffusivity of k and ω, respectively; D ω represents the cross-diffusion term; and S k and S ω are source terms.Saemi et al. [55] simulated 2-D axisymmetric calculations for laminar as well as turbulent water hammer flows using unsteady Reynoldsaveraged Navier-Stokes equations.They employed the low Reynolds SST k-ω model and high Reynolds RNG k-ε model and found that the low Reynolds SST k-ω model can produce better pressure variations that match with experimental results.Kalantar et al. [61] and Saemi et al. [62] applied the SST k-ω model for the 3-D simulation of water hammer problems.