Next Article in Journal
A Game for Learning How to Model in Graph Theory
Next Article in Special Issue
Frequency Analysis of the Railway Track under Loads Caused by the Hunting Phenomenon
Previous Article in Journal
A Symmetric/Asymmetric Bimodal Extension Based on the Logistic Distribution: Properties, Simulation and Applications
Previous Article in Special Issue
On Periods of Interval Exchange Transformations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Calculating Column Separation in Liquid Pipelines Using a 1D-CFD Coupled Model

1
Department of Civil Engineering, University of North Dakota, Grand Forks, ND 58202, USA
2
Innovative Hydraulic Group, 89 Loire Valley Ave, Thornhill, ON L4J 8V7, Canada
*
Author to whom correspondence should be addressed.
Mathematics 2022, 10(12), 1960; https://doi.org/10.3390/math10121960
Submission received: 1 May 2022 / Revised: 25 May 2022 / Accepted: 5 June 2022 / Published: 7 June 2022
(This article belongs to the Special Issue Theory and Application of Dynamical Systems in Mechanics)

Abstract

:
This paper proposes a coupled 1D-CFD model for calculating column separation in liquid pipelines. ANSYS Fluent is utilized to calculate two-phase flow analysis. Method of Characteristics and Discrete Gas Cavity Model (DGCM) are both employed to conduct 1D transient analysis. The results show that the proposed model, with both 2D and 3D CFD analysis, captures the transient responses of the system that have nearly identical accuracy and are both consistent with the results of a physical experiment. The results of a pure CFD analysis are employed to evaluate the performance of the proposed model in capturing the shape of the vapor cavity. The comparison shows that the results are generally consistent. However, the vapor cavity in the pure CFD model is established and grown on top of a film of liquid, while the cavity in some places fills the whole pipe cross-section in the proposed model. In addition, the results obtained from using Modified Two-Component Pressure Approach (MTPA), an open-channel based model proposed by the authors, also confirm the results obtained from the pure CFD analysis. Some minor discrepancies are found, which may be attributed to the uniform velocity distribution considered at the interface between 1D and CFD zones.

1. Introduction

Negative pressures may occur in liquid pipe systems under a variety of circumstances, examples of which include pumping power failure, hydropower load rejection, rapid closing of an inline valve, etc. [1,2]. When the pressure at a particular point of the pipe drops to the vapor pressure of the liquid, the liquid column locally ruptures, and the resulting vapor cavity causes the liquid to be separated. During the column separation, the pressure at the cavitating zone remains constant at the vapor pressure, allowing the liquid water columns on either side of the cavity to move independently and gain different velocities, whereby the cavity size changes [3]. Eventually, a positive pressure condition in the system causes the cavity to collapse and let the separated water columns rejoin. Rejoining the liquid column can induce significant water hammer pressure with an intensity that depends on the difference in the liquid column’s velocities as well as the acoustic wave speed of the pipe. The induced water hammer pressure can be strong enough to rupture the pipe and associated joints.
Due to its destructive nature, column separation has been the subject of a handful of research, both numerically [4,5,6,7,8,9,10,11,12] and experimentally [7,9,13,14,15], since it was first identified by Joukowsky (2006). Several numerical models have been proposed for calculating column separation, among which Discrete Vapor Cavity Model (DVCM) and Discrete Gas Cavity Model (DGCM) are the most popular approaches used in research and industry [3]. In DVCM, a vapor cavity is inserted in the computation node when the pressure drops to the vapor pressure. By keeping the pressure of the cavity constant at the vapor pressure, the size of the cavity at a given time is calculated by a simple mass balance approach. If the size of the cavity becomes nil or negative, the cavity is removed, and the liquid column rejoining is replicated. In DGCM, a permeant-infinitesimal gas cavity assumed in each computation allows for automatically capturing column separation and rejoining everywhere across the pipe length. Although the validity of these models has been widely approved by experimental studies, they cannot provide accurate results if a large cavity is extended over a long distance across a pipe with a steep slope. Khani et al. [11] presented a hypothetical test case in which the system is working under a steady-state-cavitating flow regime and showed that in such cases, extensive energy dissipation occurs in the cavitating zone that cannot be accounted for by either the DVCM or the DGCM. Khani also demonstrated that in this case, the DGCM captures a steady-state flow which is far off the true solution.
To improve the performance of DVCM and DGCM, different types of open-channel-based models have been proposed [16,17,18,19,20]. At each computational time level, the models utilize a shock-fitting approach to trace the interface between the cavitating flow and water hammer zones. Once the location of the interface is calculated, the flow parameters on either side of the interface are calculated using different sets of governing equations, the open channel flow equations for the cavitating zone, and the water hammer equations for the water hammer region. However, the results were not consistent with the experiments due to the challenge associated with calculating the location of the interface. Khani et al. [10] proposed the Modified Two-Component Pressure Approach (MTPA) that utilizes the open channel flow equations in both cavitating and water hammer flow regions. The model employs a shock-capturing approach which can automatically calculate both the interface between the water hammer and cavitating flow regions as well as the spread of the cavity across the pipe. The model results were shown to be in excellent agreement with both experiments and analytical solutions, implying that the MTPA can provide reasonable results even when the DVCM and DGCM provide poor results. However, when cavitating flow occurs in a highly three-dimensional flow field such as a large Turbine’s draft tube, the MTPA no longer provides reasonable results; in such cases, using Computational Fluid Dynamics (CFD) analysis appears to be inevitable.
CFD analysis has not received significant attention in the realm of column separation, probably because 1D models can alternatively cover a wide range of applications. Nevertheless, a few researchers tried to investigate the performance of CFD in capturing the key aspects of cavitating flow and column separation in pipe systems. Tang et al. [21] and Wang et al. [22] utilized 2D CFD analysis for calculating column separation that occurs following the rapid closing of a valve in a simple reservoir-pipe-valve system and found that not only can CFD accurately capture the time history of the resulting transient pressures, but it also can successfully account for the spread of the vapor cavity across the pipe. Warda et al. [23] employed a 3D CFD analysis for calculating column separation in a laboratory-scale pipe system and concluded that CFD can capture key features of cavitating flow and column separation quite accurately. However, CFD analysis is quite expensive, at least at present, in terms of computational resources. Wang et al. [22] demonstrated that it takes a few days for a personal computer to complete 0.6 s of transient analysis in a pipe with an internal diameter and length of 0.02 and 30 m, respectively. This implies that CFD is inefficient in practice due to the scale of pipe systems and the fact that an engineering design is iterative and may need numerous simulations.
Coupling CFD with 1D analysis appears to be an efficient solution to significantly reduce the computational time. This idea is supported by the fact that mass, momentum, and energy are dominantly transferred in the longitude dimension of the pipe, and the information exchanged in the radial dimension of the pipe is not of insignificant importance. In CFD-1D analysis, the limited part of the system with the significant 3D feature is calculated by the CFD analysis, and the rest of the system is handled by the 1D analysis; the two computational approaches exchange the information at each computational time level through the interface connecting 1D and CFD domains. This method has been successfully employed for calculating transient flow in liquid pipe systems [24]. For example, Zhang and Cheng [25] applied 1D-3D analysis to different components of hydropower systems. Zhang et al. [26] successfully utilized 1D-3D analysis in transient analysis of a pumped-storage system. Wu et al. [27] applied 1D-3D coupling for steady-state and transient flow analysis in a pumped pipeline. Maddahian et al. [28] employed 1d-3D analysis in rapid pressurization of a pipe system with an entrapped air pocket and showed that the modeling results are consistent with experiments. Different CFD-1D coupling analysis approaches, including the Picard-based and Newton-based methods, have been successfully utilized in other realms of engineering [29,30,31]. Nevertheless, to the best of the authors’ knowledge, the 1D-3D coupling approach has not been utilized for calculating column separation so far.
Thus, this paper aims to shed light on the performance of 1D-CFD coupling in simulating liquid column separation analysis in pipe systems including water pipe systems. The organization of this paper is as follows: the theoretical background is first presented, and the numerical results are then discussed. Finally, the paper is summarized, and conclusions are drawn.

