Next Article in Journal
Effects of Drip Irrigation Design on a Lemon and a Young Persimmon Orchard in Semi-Arid Conditions
Previous Article in Journal
Retrieving and Verifying Three-Dimensional Surface Motion Displacement of Mountain Glacier from Sentinel-1 Imagery Using Optimized Method
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Water Hammer Simulation Method in Pressurized Pipeline with a Moving Isolation Device

1
College of Naval Architecture and Ocean Engineering, Dalian Maritime University, Dalian 116026, China
2
School of Civil, Environmental and Mining Engineering, University of Adelaide, Adelaide, SA 5005, Australia
3
College of Mechanical and Transportation Engineering, China University of Petroleum-Beijing, Beijing 102249, China
4
Power China Kunming Engineering Corporation Limited, Kunming 650051, China
*
Author to whom correspondence should be addressed.
Water 2021, 13(13), 1794; https://doi.org/10.3390/w13131794
Submission received: 28 May 2021 / Revised: 23 June 2021 / Accepted: 28 June 2021 / Published: 29 June 2021

Abstract

:
Smart isolation devices (SIDs) are commonly used in pressurized subsea pipelines that need to be maintained or repaired. The sudden stoppage of the SID may cause large water hammer pressures, which may threaten both the pipeline and the SID. This paper proposes a simulation method by using a coupled dynamic mesh technique to simulate water hammer pressures in the pipeline. Unlike other water hammer simulations, this method is the first to be used in the simulation in pipelines with a moving object. The implicit method is applied to model the moving SID since it has the mutual independence between the space step and the time step. The movement of the SID is achieved by updating the size of the computational meshes close to the SID at each time step. To improve the efficiency of the simulation and the ability of handling complex boundary conditions, the pipe sections far away from the SID can also be simulated by using the explicit Method of Characteristics (MOC). Verifications were conducted using the simulated results from the Computational Fluid Dynamics (CFD) numerical simulation. Two scenarios have been studied and the comparisons between the simulated results by using the dynamic meshes in 1D methods and those by the CFD simulation show a high correlation, thus validating the new method proposed in this paper.

1. Introduction

Corrosions and cracks are common in aged water, oil and gas subsea pipelines. Strategically targeted pipeline maintenance, replacement and rehabilitation are of critical importance in these pipelines. In recent years, subsea pipeline maintenance technologies using smart isolation devices (SIDs) have been widely used [1,2,3]. The SIDs are usually used in pairs. The first SID driven by the flowing water moves in the pipeline with a similar velocity as for the flowing water, as shown in Figure 1a. When it approaches the targeted section and is allowed to stop, it starts to decelerate suddenly by extending some external slips to press against the pipeline internal wall (Figure 1b). After a short deceleration process (normally less than 1 s), the SID will stop at the desired location. By opening the valve (SID valve) in the middle of the SID, the SID then becomes hollow to keep the water flowing in the pipe (Figure 1c). Another SID can be then injected into the pipeline with the same procedure to be stopped at a location upstream of the original SID (Figure 1d). Finally, the valve of the original SID is closed and the packers around two SIDs will expand to seal the gap between the SID and the pipeline internal wall to fully stop the water flow (Figure 1e). Once the pair of SIDs are in place, the pipeline section between two SIDs is then isolated and can be repaired.
During the isolation process, the sudden deceleration of the SID may cause large pressure surges (water hammer pressures) at the front and back surfaces of the SID in the pipeline. This may lead to potential damage to the SID and the pipeline and thus the failure of the isolation. In addition, the induced water hammer pressures and the friction between the slips of the SID and the pipeline wall will determine the deceleration process and thus the final stopping position. An external communication link (ECL) needs to be placed close enough to the SID to ensure successful signal transmission with the SID. Therefore, the SID stopping position needs to be accurately predicted. This paper describes the development of a simulation method to determine the water hammer pressures in the isolation process to improve the operation safety of the SID and the reliability of data communication.
The simulation of water hammer pressures in pipelines has been studied by many researchers [4]. As one kind of explicit method, the Method of Characteristics (MOC) [5,6] is the most popular method to simulate water hammer pressures because of its simplicity in coding, high efficiency of simulation and its capacity to incorporate a variety of boundary conditions. There are also many other improved explicit methods that have been proposed for water hammer simulation. Chaudhry and Hussaini [7] solved the water hammer equations by using the MacCormack, Lambda and Gabutti explicit Finite Difference (FD) scheme. Compared with simulation results from the first-order MOC method, it is shown that for the same accuracy, second-order schemes require fewer computational nodes and less computer time as compared to those required by the first-order characteristic method. In addition, the first-order and second-order explicit Finite Volume (FV) methods of Godunov-type for solving water hammer equations were formulated by Hwang and Chung [8], Zhao and Ghidaoui [9] and Ghidaoui et al. [4] The simulation results show that the first-order FV Godunov-type scheme can achieve very similar outcomes to that from MOC, and the second-order Godunov-type scheme can achieve the same accuracy with MOC but with less memory storage and execution time.
The explicit methods mentioned above are efficient and robust when solving hydraulic transient problems in long pipelines or large pipe networks with complex boundary conditions. However, the space step and time step of these methods are dependent and have to meet the Courant-Friedrichs-Lewy stability condition (CFL = a · Δ t / Δ x ) to achieve the stability of the simulation [5]. Thus, the dynamic meshes with adjustable space steps, which can be used to simulate the moving SID, are difficult to achieve in these explicit methods.
A spatial four-point finite-difference implicit method was presented by Jin et al. [10] to simulate the unsteady flow process of a drainage system with accurate results achieved even using a large space step. Another implicit method of characteristics (IMOC) was proposed by Afshar and Rohani [11] with the simulation validated by the results from the MOC. Since the simulation using the implicit method can achieve unconditional convergence and mutual independence between the space step and the time step, the space step can be flexibly changed without changing the time step during the simulation. Thus, the implicit method is predicted to have the potential to be used with dynamic meshes (where the mesh size changes with time step), which is essential to simulate a moving object in pipes.
To further improve the efficiency and robustness of the hydraulic transient simulation, Wang and Yang [12] combined the explicit MOC with the Preissmann four-point finite-difference implicit method [13] to simulate the unsteady flow in pipelines and the transient processes in hydropower systems. The results show that the coupling method is effective. This method, however, has only been applied to pipe systems without a moving object.
In this paper, a dynamic mesh technique in the implicit method has been developed to simulate the hydraulic transient process with a moving SID in a pipeline. The pipeline was divided into five sections. The lengths of two sections change during the simulation, and so do the sizes of the meshes in these sections. A linear interpolation method has been used to obtain values at the updated mesh nodes. The explicit MOC has been also used to be coupled with the implicit method (referred to reference [12]) to improve the computing efficiency and the capability of handling complex boundary conditions. To validate this novel simulation approach, numerical verifications using Computational Fluid Dynamics (CFD) and numerical simulation have been conducted. The comparison between the CFD simulation and the simulation using the new method illustrates the effectiveness of the proposed approach for the water hammer simulation of pipelines with a moving object, such as an SID.

