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

: 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 ﬁrst 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 efﬁciency 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). Veriﬁcations 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.


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. 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 finitedifference 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.

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.
conditions. However, the space step and time step of these methods are dependent an have to meet the Courant-Friedrichs-Lewy stability condition (CFL= a • ∆t/∆x) to achiev 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 explici 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 achieve even using a large space step. Another implicit method of characteristics (IMOC) was pro posed by Afshar and Rohani [11] with the simulation validated by the results from th MOC. Since the simulation using the implicit method can achieve unconditional conver gence and mutual independence between the space step and the time step, the space ste can be flexibly changed without changing the time step during the simulation. Thus, th implicit method is predicted to have the potential to be used with dynamic meshes (wher the mesh size changes with time step), which is essential to simulate a moving object i pipes.
To further improve the efficiency and robustness of the hydraulic transient simula tion, Wang and Yang [12] combined the explicit MOC with the Preissmann four-poin finite-difference implicit method [13] to simulate the unsteady flow in pipelines and th transient processes in hydropower systems. The results show that the coupling method i effective. This method, however, has only been applied to pipe systems without a movin object.
In this paper, a dynamic mesh technique in the implicit method has been develope to simulate the hydraulic transient process with a moving SID in a pipeline. The pipelin 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 bee used to obtain values at the updated mesh nodes. The explicit MOC has been also used t be coupled with the implicit method (referred to reference [12]) to improve the computin efficiency and the capability of handling complex boundary conditions. To validate thi novel simulation approach, numerical verifications using Computational Fluid Dynamic (CFD) and numerical simulation have been conducted. The comparison between the CFD simulation and the simulation using the new method illustrates the effectiveness of th proposed approach for the water hammer simulation of pipelines with a moving objec such as an SID.

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 movin cylinder. In order to achieve 1D simulation, the SID inside the pipeline has been simplifie to a diameter-reduced section, and the movement of the SID is considered to be equivalen to the change of the lengths of the pipe Section S2 and S4 as shown in the schematic i Figure 3.    The area (A3) of the simplified section in Figure 3 is assumed to be equal to the clearance area (A1 − A2) between the SID and the pipe wall when the SID is moving in Figure  2, as expressed by

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.

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]: in which = piezometric head; = discharge; = distance; = time; = wave speed; = gravitational acceleration; = Darcy-Weisbach friction factor; = pipe diameter; and = 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  Figure 3 is assumed to be equal to the clearance area (A 1 − A 2 ) between the SID and the pipe wall when the SID is moving in Figure 2, as expressed by

The Division of the Pipeline with Corresponding Methods
In order to simulate the water hammer pressures in pipelines with a moving diameterreduced 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.

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
∂H ∂x 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 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  (6) and (7) and H j i+1 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).
which are valid along the positive ( ) and negative ( ) characteristic lines, as shown in Figure 4. The subscripts i − 1, i and i + 1 represent the space steps, and the superscripts and j + 1 represent the time steps. By rearranging Equations (4) and (5), the following compatibility equations can be obtained According to the space-time plane shown in Figure 4, and at time step j + 1 can be calculated by using Equations (6) and (7) if , , and 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).

Discrete Implicit Method
The Preissmann four-point finite-difference implicit method (referred to as the "im plicit method" in the following) is introduced in this section. The discretization scheme o this implicit method with a temporal weighting coefficient is represented by [13] in which y can stand for the piezometric head or the flow rate. Replacing the partial de rivatives of Equations (2) and (3) by the above finite-difference scheme yields A further rearrangement of Equations (10) and (11) gives

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] 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 A further rearrangement of Equations (10) and (11) gives Water 2021, 13, 1794 All the coefficients are related to the heads and discharges at the time step j. ∆H and ∆Q are defined by 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 where EE 1 and FF 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 where EE i and FF 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.

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.
coefficients are related to the heads and discharges at the time step j. ∆ and ∆ fined by and denote the increments of piezometric head and discharge over the paramete previous time step. By combining Equations (12) and (13) with the boundary co the relationship of the discharge increment and piezometric head increment (re as the increment equation) at the first node can be defined as where and are two terms that are related to the boundary conditions, th metric head and the discharge at the first node at the time step j. The increment at the second node at time step j + 1 can be then obtained by incorporating Equat and (13) into (16). By repeating the processes, the increment equations at all the time step j+1 can be obtained as where and are two terms that are related to the boundary conditions, th metric head and the discharge at the node i at the time step j. This process is know forward scan of the chasing method [12]. By combining the increment equation a node with the boundary conditions at this node, the piezometric head and disc the end node can be obtained. Similarly, a backward scan of the chasing method w by combining the increment equation at node i with the solved piezometric head charge at node i+1 to solve the piezometric heads and discharges at all the nodes