2. Theoretical Backgrounds

2.1. DGCM

Transient flows in liquid pipe systems are governed by the following 1D equations, which are also known as water hammer equations [2]:
V t + g H x + f 2 D | V | V = 0
H t + a 2 g V x = 0
where V = velocity, H = piezometric head, D = pipe diameter, a = elastic wave velocity, g = gravitational acceleration, f = friction factor in the Darcy–Weisbach equation, and x, t = distance and time independent variables, respectively.
These equations hold as long as the pressure exceeds the vapor pressure of the liquid. Otherwise, cavitating flow and column separation occur, and the equations fail to replicate the physics of the flow. There are several approaches to include the column separation feature in the water hammer equations, among which DGCM is the more popular and accurate one. In DGCM, a small gas pocket assumed in each computational cell can account for column separation in the system.
To numerically implement DGCM, the well-known Method of Characteristics (MOC) is utilized. The first step in the MOC is to discretize the computational domain into numerous computational cells, as shown in Figure 1, with identical spatial and temporal increments, Δ x and Δ t , respectively. The unknown flow parameters at the computational points at time t n + 1 can then be explicitly calculated based on the information from the previous timeline.
The four unknowns at a given point (P) include the flow rates on either side of the gas pocket, the piezometric head, and the volume of the gas pocket and can be calculated by the following discretized equations, which are called positive and negative characteristic equations, continuity equation and the gas equation of state.
C + :     H P u = H P d = C 1 C 2 Q P u
C :   H P d = C 3 + C 4 Q P d
( P ) t + Δ t = ( P ) t + [ ψ ( Q P d Q P u ) t + Δ t + ( 1 ψ ) ( Q M d Q M u ) t ] × Δ t
( P ) t + Δ t = P 0 * α 0 R ρ g ( H P d Z P d H V )                      
with
C 1 = H L d + a g A Q L d   ;     C 2 = a g A + f Δ x 2 g D A 2 | Q L d |
C 3 = H R u a g Q R u ;   C 4 = a g A + f Δ x 2 g D A 2 | Q R u |
where Q P u , Q P d = flow rates at the upstream and downstream side of the gas pocket at point P , respectively, H P u = H P d = piezometric head at the upstream and downstream side of the gas pocket at point P , respectively, Z P u = Z P d = pipe elevation upstream and downstream side of the gas pocket at point P , respectively, Q L d = flow rate at the downstream side of the gas pocket at point L , H L d = piezometric head at the downstream side of the gas pocket at point L , Q R u = flowrate at the upstream side of the gas pocket at point R , H R u = piezometric head at the upstream side of the gas pocket at point R , A = pipe cross sectional area, D = pipe diameter, ρ = water density, P = gas pocket volume at point P , H V = gage vapour pressure head, α 0 = initial gas void fraction, R = computational cell volume, P 0 * = reference pressure, Δ t = computational time step, S 0 = pipe slope, S f = energy slope, g = gravitational acceleration, and t ,   t + Δ t   = indexes refer to the previous and current time lines, respectively.
By combining Equations (3)–(6) and doing some algebra, the following equation is reached:
( H P d Z P d H V ) 2 + 2 B 1 ( H P d Z P d H V ) B 4 = 0
The analytical solution of the above equation is as follows.
H P d Z P 1 d H V = B 1 ( 1 + 1 + B B )     i f     B 1   0
H P d Z P 1 d H V = B 1 ( 1 1 + B B )     i f     B 1 > 0
where
B 1 = B 2 C 4 C 2 B V + 0.5 ( Z P d + H V ) B 2 ( C 2 C 3 + C 4 C 1 ) ;       B B = B 4 B 1 2
B 2 = 0.5 C 2 + C 4 ;     B V = ( P ) t ψ Δ t + 1 ψ ψ ( Q P d Q P u ) t ;     B 4 = P 0 * α 0 R B 2 C 4 C 2 0.5 ρ g ψ Δ t
In cases with high gas volume and very low pressure or very low volume but high pressure, indicated by | B B | 1 , Equations (8) or (9) may provide inaccurate results due to miscalculation of the radical. The appropriate results in such cases can be achieved through the linearization of the original equations [2]:
H P d Z P d H V = 2 B 1 B 4 2 B 1                           i f     B 1   0
H P d Z P d H V = B 4 2 B 1                                                         i f     B 1   > 0
Having H P d calculated, Q P u , Q P d , and P can be calculated by Equations (3)–(5), respectively.