2. Methods to Model the SID Movement in Pipelines

2.1. The Method to Simplify the Moving SID in the Pipeline

As shown in Figure 2, the moving SID in the pipeline can be regarded as a moving cylinder. In order to achieve 1D simulation, the SID inside the pipeline has been simplified to a diameter-reduced section, and the movement of the SID is considered to be equivalent to the change of the lengths of the pipe Section S2 and S4 as shown in the schematic in Figure 3.
The area (A3) of the simplified section in Figure 3 is assumed to be equal to the clearance area (A1A2) between the SID and the pipe wall when the SID is moving in Figure 2, as expressed by
A 3 = A 1 A 2

2.2. The Division of the Pipeline with Corresponding Methods

In order to simulate the water hammer pressures in pipelines with a moving diameter-reduced section, as mentioned above, the whole pipeline has been divided into five sections, as shown in Figure 3. Section S2, S3 and S4 are solved using the implicit method, and the space steps are independent of the time step, making the space steps flexible to change during the simulation. The explicit MOC method has been used in Section S1 and S5 due to its high computational efficiency and easy boundary condition setting at the inlet and outlet nodes. It should be noted that the implicit method was also applied in Section S1 and S5 in the following case studies for comparison. Overall, the major challenges to achieve the simulation are: (a) the implicit scheme with dynamic meshes in Section S2, S3 and S4 with the moving boundaries at the diameter-reduced interfaces (Interface 2 and 3 in Figure 3) and (b) the explicit–implicit coupling method at Interface 1 and 4 in Figure 3.

3. Discretization of the Governing Equations

3.1. Explicit Method of Characteristics

The 1D continuity and momentum equations governing the unsteady flow in pipelines, neglecting the convective acceleration terms, can be described as follows [6]:
H t + a 2 g A Q x = 0
H x + 1 g A Q t + f Q | Q | 2 g D A 2 = 0
in which H = piezometric head; Q = discharge; x = distance; t = time; a = wave speed; g = gravitational acceleration; f = Darcy–Weisbach friction factor; D = pipe diameter; and A = cross-sectional area.
The explicit MOC can be used for solving the governing equations. The partial differential equations (Equations (2) and (3)) can be transformed to ordinary differential equations by
C + :   ( Q i j + 1 Q i 1 j ) + g A a ( H i j + 1 H i 1 j ) + f 2 D A Q i 1 j | Q i 1 j | Δ t = 0
C :   ( Q i j + 1 Q i + 1 j ) g A a ( H i j + 1 H i + 1 j ) + f 2 D A Q i + 1 j | Q i + 1 j | Δ t = 0
which are valid along the positive ( C + ) and negative ( C ) characteristic lines, as shown in Figure 4. The subscripts i − 1, i and i + 1 represent the space steps, and the superscripts j and j + 1 represent the time steps. By rearranging Equations (4) and (5), the following compatibility equations can be obtained
C + :   Q i j + 1 = C P B · H i j + 1
C :   Q i j + 1 = C M + B · H i j + 1
in which C P = Q i 1 j + B · H i 1 j R · Q i 1 j | Q i 1 j | Δ t ; C M = Q i + 1 j B · H i + 1 j R · Q i + 1 j | Q i + 1 j | Δ t ; B = g A / a ; R = f / 2 D A . According to the space-time plane shown in Figure 4, Q i j + 1 and H i j + 1 at time step j + 1 can be calculated by using Equations (6) and (7) if Q i 1 j , Q i + 1 j , H i 1 j and H i + 1 j at time step j have already been obtained from the calculation in the time step j-1 or are already known (e.g., steady state conditions).

3.2. Discrete Implicit Method