Dynamic Meshes in the Implicit Method
The diameter-reduced Section S3 representing the SID as shown in Figure 5 the velocity of which is equal to the initial flow velocity. The size of the meshe tion S3 is fixed since the length of Section S3 remains unchanged during its mo The movement of Section S3, however, changes the lengths of Section S2 and S4 i 5, increasing for Section S2 and decreasing for Section S4. Lengths of S2 and S4 s longer than the moving distance of the SID during the deceleration process. B fixed mesh numbers for Section S2 and S4, the size of the meshes in these section be able to be altered during the simulation to enable the length change of Sectio 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 L 3 and L 4 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. The number of nodes of Section S2 and S4 are assumed to be and , 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 ∆ = ( − 1) ⁄ and ∆ = ( − 1) ⁄ with uniform meshes used in these sections. During one time step ∆ , the distance that S3 moves can be calculated as ∆ = ( + )∆ /2 with representing the velocity of the SID, and the lengths of Section S2 and S4 change to + ∆ and − ∆ , respectively. As a result, the mesh sizes of Section S2 and S4 need to be updated to ∆ = ( + ∆ )/( − 1) and ∆ = ( − ∆ )/( − 1). A simple case with only 4 nodes is presented in Figure 6. Since the distance moved of ∆ 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 can be expressed by The piezometric heads and discharges at the next time step + 1 can be then calculated using the parameters at the updated mesh nodes in the implicit method.

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 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  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 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.

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 S2,n 1 = EE S2,n 1 ·∆H S2,n 1 + FF S2,n 1 (20) 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, n 1 ) means the n 1th node of Section S2. At the connection node of Section S2 and S3 (Interface 2 in Figure 5) 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 S3 /A S2 .
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 Therefore, the increment equation at the first node in Section S3 can be then obtained by combining Equations (20), (21) and (23), which gives where 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 S3,n = EE S3,n ·∆H S3,n + FF S3,n where ξ 2 is the minor loss coefficient at Interface 3 for a sudden pipeline expansion and it can be calculated through The increment equation at the first node in Section S4 can be then obtained by substituting Equations (25) and (26) into (27), which gives where EE S4,1 = 1 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 Water 2021, 13, 1794 9 of 17 (shown in Figure 3) given in the following section and the range of the scan is reduced to Section S2, S3 and S4.

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: H j+1 S1,n 3 + Q j+1 S1,n 3 ·Q j+1 S1,n 3 2gA 2 S1,n 3 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.
Water 2021, 13, x FOR PEER REVIEW 9 of 1 (shown in Figure 3) given in the following section and the range of the scan is reduced t Section S2, S3 and S4.

Explicit-Implicit Coupling Method
In practice, there may be some complex boundary conditions at the inlet and outle of the pipe, which may be difficult to formulate using the implicit method. These bound ary conditions can be easily handled in the explicit MOC, which is numerically stable an computationally efficient. Based on the pipeline properties and the surrounding environ ment in the field, the unsteady friction loss in the pipe [14], the viscoelasticity of the pip wall [15,16], the fluid-structure interaction [17] and water-column separation [18] ma need to be considered, and the strategies to incorporate them into the MOC have bee extensively investigated and can be utilized if needed. Thus, the transient simulation o pipelines with a moving object inside can be extended to more complicated pipeline con figurations if it can be coupled with the MOC. To achieve that, the explicit-implicit cou pling boundaries have been applied and illustrated in the following way.
The explicit-implicit coupling boundaries are shown in Figure 7. The compat bility equation for the end node of Section S1, the continuity equation and the energ equation at Interface 1 as shown in Figure 7 can be written as follows: where is the minor loss coefficient for a sudden contraction between Section S1 and S2 In the case shown in Figure 7, = 0 since the diameters in these two sections are th same.  By combining Equations (30)-(32), the increment equation at the first node of Section S2 can be expressed as where EE S2,1 = −B, FF S2,1 = C P − B·H j S2,1 − Q j S2,1 . 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 S4,n 2 = EE S4,n 2 ·∆H S4,n 2 + FF S4,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. At Interface 4, Equation (34) can be combined with the continuity equation and energy equation at the interface and the 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.

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 = 0.1 m, the diameter of the SID is = 0.08 m. Thus the equivalent diameter of the simplified section in the 1D model in Figure 3 is = 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.

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. 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

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.

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 Fig  ure 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 simu lation model, the initial position of the SID is 150 m (3 m/s × 50 s) forward and thus = 150 m and = 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 ∆ of 0.001 m. The time step, which can be estimated by ∆ ≈ ∆ ⁄ 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 discretiza tion 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 sim ulation 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 move ment 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 cou pling 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 decel erate with a constant acceleration, and it stopped completely after a one-second decelera tion 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. Figures 10 and 11 show the differen 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 o the numerical precision, with deviations lower than 0.5%. 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. Figures 10 and 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%.

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 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 where = 4 with = ⁄ representing the hydraulic radius and 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

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 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 where = 4 with = ⁄ representing the hydraulic radius and 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

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 For Section S3, the friction loss in the 1D model is calculated based on the nonsimplified model as shown in Figure 2. Thus, it is where D e = 4R 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 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 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 To make the minor loss the same in both 1D and CFD simulations, the coefficient ϕ in Equation (38) can be calculated by 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.

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 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.
Water 2021, 13, x FOR PEER REVIEW 14 of 18 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] 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. To further validate the simulation, the Joukowsky Equation [24] 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.

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 f or 50 s < t ≤ 51 s 0 f or t > 51 s Water 2021, 13, x FOR PEER REVIEW 15 of 18 Figure 13. The variation of head rises at point A with different clearance areas.

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 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. 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.

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 Figures 5 and 7) are = Figure 14. Comparison of the simulated piezometric heads at the back surface of the SID for Scenario 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 (L E1 = L E2 = 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 Figures 5 and 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.

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.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.