2.2. MTPA

The Two-Component Pressure Approach (TPA) was proposed by Vasconcelos et al. [32] to address the limitation of the Preissmann Slot Method (PSM) in capturing negative pressures [33,34]. In this method, the conservative form of the equations that govern unsteady flow in open channels are rearranged as follows:
A t + Q x = 0                                                                                                                                     Continuty   Equation  
Q t + ( Q 2 A + A g [ h c + h s ] ) x = g A ( S 0 S f )                         Momentum   Equation
where A   = flow cross sectional area, Q = flow rate, h c = the distance between the free surface and the centroid of the flow cross-sectional, h s = pressure head, g = gravitational acceleration, S 0 = pipe slope, and S f = frictional slope, which can be calculated by the following equation:
S f = Q 2 n m 2 A 2 R h 4 3
where n m = Manning coefficient and R h = hydraulic radius of the flow.
Splitting the pressure term in the momentum equation allows TPA to account for both pressurized and open channel flows, with h s representing the pressure of the conduit in pressurized flow and h c measuring the flow depth in the open channel flow regime. In TPA, the pressure head and pipe acoustic wave speed are related by the following equation [29]:
h s = [ A A f A f ] a 2 g
where a = acoustic wave speed, and A f = pipe cross-sectional area.
Many researchers have successfully utilized the TPA for capturing mixed-flow and negative pressures in close conduit systems [32,35,36,37,38]. However, TPA still cannot replicate the physics of the flow when column separation occurs in the system. In such cases, the TPA generates negative pressures that are far below the vapour pressure of the liquid, which is not physically sound.
To fill this gap, Khani et al. [11] proposed a modified version of the TPA (MTPA) which can account for cavitating flow and column separation. In the MTPA, Equations (1) and (2) can calculate the transient pressure heads and flow rates across the system as long as the pressures are above the vapor pressure. Otherwise, Equation (2) no longer holds, and the following equation should be utilized instead.
Q t + ( Q 2 A + g [ A h c + A f h v ] ) x = g A ( S 0 S f )
where h v = vapor pressure head of the liquid.
To solve the governing equations, the Godunov scheme is utilized in this research. The Godunov scheme is a finite volume numerical scheme widely used in solving hyperbolic partial differential equations [39,40]. The spatial domain in this approach is broken down into equal size computational cells with a spatial distance of Δ x , and the temporal domain is discretized with a constant time step, Δ t . The first-order accuracy of the scheme suits the current application as the flow transition from open channel to pressurized flow may induce significant spurious numerical oscillations that can be better suppressed by a first-order accuracy scheme. Applying piecewise constant data reconstruction at the computational cells results in first-order accuracy.
By discretizing the governing equations, unknowns at the current time level can be explicitly calculated based on the data retrieved from the previous timeline using the following equation:
A i n + 1 = A i n Δ t Δ x ( F i + 1 2 n F i 1 2 n )
Q i n + 1 = Q i n Δ t Δ x ( G i + 1 2 n G i 1 2 n ) + Δ t g A   ( S 0 S f i n )
where i is computational cell number, i + 1 2 and i 1 2 refer to the downstream and upstream boundaries of the i t h cell, respectively, n and n + 1 refer to the previous and current timelines, respectively, and F and G are mass and momentum fluxes, respectively.
In the Godunov scheme, the mass and momentum fluxes at cell boundaries are calculated by solving the Riemann solution [39,40]. The exact Riemann solution can be obtained through an iterative procedure, but as shown by Malekpour and Karney [41], it can produce significant numerical oscillation during flow transition. Khani et al. [38] adopted the HLL solver that can automatically increase the numerical viscosity of the scheme during the flow regime transient and showed that the scheme provides a non-oscillatory solution at a wide range of acoustic pipe speeds.
The fluxes in the HLL solver are generally calculated through a procedure summarized in Equation (19).
Γ * = Γ L       i f   S L > 0 Γ * = S R Γ L S L Γ R + S L S R ( U R U L ) S R S L         i f   S L 0   a n d   S R 0   Γ * = Γ R       i f   S R < 0
where S L and S R are the left and right shock wave speeds, respectively, U L = [ A L Q L ] and U R = [ A R Q R ] are dependent variables at the computational cells i 1 and i + 1 , respectively, and Γ L = [ F L G L ] and Γ R = [ F R G R ] are the fluxes at the computational cells i 1 and i + 1 , respectively.
Khani et al. [37,38] showed that the numerical viscosity admitted to the numerical scheme is controlled by the magnitude of the right and left wave velocities. A simple strategy for calculating the left and right wave speeds in both water hammer and open channel flow regions was proposed. The left and right wave speeds are calculated using the following equations:
S L = V L Ω L   ;   S R = V R + Ω R
Ω K ( L , R ) = g [ Y G A G ( h C + | h s | ) K A K ] A G A K ( A G A K )
where variables with sub-index G are the function of Y G that needs to be estimated.
As shown by Malekpour and Karney [41], Equation (21) significantly increases wave speeds only when the water depth is close to the pipe crown. By using this equation, the artificial viscosity of the scheme is increased whenever required, i.e., when the flow transient is proximate. Khani et al. [37] found that the optimal amount of artificial viscosity is admitted to the scheme if Y G is calculated by the following equation:
Y G = K a × M A X   [ d i N S , d i N S + 1 , , d i , d i + 1 , d N S ]
where d = h c + | h s | .
Equation (22) ensures that artificial viscosity is distributed within NS computational cells on either side of the cell for which the wave velocity is calculated. Numerical exploration reveals that N S should be selected such that the numerical viscosity is spread over a length 3 or 4 times higher than the pipe height. In any case, the NS should not be less than 3. Numerical exploration also reveals that in the computational cells carrying free surface flow, K a can take a value of around 1.4 if a pressurization front exists within N S computational cells on either side of the target computational cell. Otherwise, N S can be considered 1.001.