The Preissmann four-point finite-difference implicit method (referred to as the “implicit method” in the following) is introduced in this section. The discretization scheme of this implicit method with a temporal weighting coefficient θ is represented by [13]
y x = θ y i + 1 j + 1 y i j + 1 Δ x + ( 1 θ ) y i + 1 j y i j Δ x
y t = y i + 1 j + 1 y i + 1 j + y i j + 1 y i j 2 Δ t
in which y can stand for the piezometric head or the flow rate. Replacing the partial derivatives of Equations (2) and (3) by the above finite-difference scheme yields
H i + 1 j + 1 H i + 1 j + H i j + 1 H i j 2 Δ t + a 2 g A [ θ Q i + 1 j + 1 Q i j + 1 Δ x + ( 1 θ ) Q i + 1 j Q i j Δ x ] = 0
θ H i + 1 j + 1 H i j + 1 Δ x + ( 1 θ ) H i + 1 j H i j Δ x + 1 g A Q i + 1 j + 1 Q i + 1 j + Q i j + 1 Q i j 2 Δ t + f 2 g D A 2 ( Q i + 1 j | Q i + 1 j | + Q i j | Q i j | 2 ) = 0
A further rearrangement of Equations (10) and (11) gives
I 1 · Δ H i + 1 + J 1 · Δ Q i + 1 = L 1 · Δ H i + M 1 · Δ Q i + N 1
I 2 · Δ H i + 1 + J 2 · Δ Q i + 1 = L 2 · Δ H i + M 2 · Δ Q i + N 2
where I 1 = 1 ; J 1 = 2 a 2 · θ · Δ t g · A · Δ x ; L 1 = 1 ; M 1 = 2 a 2 · θ · Δ t g · A · Δ x ; N 1 = 2 a 2 · Δ t g · Δ x ( Q i + 1 j Q i j A ) ; I 2 = θ · Δ t Δ x ; J 2 = 1 2 g · A ; L 2 = θ · Δ t Δ x ; M 2 = 1 2 g · A ; N 2 = Δ t Δ x ( H i + 1 j H i j ) f · Δ t 4 g ( Q i + 1 j | Q i + 1 j | D · A 2 Q i j | Q i j | D · A 2 ) . All the coefficients are related to the heads and discharges at the time step j. Δ H and Δ Q are defined by
Δ H i = H i j + 1 H i j
Δ Q i = Q i j + 1 Q i j
and denote the increments of piezometric head and discharge over the parameters at the previous time step. By combining Equations (12) and (13) with the boundary conditions, the relationship of the discharge increment and piezometric head increment (referred to as the increment equation) at the first node can be defined as
Δ Q 1 = E E 1 · Δ H 1 + F F 1
where E E 1 and F F 1 are two terms that are related to the boundary conditions, the piezometric head and the discharge at the first node at the time step j. The increment equation at the second node at time step j + 1 can be then obtained by incorporating Equations (12) and (13) into (16). By repeating the processes, the increment equations at all the nodes at time step j+1 can be obtained as
Δ Q i = E E i · Δ H i + F F i
where E E i and F F i are two terms that are related to the boundary conditions, the piezometric head and the discharge at the node i at the time step j. This process is known as the forward scan of the chasing method [12]. By combining the increment equation at the end node with the boundary conditions at this node, the piezometric head and discharge at the end node can be obtained. Similarly, a backward scan of the chasing method was used by combining the increment equation at node i with the solved piezometric head and discharge at node i+1 to solve the piezometric heads and discharges at all the nodes.

4. Coupling Method with Dynamic Meshes

4.1. Dynamic Meshes in the Implicit Method

The diameter-reduced Section S3 representing the SID as shown in Figure 5 moves at the velocity of V which is equal to the initial flow velocity. The size of the meshes in Section S3 is fixed since the length of Section S3 remains unchanged during its movement. The movement of Section S3, however, changes the lengths of Section S2 and S4 in Figure 5, increasing for Section S2 and decreasing for Section S4. Lengths of S2 and S4 should be longer than the moving distance of the SID during the deceleration process. By setting fixed mesh numbers for Section S2 and S4, the size of the meshes in these sections should be able to be altered during the simulation to enable the length change of Section S2 and S4.
The number of nodes of Section S2 and S4 are assumed to be n 1 and n 2 , respectively. The lengths of these two sections are L3 and L4 as shown in Figure 5. Thus, the mesh sizes for Section S2 and S4 are Δ x 1 = L 3 / ( n 1 1 ) and Δ x 2 = L 4 / ( n 2 1 ) with uniform meshes used in these sections. During one time step Δ T , the distance that S3 moves can be calculated as Δ L = ( V j + V j + 1 ) Δ T / 2 with V representing the velocity of the SID, and the lengths of Section S2 and S4 change to L 3 + Δ L and L 4 Δ L , respectively. As a result, the mesh sizes of Section S2 and S4 need to be updated to Δ x 1 = ( L 3 + Δ L ) / ( n 1 1 ) and Δ x 2 = ( L 4 Δ L ) / ( n 2 1 ) . A simple case with only 4 nodes is presented in Figure 6.
Since the distance moved of Δ L is much smaller than the initial mesh sizes at each time step, the piezometric head and discharge between two adjacent nodes (such as node i and node i + 1 in Figure 6) can be assumed to be varying linearly. Thus, the values at the updated nodes in Section S2 and S4 can be obtained by a linear interpolation method, as shown in Figure 6. The piezometric head and discharge at the updated node i can be expressed by
H i j = ( H i + 1 j H i j ) x i x i x i + 1 x i + H i j
Q i j = ( Q i + 1 j Q i j ) x i x i x i + 1 x i + Q i j
The piezometric heads and discharges at the next time step j + 1 can be then calculated using the parameters at the updated mesh nodes in the implicit method.

4.2. Implicit Method in Pipelines with a Moving Diameter-Reduced Section