2.3. CFD Analysis

The ANSYS Fluent is utilized to perform CFD analysis in the portion of the pipeline experiencing cavitating flow and column separation. Fluent can simulate two-phase flow by numerous approaches, among which the Volume of Fluid (VOF) is widely used for capturing cavitating flow [42]. In the VOF method, a single set of Unsteady Random-averaged Navier-Stokes (URANS) equations presented below are solved for the mixture of vapor and liquid with a weighted average density and viscosity of ρ m and μ m , respectively.
ρ m t + x i ( ρ m u i ) = 0
t ( ρ m u i ) + x j ( ρ m u i u j ) = P x i + x j [ μ m ( u i x j + u j x i 2 3 δ i j u k x k ) ]   + x j ( ρ m u i u j ¯ ) + ρ m g i
with
ρ m = α v ρ v + ( 1 α v ) ρ l
μ m = α v μ v + ( 1 α v ) μ l
where u i = flow velocity vector, P = pressure, ρ v = vapor density, ρ l = liquid density,   μ v = vapor viscosity, μ l = liquid viscosity, δ i j = Kronecker delta, α v = vapor void fraction, g i = gravitational acceleration, x j = spatial independent vector, and t = temporal independent variable.
The above equations are not in closed form as the time averaging included an extra term, ρ m u i u j ¯ , which is well known as the Reynolds Stresses term. A few common methods used for approximating the Reynolds stresses are based on the Boussinesq hypothesis that relates the Reynolds stresses to the mean velocity gradients [42]. In this research, the standard k ε method is employed to calculate the Reynolds stresses. Since only the velocity and pressure are communicated at the boundary, the turbulence information at the boundary is not changed during the transient analysis, and it is assumed that it remains constant at what is set at the initial condition.
Mass transfer between liquid and vapor phases during cavitating flow is calculated by including the mass transfer terms on the right-hand side of the following partial differential equation utilized for tracking the interface of vapor and liquid.
t ( α ρ v ) + x i ( α ρ v u v i ) = R e R c
where R e and R c are mass transfer source terms connected to the growth and collapse of the vapor bubbles, respectively.
In ANSYS Fluent, the Singhal et al., Zwart–Gerber–Belamri, and Schnerr and Sauer models can be used for calculating the source terms in Equation (25), the latter of which is employed in this study. The Schnerr and Sauer model calculates the source terms using the following equations:
R e = ρ v ρ l ρ m α ( 1 α ) 3 B 2 3 ( P v P ) ρ l
when P v     P
R c = ρ v ρ l ρ m α ( 1 α ) 3 B 2 3 ( P P v ) ρ l
where B = radius of the bubbles and P v = vapor pressure of the liquid.
The vapor void fraction and radius of the bubbles can be calculated by the following equations:
α = n b 4 3 π B 3 1 + n b 4 3 π B 3
B = ( α 1 α 3 4 π 1 n b ) 1 3  
The above equations show that both vapor void fraction and the radius of the bubbles depend on the number of spherical bubbles per volume of liquid ( n b ). Assuming that no bubbles are created or destroyed, the initial conditions for the nucleation site volume fraction and the equilibrium bubble radius would therefore be sufficient to specify the bubble number density ( n b ).
It is worth noting that to accurately calculate column separation and resulting water hammer pressures, the liquid phase needs to be considered compressible. The compressibility of the liquid is included through a User-Defined Function (UDF). The UDF function relates the density of the liquid to its pressure with the aid of the following equation:
ρ l = ρ r e f ( 1 + P P r e f K e q )
where ρ r e f = density of the liquid at a reference pressure ( P r e f ) and K e q = equivalent bulk modulus elastic of the liquid that can be calculated based on the real acoustic speed of the pipe by the following equation:
K e q = a 2 ρ r e f

2.4. CFD Coupling

In 1D-CFD analysis, the surface(s) separating the CFD and 1D computational zones are the places where the data are exchanged between the CFD and 1D models. In this study, 1D-CFD modeling is carried out in the context of a simple pipe system in which a vapor cavity is concentrated at a local point. Therefore, a small part of the pipe in which the cavity is expanded and contracted is simulated by CFD and the rest is treated with the aid of MOC and DGCM; the data are exchanged between 1D and 3D computational zones at a single pipe section during the simulation.
There are a handful of approaches presented for exchanging data between 1D and CFD computational regions [24], but they can be generalized into two categories: partially and fully coupled methods. In the partially coupled method, the data are exchanged in every time step, while in the other one it is done in every iteration made by CFD toward numerical convergence. The former method can provide numerically stable results as long as the shared flow parameters do not change significantly in each time step or when a very small time step is used. Since strong water hammer pressures would be exchanged at the interface in the current problem, the latter method is utilized.
To facilitate the explanation of the fully coupled method, consider Figure 2. The simulation starts from a steady-state flow condition in which both CFD and 1D flow regions agree on both flow velocity and pressure at the interface. In each CFD iteration toward convergence, the total pressure calculated at the interface in the CFD zone is employed to calculate the resulting new flow velocity of the 1D using the positive characteristic equation. This method can produce some over- and undershooting around the solution and the convergence is therefore hardly achieved. A relaxation method proposed by Wu et al. [27] is used to resolve the convergence issue. In this method, the weighted average of the velocities obtained in two subsequent iterations is used to update the velocity at the CFD boundary. As shown in Figure 3, the convergence is assumed to be reached when the difference between the velocities (e) becomes less than a predefined value. Note that by reducing the relaxation factor ( R f ), a better convergence condition is achieved but at the expense of more iterations.

3. Numerical Results

The 1D-CFD analysis is performed in the context of a pipe system utilized by Bergant and Simpson [9] for experimenting with column separation. As shown in Figure 4, the test rig consists of a copper pipeline with an internal diameter, length, wave velocity, and upward slope of 22.1 mm, 37.23 m, 1319 m/s, and 5.6%, respectively. The pipe connects two pressurized tanks with constant pressures. Column separation is induced by closing the upstream valve in 0.009 s when the system is running in a steady-state condition at a flow velocity of 1.5 m/s. During transient, the pressure head at Tank 2 is kept constant at 22 m.
Warda et al. [23] utilized a 3D-CFD analysis to calculate column separation in the above-described system and showed that a local vapor cavity would form at the upstream end of the pipe right after the valve. This local vapor cavity was shown to remain in the upstream 30 cm of the pipe. Thus, in the 1D-CFD coupling, only the upstream 50 cm of the pipe is considered the CFD zone, and the rest (36.73 m) is analyzed by the DGCM. Numerical explorations show that a time step = 0.0001 s provides a numerical stable condition when the CFD zone is covered with a fine mesh with an average size of 1 mm. Both 2D- and 3D-CFD analyses use the same quadrilateral mesh type and 40,000 and 600,000 computational cells, respectively. The number of computation cells in the 1D analysis is selected in such a way that (1) the time steps in both 1D and CFD become identical, and (2) the Courant number in the 1D zone remains at 1. This results in 278 computational cells with a spatial increment length of 0.1321 m. Note that when the Courant number approaches 1, the interpolation error tends toward zero and the MOC provides quite accurate results [43].
The 1D-CFD analysis is performed for two distinct cases: 2D- and 3D-CFD analysis. A UDF is developed to linearly reduce the flow velocity at the upstream face of the CFD zone from 1.5 m/s to 0 in 0.009 s. Figure 5 and Figure 6 compares the pressure time histories obtained from both the DGCM and the 1D-2D coupling method with the experimental results at the valve location and the midpoint of the pipe. As can be seen, the results are consistent with the experiments. Figure 7 also demonstrates the evolution of the vapor cavity in different computational times.
Figure 8 compares the pressure time histories obtained from the 3D-1D coupling method with the experimental results at the valve location and the midpoint of the pipe. As can be seen, the numerical results are consistent with the experiment. Comparing Figure 8 and Figure 6 reveals that both the 1D-2D and 1D-3D coupling approaches provide almost identical results. Figure 9 also shows the evolution of the vapor cavity in different computational times. Comparing Figure 9 and Figure 7 reveals that the shapes of the cavity captured by the 1D-2D and 1D-3D methods are almost identical.
Although the experimental data confirm the validity of the results captured by 1D-CFD coupling, it is not known if 1D-CFD coupling succeeded in capturing the shape of the vapor cavity. The CFD results obtained from a CFD analysis are utilized to verify this aspect of the model. Warda et al. [23] captured the evolution of the cavity for the current problem using 3D CFD analysis. Figure 10 shows a few snapshots of the cavity shape during its contraction and expansion. Note that Figure 10 shows the first 30 cm of the pipe while Figure 7 and Figure 9 represent the whole CFD domain, which is 50 cm. Comparing Figure 7 and Figure 9 with Figure 10 reveals that the shape of the cavity in 1D-CFD analysis differs from what is captured by the CFD analysis. The figures demonstrate that in the CFD analysis, the cavity is formed on top of a film of liquid, but in the 1D-CFD analyses, the vapor cavity partially occupied the whole section of the pipe. In the absence of an experimental study, it is not easy to identify which one better represents reality. Nevertheless, the physics of the flow shows that the results from the CFD analysis can better represent reality. Due to the gravity effect, the pipe obvert receives higher pressure than the bottom, so a vapor cavity starts forming at the top and growing toward the bottom. This explanation may explain why the vapor cavity grows over a film of liquid.
An independent way of justifying the shape of the vapor cavity is to simulate the current problem with the aid of the MTPA. Figure 11 and Figure 12 show the pressure time histories and the cavity evolution calculated by this method. As shown in Figure 11, the pressure time histories calculated by MTPA match the experiment. Figure 12 also reveals that like pure CFD results, the cavity is expanding and growing over a film of liquid, which is another justification for the CFD’s ability to represent the cavity evolution better than the coupled 1D-CFD model. That the 1D-CFD model distorts the shape of the cavity may be attributed to the issue of considering a uniform velocity at the interface between the CFD and 1D domains. In reality, the velocity distribution is not uniformly distributed on the pipe cross-section, and this may explain the discrepancy between the CFD and coupled CFD-1D models. This is of course just a hypothesis that is subject to verification.
It is worth mentioning that the proposed coupled 1D-CFD analysis can significantly decrease the computational time. For the test case considered in this paper, it takes 4 h for the coupled 1D-CFD model to complete 1.5 s of transient flow, while a full CFD analysis may take around a few days to perform the same analysis. Since, in practice, the design is iterative, the proposed model can help the designer present the state-of-the-CFD knowledge in a reasonable time. Nevertheless, it should be pointed out that although the proposed model provides reasonable accuracy, further investigation is required to shed light on how the scale of the system affects the accuracy of calculations.