In this section, the implicit method has been applied to the whole pipeline and the boundary conditions at two ends of the pipe have been set according to the real field situation. Examples include constant heads at both upstream and downstream ends of the pipe. These values have been used in the case studies in the paper.
According to the forward scan of chasing method with the upstream boundary conditions, the increment equation at the last node of Section S2 can be obtained as
Δ Q S 2 , n 1 = E E S 2 , n 1 · Δ H S 2 , n 1 + F F S 2 , n 1
where the subscripts S1, S2, S3, S4 and S5 refer to the pipeline sections as defined in Figure 3, and the numbers or characters following these subscripts are the node number. For example, the subscript (S2, n1) means the n1th node of Section S2.
At the connection node of Section S2 and S3 (Interface 2 in Figure 5), the energy equation is
H S 2 , n 1 j + 1 + Q S 2 , n 1 j + 1 · Q S 2 , n 1 j + 1 2 g A S 2 , n 1 2 = H S 3 , 1 j + 1 + Q S 3 , 1 j + 1 · Q S 3 , 1 j + 1 2 g A S 3 , 1 2 + ξ 1 | Q S 3 , 1 j + 1 | · Q S 3 , 1 j + 1 2 g A S 3 , 1 2
in which
ξ 1 = ( 1 C c 1 ) 2
Equation (22) is the minor loss coefficient at Interface 2 in Figure 5 for a sudden pipeline contraction. C c is the contraction coefficient [5], which can be determined according to the area ratio A S 3 / A S 2 .
Since the diameter-reduced section is moving with the fluid, the continuity equation at the connection node needs to include both the flowing fluid and the moving section. Thus, it can be expressed as
Q S 2 , n 1 j + 1 = Q S 3 , 1 j + 1 + V j + 1 · A 2
Therefore, the increment equation at the first node in Section S3 can be then obtained by combining Equations (20), (21) and (23), which gives
Δ Q S 3 , 1 = E E S 3 , 1 · Δ H S 3 , 1 + F F S 3 , 1
where E E S 3 , 1 = 1 Y ; F F S 3 , 1 = Z Y · Q S 3 , 1 j Y ; Y = Q S 3 , 1 j 2 g A S 2 , n 1 2 ( 1 + ξ 1 ) Q S 3 , 1 j 2 g A S 3 , 1 2 + 1 E E S 2 , n 1 + V j + 1 · A 2 g A S 2 , n 1 2 ;
Z = H S 3 , 1 j H S 2 , n 1 j + Q S 2 , n 1 j + F F S 2 , n 1 V j + 1 · A 2 E E S 2 , n 1 V j + 1 · V j + 1 · A 2 2 2 g A S 2 , n 1 2 ;
Along with Equation (24), the increment equations at other remaining nodes of Section S3 can be obtained by the forward scan of chasing method.
Similarly, the increment equation at the last node of Section S3, the continuity equation and energy equation and at Interface 3 in Figure 5 can be expressed by
Δ Q S 3 , n = E E S 3 , n · Δ H S 3 , n + F F S 3 , n
Q S 3 , n j + 1 = Q S 4 , 1 j + 1 V j + 1 · A 2
H S 3 , n j + 1 + Q S 3 , n j + 1 · Q S 3 , n j + 1 2 g A S 3 , n 2 = H S 4 , 1 j + 1 + Q S 4 , 1 j + 1 · Q S 4 , 1 j + 1 2 g A S 4 , 1 2 + ξ 2 | Q S 3 , n j + 1 | · Q S 3 , n j + 1 2 g A S 3 , n 2
where ξ 2 is the minor loss coefficient at Interface 3 for a sudden pipeline expansion and it can be calculated through
ξ 2 = ( 1 A S 3 , n A S 4 , 1 ) 2
The increment equation at the first node in Section S4 can be then obtained by substituting Equations (25) and (26) into (27), which gives
Δ Q S 4 , 1 = E E S 4 , 1 · Δ H S 4 , 1 + F F S 4 , 1
where E E S 4 , 1 = 1 Y 1 ; F F S 4 , 1 = Z 1 Y 1 · Q S 4 , 1 j Y 1 ; Y 1 = Q S 4 , 1 j 2 g A S 3 , n 2 ( 1 + ξ 2 ) Q S 4 , 1 j 2 g A S 4 , 1 2 + 1 E E S 3 , n V j + 1 · A 2 g A S 3 , n 2 ;
Z 1 = H S 4 , 1 j H S 3 , n j + Q S 3 , n j + F F S 3 , n + V j + 1 · A 2 E E S 3 , n V j + 1 · V j + 1 · A 2 2 2 g A S 3 , n 2 ;
Along with Equation (29), the increment equations at all nodes in Section S4 and S5 can be obtained by the forward scan of chasing method. Once all the increment equations are known, the backward scan of the chasing method can be applied, starting from the end node of Section S5 (combined with the downstream boundary conditions) to the first node of Section S1, to solve the discharges and piezometric heads at all nodes.
If Section S1 and S5 in Figure 3 are simulated using the MOC, the upstream and downstream boundary conditions of the implicit sections used in the forward and backward scan of chasing methods, respectively, are no longer the boundary conditions at the pipe ends. They change to the explicit–implicit coupling boundaries at Interface 1 and 4 (shown in Figure 3) given in the following section and the range of the scan is reduced to Section S2, S3 and S4.

4.3. Explicit–Implicit Coupling Method

In practice, there may be some complex boundary conditions at the inlet and outlet of the pipe, which may be difficult to formulate using the implicit method. These boundary conditions can be easily handled in the explicit MOC, which is numerically stable and computationally efficient. Based on the pipeline properties and the surrounding environment in the field, the unsteady friction loss in the pipe [14], the viscoelasticity of the pipe wall [15,16], the fluid-structure interaction [17] and water-column separation [18] may need to be considered, and the strategies to incorporate them into the MOC have been extensively investigated and can be utilized if needed. Thus, the transient simulation of pipelines with a moving object inside can be extended to more complicated pipeline configurations if it can be coupled with the MOC. To achieve that, the explicit–implicit coupling boundaries have been applied and illustrated in the following way.
The explicit–implicit coupling boundaries are shown in Figure 7. The C + compatibility equation for the end node of Section S1, the continuity equation and the energy equation at Interface 1 as shown in Figure 7 can be written as follows:
Q S 1 , n 3 j + 1 = C P B · H S 1 , n 3 j + 1
Q S 1 , n 3 j + 1 = Q S 2 , 1 j + 1
H S 1 , n 3 j + 1 + Q S 1 , n 3 j + 1 · Q S 1 , n 3 j + 1 2 g A S 1 , n 3 2 = H S 2 , 1 j + 1 + Q S 2 , 1 j + 1 · Q S 2 , 1 j + 1 2 g A S 2 , 1 2 + ξ 3 | Q S 2 , 1 j + 1 | · Q S 2 , 1 j + 1 2 g A S 2 , 1 2
where ξ 3 is the minor loss coefficient for a sudden contraction between Section S1 and S2. In the case shown in Figure 7, ξ 3 = 0 since the diameters in these two sections are the same.
By combining Equations (30)–(32), the increment equation at the first node of Section S2 can be expressed as
Δ Q S 2 , 1 = E E S 2 , 1 · Δ H S 2 , 1 + F F S 2 , 1
where E E S 2 , 1 = B , F F S 2 , 1 = C P B · H S 2 , 1 j Q S 2 , 1 j .
Along with Equation (33), the increment equation at the last node of Section S2 can be obtained according to the forward scan of chasing method. Then, by conducting the same procedures with those presented in Section “Implicit method in pipelines with a moving diameter-reduced section”, the increment equations at all the nodes in the implicit sections can be obtained, including the one at the end node of Section S4, as shown in Equation (34).
Δ Q S 4 , n 2 = E E S 4 , n 2 · Δ H S 4 , n 2 + F F S 4 , n 2
At Interface 4, Equation (34) can be combined with the continuity equation and energy equation at the interface and the C compatibility equation at the first node of Section S5. The piezometric head and discharge at the end node of Section S4 can be then obtained. A similar backward scan of chasing method can be applied to solve the piezometric heads and discharges at all the nodes of the implicit sections.
The calculation flow chart of the simulation using the dynamic mesh technique in the coupling method with a moving diameter-reduced section is shown in Figure 8.

5. CFD Validation of the New Moving Mesh Method

5.1. System Configuration

A horizontal water pipeline with a moving SID shown in Figure 2 is considered in this section to verify the proposed method. The inner diameter of the pipeline is D 1 = 0.1 m, the diameter of the SID is D 2 = 0.08 m. Thus the equivalent diameter of the simplified section in the 1D model in Figure 3 is D 3 = 0.06 m. Other properties and pipe parameters are presented in Table 1. The friction losses along the sections (S1, S2, S4 and S5) in both models are the same by choosing an equivalent roughness height in the CFD model corresponding to the Darcy–Weisbach friction factor f in the 1D model. The boundary conditions at the inlet and outlet of the pipe for all the simulations are set to be constant pressure heads according to Table 1.

5.2. Basic Information of the CFD Model

The CFD simulation was conducted on the ANSYS Fluent with an unsteady and two-dimensional axisymmetric model. The geometry of the numerical model is shown in Figure 9. Before the simulation of the SID retardation process, the process that the SID flows with the water at 3 m/s was simulated for a period of 50 s, which is enough to initialize the flow field. To make the SID begin to stop at the same position as that in the 1D simulation model, the initial position of the SID is 150 m (3 m/s × 50 s) forward and thus L 1 = 150 m and L 2 = 450 m in the CFD model, as shown in Figure 9.
Two different simulation zones are included in this model. The first one is a fluid zone which contains the fluid (water) in the pipeline. The other one is a moving solid zone, which is used to represent the moving SID. The inlet and outlet boundary conditions were each set with a constant pressure value. Structured square meshes have been used in this model with the size Δ l of 0.001 m. The time step, which can be estimated by Δ t Δ l / a , was set to be 10−6 s. The RNG k- ε turbulence model has been used in the simulation, which has been tested to be the most suitable model for turbulent flow simulation [19,20,21]. The coupling of pressure and velocity was achieved by using the first-order implicit method for the pressure-linked equation (SIMPLE) algorithm. Second-order upwind discretization was used for the momentum equation.
The motion of the SID was implemented by means of the user-defined function (UDF) with dynamic meshes (it should be noted that the dynamic meshes mentioned here are in the CFD model, which are different to those of the dynamic meshes in the 1D simulation model). Due to the superior ability to deal with structured meshes under linear movement in the CFD model, the layering method combined with the spring-based smoothing method was applied to update the structured meshes in the deforming regions (Region R2 and R4 in Figure 9) subjected to the wall motion. In essence, with the movement of SID, the front and back deforming regions advanced a prescribed threshold, and therefore, the complete rows of meshes would be created or collapsed [22,23].
The area-weighted average pressure head at the back surface of the SID (as shown in Figure 9) will be compared with the results from the implicit method as well as the coupling method to validate the proposed methods. In the paper, the following two scenarios have been simulated: Scenario 1: The SID moves in the pipeline and stops immediately (in one time step); Scenario 2: The SID moves in the pipeline and stops with a deceleration process.
The independence analysis of time step and space step size on the CFD simulations have been conducted to eliminate the influence of the time step and space step size on the results. In this simulation, the initial velocity of the SID was 3 m/s and it was set to decelerate with a constant acceleration, and it stopped completely after a one-second deceleration process. The pressure as the back surface of the SID has been monitored to analyze the influence of the time step and space step size. Figure 10 and Figure 11 show the different values employed in the study. It can be concluded that the time step and space step size used in this paper are found sufficiently independent to give excellent results in terms of the numerical precision, with deviations lower than 0.5%.

5.3. Steady State Flow

Before the transient simulation, the steady state flow for the 1D model is calculated and the friction losses are calibrated to make the steady state same in the CFD simulation.
The friction losses in Section S1, S2, S4 and S5 (as shown in Figure 3) should be the same in both 1D model and the CFD model, and can be calculated by
h f 1 = f · ( L 1 + L 2 ) D 1 · Q S 1 2 2 g A 1 2
For Section S3, the friction loss in the 1D model is calculated based on the non-simplified model as shown in Figure 2. Thus, it is
h f 2 = f · L s D e · Q S 3 2 2 g A 3 2
where D e = 4 R with R = A 3 / P representing the hydraulic radius and P is the wetted perimeter for pipes with a non-circular cross-section.
Another factor that affects the steady state is the minor loss when the water flows across the SID or the diameter-reduced section. The total minor loss in the 1D model for both the sudden contraction and sudden expansion can be calculated by
h L = ξ 1 | Q S 3 | · Q S 3 2 g A 3 2 + ξ 2 | Q S 4 | · Q S 4 2 g A 1 2
in which ξ 1 and ξ 2 can be estimated using Equations (22) and (28) with the SID simplified to the diameter-reduced section. Due to the simplification, the calculated minor loss may be different to that in the original model (non-simplified model) used in CFD simulation. Thus, a coefficient φ has been applied to the Equation (37) in the 1D model as
h L 1 = φ ( ξ 1 | Q S 3 | · Q S 3 2 g A 3 2 + ξ 2 | Q S 4 | · Q S 4 2 g A 1 2 )
which is to compensate for the error of the minor loss calculation caused by the simplification.
The steady state of the CFD model was simulated using the CFD tool with the SID immobile at the initial position. The total head loss H f in the SID section in the CFD simulation is equal to the pressure head difference between the back surface and the front surface of the SID (as shown in Figure 9), and the result is 4.2 m. By assuming the simulated friction loss along the SID section (excluding the minor losses) is equal to the value calculated by Equation (36), the total minor loss in the SID section in the CFD simulation is
h L = H f h f 2
To make the minor loss the same in both 1D and CFD simulations, the coefficient φ in Equation (38) can be calculated by
φ = H f h f 2 h L
which is equal to 1.1 in this case. With the minor loss calibrated in the 1D model, the steady state will be the same as for the CFD model.