4. Summary and Conclusions

In this section, ANSYS Fluent is utilized to calculate cavitating flow and column separation, and the Method of Characteristics and the Discrete Gas Cavity Model (DGCM) are employed to conduct 1D transient analysis. The simulations are conducted in the context of a lab-scale pipe system for which experimental data are available. The small portion of the pipe containing the vapor cavity is treated by CFD analysis, while the rest of the system is analyzed by the 1D model. The model uses a fully coupled approach in which the data are exchanged in every CFD iteration through the pipe section separating 1D and CFD zones. A relaxation approach is utilized to ensure the convergence of the CFD analysis in each time step. Both 2D- and 3D-CFD analyses are independently employed to explore to what extent the 3D effect can affect the solution.
The results show that the proposed model with both 2D- and 3D-CFD analysis captures the transient responses of the system that have almost identical accuracy are are both consistent with the experiment. Since the evolution of the cavity has not been experimented with, the results of a pure CFD analysis published in the literature are employed to evaluate the performance of the proposed model in capturing the shape of the vapor cavity. The comparison shows that although the results are in general agreement, there is still some discrepancy that needs to be explained. The vapor cavity in the pure CFD model is established and grown on top of a film of liquid, while in the proposed model the cavity in some places fills the whole pipe cross-section. As a further test, the proposed model results are compared to the results obtained from the Modified Two-Component Pressure Approach (MTPA), a novel open channel-based 1D model recently proposed by the author for calculating column separation. The MTPA results also confirm the results from the CFD analysis, according to which the cavity is formed and evolved on top of a film of liquid. What makes the proposed coupled 1D-CFD model distort the shape of the cavity may be attributed to the way the data exchange between the 1D and CFD models. In each time step, the flow velocity calculated by the 1D model is applied to the CFD boundary with a uniform distribution. However, in reality, the velocity distribution at that section of the pipe is not uniform and this may affect the shape of the cavity. Nevertheless, this is just a hypothesis that warrants a detailed investigation.
Finally, the proposed model has great potential to be used in the context of a wide range of engineering problems, including transient analysis in hydropower systems and analysis of pipe collapse, which may result in extensive negative pressure and column separation.

Author Contributions

D.K. conducted the research and prepared the manuscript of the paper, while the Y.H.L. and A.M. reviewed the paper and provided their editorial and technical comments. All authors have read and agreed to the published version of the manuscript.

Funding

This material is based upon work supported by the U.S. Geological Survey under Grant/Cooperative No. (G16AP00075).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The experimental data used in this paper are part of the research conducted by Bergant and Simpson (1999, and 1999a).

Acknowledgments

The authors would like to thank Anton Bergant, who generously provided us with the experimental data utilized for the validation of the proposed model.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Chaudhry, M. Applied Hydraulic Transients; Van Nostrana Reinhold Co.: New York, NY, USA, 1987. [Google Scholar]
  2. Wylie, E.B.; Streeter, V.L. Fluid Transients in Systems; Prentice-Hall: Upper Saddle River, NJ, USA, 1993. [Google Scholar]
  3. Bergant, A.; Simpson, A.R.; Tijsseling, A.S. Water hammer with column separation: A historical review. J. Fluids Struct. 2006, 22, 135–171. [Google Scholar] [CrossRef] [Green Version]
  4. Wylie, E. Simulation of vaporous and gaseous cavitation. ASME J. Fluids Eng. 1984, 106, 307–311. [Google Scholar] [CrossRef]
  5. Bergant, A.; Simpson, A.R. Range of Validity of the Discrete Vapour Cavity Model for Pipeline Water-Hammer; University of Adelaide, Department of Civil and Environmental Engineering: Adelaide, Australia, 1994. [Google Scholar]
  6. Simpson, A.; Bergant, A. Numerical comparison of pipe column-separation models. J. Hydraul. Eng. 1994, 120, 361–377. [Google Scholar] [CrossRef] [Green Version]
  7. Provoost, G.A.; Wylie, E.B. Discrete gas model to represent distributed free gas in liquids. In Proceedings of the 5th International Symposium on Water Column Separation, Obernach, West Germany, 1981; Volume 1, pp. 28–30. [Google Scholar]
  8. Adamkowski, A.; Lewandowski, M. Investigation of hydraulic transients in a pipeline with column separation. J. Hydraul. Eng. 2012, 138, 935–944. [Google Scholar] [CrossRef]
  9. Bergant, A.; Simpson, A.R. Pipeline column separation flow regimes. J. Hydraul. Eng. 1999, 125, 835–848. [Google Scholar] [CrossRef] [Green Version]
  10. Khani, D.; Lim, Y.H.; Malekpour, A. An Innovative Open Channel Based Model for Calculating Column Separation in Conduit Systems. J. Hydraul. Eng. 2022; under review. [Google Scholar]
  11. Khani, D.; Lim, Y.H.; Malekpour, A. Investigating the Performance of the Modified Two-Component Pressure Approach in Capturing Column Separation with Large Vapor Cavities. J. Hydraul. Res. 2022; under review. [Google Scholar]
  12. Vasconcelos, G.J.; Marwell, T.B. Innovative simulation of unsteady low-pressure flows in water mains. J. Hydraul. Eng. 2011, 137, 1490–1499. [Google Scholar] [CrossRef] [Green Version]
  13. Simpson, A.; Wylie, E. Large water-hammer pressures for column separation in pipelines. J. Hydraul. Eng. 1991, 117, 1310–1316. [Google Scholar] [CrossRef] [Green Version]
  14. Martin, C. Experimental investigation of column separation with rapid closure of downstream valve. In Fourth International Conference on Pressure Surges; BHRA Fluid Engineering: Cranfield, UK, 1983; pp. 77–88. [Google Scholar]
  15. Autrique, R.; Rodal, A.E. Laboratory studies of water column separation. IOP Conf. Ser. Mater. Sci. Eng. 2013, 52, 022022. [Google Scholar] [CrossRef] [Green Version]
  16. Baltzer, R.A. Column separation accompanying liquid transients in pipes. J. Basic Eng. 1967, 89, 837–846. [Google Scholar] [CrossRef]
  17. Baltzer, R. A Study of Column Separation Accompanying Transient Flow of Liquids in Pipes. Ph.D. Thesis, University of Michigan, Ann Arbor, MI, USA, 1967. [Google Scholar]
  18. Siemons, J. The phenomenon of cavitation in a horizontal pipe-line due to a sudden pump-failure. J. Hydraul. Res. 1967, 5, 135–152. [Google Scholar] [CrossRef]
  19. Kalkwijk, J.; Kranenburg, C.; Vreugdenhil, C.; De Vries, A. Cavitation Caused by Water Hammer in Horizontal; Delft Hydraulics Laboratory: Delft, The Netherlands, 1972; Publication No. 97. [Google Scholar]
  20. Marsden, N.; Fox, J. An alternative approach to the problem of column separation in an elevated section of pipeline. In Second International Conference on Pressure Surges; BHRA: London, UK, 1976; pp. 1–13. [Google Scholar]
  21. Tang, X.; Duan, X.; Gao, H.; Li, X.; Shi, A.X. CFD investigations of transient cavitation flows in pipeline based on weakly-compressible model. Water 2020, 12, 448. [Google Scholar] [CrossRef] [Green Version]
  22. Wang, H.; Zhou, L.; Liu, D.; Karney, B.; Wang, P.; Xia, L.; Ma, J.; Xu, C. CFD Approach for Column Separation in Water Pipelines. J. Hydraul. Eng. 2016, 142, 04016036. [Google Scholar] [CrossRef] [Green Version]
  23. Warda, H.A.; Wahba, E.M.; El-Din, A.M. Computational Fluid Dynamics (CFD) simulation of liquid column separation in pipe transients. Alex. Eng. J. 2020, 59, 3451–3462. [Google Scholar] [CrossRef]
  24. Mandair, S. 1D and 3D Water-Hammer Models: The Energetics of High Friction Pipe Flow and Hydropower Load Rejection. Ph.D. Thesis, University of Toronto, Toronto, ON, Canada, 2020. [Google Scholar]
  25. Zhang, X.-X.; Cheng, Y.-G. Simulation of hydraulic transients in hydropower systems using the 1-D-3-D coupling approach. J. Hydrodyn. 2012, 24, 595–604. [Google Scholar] [CrossRef]
  26. Zhang, X.X.; Cheng, Y.G.; Xia, L.S.; Yang, J.D. Dynamic characteristics of a pump-turbine during hydraulic transients of a model pumped-storage system: 3D CFD simulation. IOP Conf. Ser. Earth Environ. Sci. 2014, 22, 032030. [Google Scholar] [CrossRef] [Green Version]
  27. Wu, D.; Yang, S.; Wu, P.; Wang, L. MOC-CFD coupled approach for the analysis of the fluid dynamic interaction between water hammer and pump. J. Hydraul. Eng. 2015, 141, 06015003. [Google Scholar] [CrossRef]
  28. Maddahian, R.; Shaygan, F.; Bucur, D.M. Developing a 1D-3D model to investigate the effect of entrapped air on pressure surge during the rapid filling of a pipe. IOP Conf. Ser. Earth Environ. Sci. 2021, 774, 012069. [Google Scholar] [CrossRef]
  29. Cui, X.; Wong, L.N. A 3D fully thermo–hydro–mechanical coupling model for saturated poroelastic medium. Computer Methods Appl. Mech. Eng. 2022, 394, 114939. [Google Scholar] [CrossRef]
  30. Failer, L.; Richter, T. A parallel Newton multigrid framework for monolithic fluid-structure interactions. J. Sci. Comput. 2020, 82, 28. [Google Scholar] [CrossRef] [Green Version]
  31. Zhang, H.; Guo, J.; Lu, J.; Niu, J.; Li, F.; Xu, Y. The comparison between nonlinear and linear preconditioning JFNK method for transient neutronics/thermal-hydraulics coupling problem. Ann. Nucl. Energy 2019, 132, 357–368. [Google Scholar] [CrossRef]
  32. Vasconcelos, J.G.; Wright, S.J.; Roe, P.L. Improved simulation of flow regime transition in sewers: Two-component pressure approach. J. Hydraul. Eng. 2006, 132, 553–562. [Google Scholar] [CrossRef]
  33. Vasconcelos, J.G.; Wright, S.J. Comparison between the two-component pressure approach and current transient flow solvers. J. Hydraul. Res. 2007, 45, 178–187. [Google Scholar] [CrossRef]
  34. Sanders, B.F.; Bradford, S.F. Network implementation of the two-component pressure approach for transient flow in storm sewers. J. Hydraul. Eng. 2010, 137, 158–172. [Google Scholar] [CrossRef]
  35. Cunge, J.A.; Wegner, M. Numerical Integration of Bane de Saint-Venant’s Flow Equations by Means of an Implicit Scheme of Finite Differences. Applications in the Case of Alternately Free and Pressurized Flow in a Tunnel. La Houille Blanche 1964, 1, 33–39. [Google Scholar] [CrossRef] [Green Version]
  36. Kerger, F.; Archambeau, P.; Erpicum, S.; Dewals, B.J.; Pirotton, M. An exact Riemann solver and a Godunov scheme for simulating highly transient mixed flows. J. Comput. Appl. Math. 2011, 235, 2030–2040. [Google Scholar] [CrossRef] [Green Version]
  37. Khani, D.; Lim, Y.H.; Malekpour, A. Hydraulic Transient Analysis of Sewer Pipe Systems Using a Non-Oscillatory Two-Component Pressure Approach. Water 2020, 12, 2896. [Google Scholar] [CrossRef]
  38. Khani, D.; Lim, Y.H.; Malekpour, A. A Mixed Flow Analysis of Sewer Pipes with Different Shapes Using a Non-Oscillatory Two-Component Pressure Approach (TPA). Modelling 2021, 2, 467–481. [Google Scholar] [CrossRef]
  39. Toro, E.F. Shock-Capturing Methods for Free-Surface Shallow Flows; John Wiley: Chichester, NY, USA, 2001. [Google Scholar]
  40. LeVeque, R. Finite Volume Methods for Hyperbolic Problems; Cambridge Press: Cambridge, UK, 2002. [Google Scholar]
  41. Malekpour, A.; Karney, B. Spurious Numerical Oscillations in the Preissmann Slot Method: Origin and Suppression. J. Hydraul. Eng. ASCE 2015, 142, 04015060. [Google Scholar] [CrossRef] [Green Version]
  42. ANSYS, Inc. ANSYS Fluent Theory Guide, Release 15.0; ANSYS, Inc.: Canonsburg, PA, USA, 2013. [Google Scholar]
  43. Ghidaoui, M.; Karney, B.W. Equivalent differential equations in fixed-grid characteristics method. J. Hydraul. Eng. 1994, 120, 1159–1175. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Schematic of computational cells in the method of characteristics and DGCM.
Figure 1. Schematic of computational cells in the method of characteristics and DGCM.
Mathematics 10 01960 g001
Figure 2. Data exchange between 1D and 3D computational zones.
Figure 2. Data exchange between 1D and 3D computational zones.
Mathematics 10 01960 g002
Figure 3. Iteration loop in the 1D-CFD coupling.
Figure 3. Iteration loop in the 1D-CFD coupling.
Mathematics 10 01960 g003
Figure 4. The schematic of the experimental apparatus. Figure adapted from Figure 2 in Ref. [9].
Figure 4. The schematic of the experimental apparatus. Figure adapted from Figure 2 in Ref. [9].
Mathematics 10 01960 g004
Figure 5. Pressure time histories calculated by DGCM at; (a) valve; (b) pipe mid-point.
Figure 5. Pressure time histories calculated by DGCM at; (a) valve; (b) pipe mid-point.
Mathematics 10 01960 g005
Figure 6. Pressure time histories calculated by 1D-2D at; (a) valve; (b) pipe mid-point.
Figure 6. Pressure time histories calculated by 1D-2D at; (a) valve; (b) pipe mid-point.
Mathematics 10 01960 g006
Figure 7. Cavity size evolution captured by 1D-2D.
Figure 7. Cavity size evolution captured by 1D-2D.
Mathematics 10 01960 g007
Figure 8. Pressure time histories calculated by 1D-3D at; (a) valve; (b) pipe mid-point.
Figure 8. Pressure time histories calculated by 1D-3D at; (a) valve; (b) pipe mid-point.
Mathematics 10 01960 g008
Figure 9. Cavity size evolution captured by 1D-3D.
Figure 9. Cavity size evolution captured by 1D-3D.
Mathematics 10 01960 g009
Figure 10. Cavity size evolution captured by CFD (Warda et al. [23]).
Figure 10. Cavity size evolution captured by CFD (Warda et al. [23]).
Mathematics 10 01960 g010
Figure 11. Pressure time histories calculated by MTPA at; (a) valve; (b) pipe mid-point.
Figure 11. Pressure time histories calculated by MTPA at; (a) valve; (b) pipe mid-point.
Mathematics 10 01960 g011
Figure 12. Cavity size evolution captured by MTPA.
Figure 12. Cavity size evolution captured by MTPA.
Mathematics 10 01960 g012
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Khani, D.; Lim, Y.H.; Malekpour, A. Calculating Column Separation in Liquid Pipelines Using a 1D-CFD Coupled Model. Mathematics 2022, 10, 1960. https://doi.org/10.3390/math10121960

AMA Style

Khani D, Lim YH, Malekpour A. Calculating Column Separation in Liquid Pipelines Using a 1D-CFD Coupled Model. Mathematics. 2022; 10(12):1960. https://doi.org/10.3390/math10121960

Chicago/Turabian Style

Khani, David, Yeo Howe Lim, and Ahmad Malekpour. 2022. "Calculating Column Separation in Liquid Pipelines Using a 1D-CFD Coupled Model" Mathematics 10, no. 12: 1960. https://doi.org/10.3390/math10121960

APA Style

Khani, D., Lim, Y. H., & Malekpour, A. (2022). Calculating Column Separation in Liquid Pipelines Using a 1D-CFD Coupled Model. Mathematics, 10(12), 1960. https://doi.org/10.3390/math10121960

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