5.4. Scenario 1: The SID Moves in the Pipeline and Stops Immediately

In this scenario, the whole pipeline has been simulated by the implicit method with the diameter-reduced section stopping immediately (in one time step). The time step and the space step, which are used in all the 1D simulations in the paper, are set to be the same as that of the CFD simulation. For the CFD simulation (as shown in Figure 9), the time-varying velocity of the SID is given by
V ( t ) = {                   3   m / s         f o r   0 < t 50   s 0                           f o r   t > 50   s
The comparison between the piezometric head at the back surface of the SID simulated by CFD and that simulated by the implicit method is shown in Figure 12.
It can be observed from Figure 12 that the results from these two kinds of simulations are in good agreement, and that the maximum deviation between the two curves occurs at the peak point A. The maximum relative error between the 1D simulation and the CFD simulation is less than 2%, which illustrates that the proposed method is valid for simulating the water hammer pressures in pipelines with a moving object inside. At the first time step after the SID stops, there is an obvious pressure spike, as indicated by point A in Figure 12. This can be ascribed to the superposition of two pressure waves induced at the back and front surfaces of the SID. Due to the sudden velocity decrease of the water in front of the back surface, a positive pressure will be generated at this position. In the front surface of the SID, the sudden velocity decrease of the water will generate a negative pressure wave. The negative pressure wave can propagate backwards to the back surface after t = L S / a (0.0006 s) and become superimposed with the positive pressure wave. Thus, a sharp spike can be generated by the superposition of the pressure waves, and the duration of the spike is expected to be 0.0006 s, which is consistent with the duration observed in Figure 12.
To further validate the simulation, the Joukowsky Equation [24]
Δ H = a · Δ V g
has been used to calculate the theoretical head rise (defined as the head increase from the initial piezometric head to the maximum piezometric head), which can be then compared with the simulated head rise.
The Joukowsky head rise is calculated to be 305.8   m using Equation (41) if the velocity decreases to 0 in the pipe. The value is quite difference in comparison to the simulated head rise of 77 m (220.7 m–297.7 m) as shown in Figure 12. The difference is due to the clearance area in which water keeps flowing after the SID stops.
By decreasing the clearance area (from the right side to the left side of the abscissa in Figure 13), the simulated head rise approaches the theoretical head rise. When the clearance area gets close to zero, the simulated head rise is then equal to the Joukowsky value calculated from Equation (41). This comparison further validates the proposed 1D implicit method.

5.5. Scenario 2: The SID Stops with a Deceleration Process

In this section, the process that the SID moves and stops with a deceleration process has been simulated using the proposed dynamic mesh technique with 1D simulation methods. Two 1D methods were used in the simulation, including the purely implicit method and the explicit–implicit coupling method (as shown in Figure 3). The moving SID was set to decelerate with a constant acceleration, and it stopped completely after a one-second deceleration process. The time-varying velocity of the SID in the pipeline is given by
V ( t ) = {                       3   m / s                       f o r   0 < t 50   s 3 3 × ( t 50 )           f o r   50   s < t 51   s     0                                       f o r   t > 51   s

5.5.1. Validation of the Dynamic Mesh Technique in the Implicit Method

The whole pipeline was simulated using the implicit method in this sub-section. The simulated result, shown as the dashed-dotted line in Figure 14, shows a reasonable agreement with the CFD simulation, and the maximum relative error between the two curves is acceptably small (less than 2%). The errors are found periodically—over prediction at the peak values and under prediction just after the peak points. They can be ascribed to the following reasons: (a) Although a coefficient φ has been applied to make the SID equivalent to a diameter-reduced section in the steady state, the simplification may still affect the transient process; (b) local turbulences are not considered in the 1D simulation but are captured in the CFD simulation; (c) the pressures at Interface 2 and 3 are assumed to be uniform across the section in the 1D simulation but should be non-uniform in the CFD simulation; and (d) the piezometric head and flowrate between two adjacent mesh nodes in the 1D simulation are assumed to change linearly in the mesh updating process.

5.5.2. Validation of the Explicit–Implicit Coupling Method with Dynamic Meshes

The explicit–implicit coupling method has been verified in this section with Section S1 and S5 (LE1 = LE2 = 270 m, as shown in Figure 3, the lengths of S2 and S4 should be longer than the moving distance of the SID during the deceleration process) modeled by the MOC method and other sections modeled by the implicit method. In this model, the number of nodes for Section S1, S2, S3, S4 and S5 (as shown in Figure 5 and Figure 7) are n = ( 0.6 / Δ l ) + 1 = 601 , n 1 / 2 = ( 30 / Δ l ) + 1 = 30 , 001 , n 3 / 4 = ( 270 / Δ l ) + 1 = 270 , 001 , respectively. The simulated piezometric head at Interface 2 has been plotted as the dashed line in Figure 14. It can be seen that the simulated result using the coupling method is almost identical to that by the implicit method and is in good agreement with that in the CFD simulation. Thus, the proposed coupling method with dynamic meshes is valid to simulate the transient pressures in pipeline systems with a moving object inside.

6. Conclusions

A novel dynamic mesh technique in the implicit method has been developed in this paper, which can achieve a simulation of the water hammer pressures in maintaining subsea pipelines with a moving object, such as an SID. The whole pipeline has been divided into five sub-sections with the SID simplified to a moving diameter-reduced section in the simulation. The movement of the SID has been achieved by updating the sizes of the meshes in sections conjoined with the SID. A linear interpolation has been applied to update the values at the updated mesh nodes. Two scenarios with the SID stopping immediately and with a deceleration process have been simulated using the proposed methods with validation using CFD numerical simulations.
The proposed implicit method with the novel dynamic mesh technique has been validated in both scenarios since the 1D simulation results are in good agreement with those from CFD simulation. Considering the errors caused by assuming uniform flow at every cross section in the 1D simulations, the proposed implicit method with the novel dynamic mesh technique has been proven to be accurate to simulate the water hammer pressures in pipelines with a moving SID.
The proposed implicit method with the novel dynamic mesh technique is also able to be coupled with the MOC, which is convenient to handle complicated boundary conditions and to incorporate factors affecting hydraulic transients, such as the unsteady friction and the viscoelasticity of the pipe wall. The coupling method has been validated by comparing the simulated result with those simulated by the purely implicit method and the CFD method. The coupling with the MOC makes the 1D simulation a potential candidate for simulating complicated pipe systems with a moving object inside.
The main limitation of this proposed method is that in order to achieve the dynamic mesh method, the moving distance Δ L of the SID at every time step should be always much smaller than the mesh size in the whole transient process. Then the linear interpolation method can be used in the proposed dynamic mesh method. In addition, the lengths of the Section S2 and S4 need to be determined in every case according to the deceleration distance of the SID, when using the explicit–implicit coupling method.

Author Contributions

Conceptualization, K.Z.; methodology, K.Z.; software, K.Z. and W.Z.; validation, K.Z., W.Z. and A.R.S.; investigation, K.Z.; writing—original draft preparation, K.Z.; writing—review and editing, A.R.S., S.Z. and C.W.; supervision, A.R.S. and S.Z.; funding acquisition, K.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China, grant number: U1908228; and the Fundamental Research Funds for the Central Universities, grant number: 3132021115.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Aleksandersen, J.; Tveit, E. The Smart Plug: A Remotely Controlled Pipeline Isolation System. In Proceedings of the Eleventh International Offshore and Polar Engineering Conference, Stavanger, Norway, 17–22 June 2001. [Google Scholar]
  2. Lie, M.R.G.; Muangsuankwan, M.N. Remote-Controlled Plugging Technology Minimizes Platform Downtime: Valve Replacement Through SmartPlug® Isolation. In Proceedings of the ASME 2015 India International Oil and Gas Pipeline Conference, New Delhi, India, 17–18 April 2015. [Google Scholar]
  3. Zhao, B.; Li, C.; Zhang, J.; Hu, Y. The Isolation Technology of Oil and Gas Pipeline in China. In Proceedings of the Twentieth International Offshore and Polar Engineering Conference, Beijing, China, 20–25 June 2010. [Google Scholar]
  4. Ghidaoui, M.S.; Zhao, M.; McInnis, D.A.; Axworthy, D.H. A Review of Water Hammer Theory and Practice. Appl. Mech. Rev. 2005, 58, 49–76. [Google Scholar] [CrossRef]
  5. Chaudhry, M.H. Applied Hydraulic Transients, 3th ed.; Springer: New York, NY, USA, 2014. [Google Scholar]
  6. Wylie, E.B.; Streeter, V.L.; Suo, L. Fluid Transients in Systems; Prentice Hall: Englewood Cliffs, NJ, USA, 1993. [Google Scholar]
  7. Chaudhry, M.H.; Hussaini, M.Y. Second-Order Accurate Explicit Finite-Difference Schemes for Waterhammer Analysis. J. Fluids Eng. 1985, 107, 523–529. [Google Scholar] [CrossRef]
  8. Hwang, Y.-H.; Chung, N.-M. A Fast Godunov Method for the Water-Hammer Problem. Int. J. Numer. Methods Fluids 2002, 40, 799–819. [Google Scholar] [CrossRef]
  9. Zhao, M.; Ghidaoui, M.S. Godunov-Type Solutions for Water Hammer Flows. J. Hydraul. Eng. 2004, 130, 341–348. [Google Scholar] [CrossRef]
  10. Jin, M.; Coran, S.; Cook, J. New One-Dimensional Implicit Numerical Dynamic Sewer and Storm Model. In Global Solutions for Urban Drainage; American Society of Civil Engineers: Reston, VA, USA, 2002; pp. 1–9. [Google Scholar] [CrossRef]
  11. Afshar, M.H.; Rohani, M. Water Hammer Simulation by Implicit Method of Characteristic. Int. J. Press. Vessel. Pip. 2008, 85, 851–859. [Google Scholar] [CrossRef]
  12. Wang, C.; Yang, J.-D. Water Hammer Simulation Using Explicit–Implicit Coupling Methods. J. Hydraul. Eng. 2015, 141, 04014086. [Google Scholar] [CrossRef]
  13. Preissmann, A. Propagation Des Intumescences Dans Les Canaux et Les Riviers. In Proceedings of the First Congress of the French Association for Computation, Grenoble, France, 15 June 1961. [Google Scholar]
  14. Bergant, A.; Ross Simpson, A.; Vìtkovsk, J. Developments in Unsteady Pipe Flow Friction Modelling. J. Hydraul. Res. 2001, 39, 249–257. [Google Scholar] [CrossRef] [Green Version]
  15. Covas, D.; Stoianov, I.; Mano, J.F.; Ramos, H.; Graham, N.; Maksimovic, C. The Dynamic Effect of Pipe-Wall Viscoelasticity in Hydraulic Transients. Part II—Model Development, Calibration and Verification. J. Hydraul. Res. 2005, 43, 56–70. [Google Scholar] [CrossRef]
  16. Covas, D.; Stoianov, I.; Ramos, H.; Graham, N.; Maksimovic, C. The Dynamic Effect of Pipe-Wall Viscoelasticity in Hydraulic Transients. Part I—Experimental Analysis and Creep Characterization. J. Hydraul. Res. 2004, 42, 517–532. [Google Scholar] [CrossRef]
  17. Tijsseling, A.A. FLUID-STRUCTURE INTERACTION IN LIQUID-FILLED PIPE SYSTEMS: A REVIEW. J. Fluids Struct. 1996, 10, 109–146. [Google Scholar] [CrossRef] [Green Version]
  18. Bergant, A.; Simpson, A.; Tijsseling, A. Water Hammer with Column Separation: A Historical Review. J. Fluids Struct. 2006, 22, 135–171. [Google Scholar] [CrossRef] [Green Version]
  19. Tang, H.; Zhong, S. 2D Numerical Study of Circular Synthetic Jets in Quiescent Flows. Aeronaut. J. 2005, 109, 89–97. [Google Scholar] [CrossRef] [Green Version]
  20. Zhang, X.; Wang, S.; You, W. 1D and CFD Co-Simulation Approach Basing on General Purpose Simulation Software. In Proceedings of the 2009 International Conference on Computational Intelligence and Software Engineering, Wuhan, China, 11–13 December 2009; Institute of Electrical and Electronics Engineers (IEEE): Piscataway, NJ, USA, 2009; pp. 1–4. [Google Scholar] [CrossRef]
  21. Zhang, X.; Wang, S.; Yu, D. 2D Axisymmetric CFD Simulation of Underwater Torpedo Launch Tube Flow. In Proceedings of the 2009 International Conference on Information Engineering and Computer Science, Wuhan, China, 19–20 December 2009; Institute of Electrical and Electronics Engineers (IEEE): Piscataway, NJ, USA, 2009; pp. 1–4. [Google Scholar] [CrossRef]
  22. Ben-Mansour, R.; Habib, M.; Shaik, A.Q. Modeling of Fluid Flow in a Tube with a Moving Indentation. Comput. Fluids 2009, 38, 818–829. [Google Scholar] [CrossRef]
  23. González, M.L.; Vega, M.G.; Oro, J.M.F.; Marigorta, E.B. Numerical Modeling of the Piston Effect in Longitudinal Ventilation Systems for Subway Tunnels. Tunn. Undergr. Space Technol. 2014, 40, 22–37. [Google Scholar] [CrossRef]
  24. Tijsseling, A.S.; Anderson, A. The Joukowsky Equation for Fluids and Solids; Eindhoven University of Technology: Eindhoven, The Netherlands, 2006. [Google Scholar]
Figure 1. Sequence of events during stoppage of the SID. (a) Moving SID inside pipe; (b) Extending external slips to decelerate; (c) Opening the valve to keep water flowing; (d) Launching another SID to the desired location; (e) Sealing the pipeline section between two SIDs.
Figure 1. Sequence of events during stoppage of the SID. (a) Moving SID inside pipe; (b) Extending external slips to decelerate; (c) Opening the valve to keep water flowing; (d) Launching another SID to the desired location; (e) Sealing the pipeline section between two SIDs.
Water 13 01794 g001
Figure 2. The moving SID in the pipeline (Non-simplified model).
Figure 2. The moving SID in the pipeline (Non-simplified model).
Water 13 01794 g002
Figure 3. The schematic of the simplified model.
Figure 3. The schematic of the simplified model.
Water 13 01794 g003
Figure 4. Characteristics lines in the x-t plane.
Figure 4. Characteristics lines in the x-t plane.
Water 13 01794 g004
Figure 5. Implicit method in the pipeline with a moving diameter-variation section representing the SID.
Figure 5. Implicit method in the pipeline with a moving diameter-variation section representing the SID.
Water 13 01794 g005
Figure 6. Node updates in dynamic meshes.
Figure 6. Node updates in dynamic meshes.
Water 13 01794 g006
Figure 7. The explicit–implicit coupling boundary condition in the pipeline.
Figure 7. The explicit–implicit coupling boundary condition in the pipeline.
Water 13 01794 g007
Figure 8. Flow chart of the transient simulation using dynamic meshes in a coupling method.
Figure 8. Flow chart of the transient simulation using dynamic meshes in a coupling method.
Water 13 01794 g008
Figure 9. Geometry of the 2D-axisymmetric pipe in CFD with a moving SID.
Figure 9. Geometry of the 2D-axisymmetric pipe in CFD with a moving SID.
Water 13 01794 g009
Figure 10. Independence analysis of the time step Δ t (s).
Figure 10. Independence analysis of the time step Δ t (s).
Water 13 01794 g010
Figure 11. Independence analysis of the space step size Δ l (m).
Figure 11. Independence analysis of the space step size Δ l (m).
Water 13 01794 g011
Figure 12. Comparison of the simulated piezometric heads at the back surface of the SID for Scenario 1.
Figure 12. Comparison of the simulated piezometric heads at the back surface of the SID for Scenario 1.
Water 13 01794 g012
Figure 13. The variation of head rises at point A with different clearance areas.
Figure 13. The variation of head rises at point A with different clearance areas.
Water 13 01794 g013
Figure 14. Comparison of the simulated piezometric heads at the back surface of the SID for Scenario 2.
Figure 14. Comparison of the simulated piezometric heads at the back surface of the SID for Scenario 2.
Water 13 01794 g014
Table 1. The properties and pipe parameters.
Table 1. The properties and pipe parameters.
PropertiesValues
Pressure head at Inlet H i n (m)240
Pressure head at Outlet H o u t (m)201.4
Initial discharge Q 0 (m3/s)0.024
Initial wave speed a (m/s)1000
Density of fluid ρ (kg/m3)998.2
Viscosity of fluid (Pa·s)1.003 × 10−3
Darcy–Weisbach friction factor f0.014
Gravitational acceleration g (m/s2)9.81
Initial velocity of the SID V (m/s)3
Length of Section S1 and S2 L 1 (m)300
Length of the SID L S (m)0.6
Length of Section S4 and S5 L 2 (m)300
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhang, K.; Zeng, W.; Simpson, A.R.; Zhang, S.; Wang, C. Water Hammer Simulation Method in Pressurized Pipeline with a Moving Isolation Device. Water 2021, 13, 1794. https://doi.org/10.3390/w13131794

AMA Style

Zhang K, Zeng W, Simpson AR, Zhang S, Wang C. Water Hammer Simulation Method in Pressurized Pipeline with a Moving Isolation Device. Water. 2021; 13(13):1794. https://doi.org/10.3390/w13131794

Chicago/Turabian Style

Zhang, Kang, Wei Zeng, Angus R. Simpson, Shimin Zhang, and Chao Wang. 2021. "Water Hammer Simulation Method in Pressurized Pipeline with a Moving Isolation Device" Water 13, no. 13: 1794. https://doi.org/10.3390/w13131794

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop