Three-dimensional CFD analysis of liquid slug acceleration and impact in a voided pipeline with end orifice

ABSTRACT This research describes the dynamic behavior of an isolated slug driven by pressurized air in a voided line with an end orifice. A three-dimensional (3D) computational-fluid-dynamics (CFD) model is used to simulate the rapid propulsion and impact at the orifice for given slug length and driving air pressure, and is validated against experimental data. New mechanisms are observed: (i) the driving air pressure at the slug tail decreases with the slug motion; (ii) when the slug arrives at the orifice, the air fraction is almost one hundred percent in the most upstream part of the pipe, then it attenuates rapidly until an invariant eighty percent is achieved with a constant mass shedding rate; (iii) the velocity distribution in the radial direction of the cross-section at the midpoint of the slug length evolves from uniform to trapezoidal and then to logarithmic during slug movement; (iv) the initial vertical slug front changes its shape due to air intrusion at the top of the slug front; (v) the slug’s acceleration decreases first and then increases under the combined effects of its decreasing mass, nonlinear attenuation of the driving pressure, and increasing skin friction; (vi) the slug length has a constant rate of decrease.


Introduction
Two-phase gas-liquid flow exists in various piping systems and its different flow regimes have been classified, among which the slug flow has been extensively investigated because of its violent behavior and the complexities encountered in various applications and structures (Mohmmed et al., 2021). Much of the work published on slug flow is related to power plants, oil recovery and chemical industries. Although various achievements have been obtained for this particular flow regime, there are still many unknowns (Mohmmed et al., 2021), among which the recent attentions are flow regime transition (Deendarlianto et al., 2019;Kong et al., 2018), liquid film thickness (Liu et al., 2022;Liu et al., 2021), slug dynamic behaviors (Li et al., 2020;Parsi et al., 2017;Schmelter et al., 2021;Wang et al., 2018), effect of liquid viscosity (Bendiksen et al., 2018;Pineda-Perez et al., 2018), and its CFD modeling (Guerrero et al., 2017;Lopez et al., 2016;Pineda-Perez et al., 2018). Parsi et al. (2017) employed different qualitative and quantitative methods to compare the huge waves of vertical upward churn and churn-annular transition flow to pseudo-slug structures in horizontal and inclined pseudo-slug flows, and considered the huge waves and pseudo-slugs as the same liquid structures with different inception mechanisms. Extensive experimental investigations were conducted by Kong et al. (2018) in small and large pipe diameters to investigate the air-water slug flow's interfacial behavior utilizing four-sensor conductive probes. Their results indicated that the pipe diameter has immense effects on the interfacial structure of the slug flow. At the same time, the bubble size increased with the increasing pipe diameter. On the other hand, the constant electric current method was utilized by Deendarlianto et al. (2018) to investigate the stratified to slug flow transition in a horizontal pipe. Bendiksen et al. (2018) explained that the slug translational velocity is independent of slug length, while the slug front velocity is relatively dependent on slug length and it might vary significantly along the pipe. For liquid film thickness, when the Reynolds number of slug bubble is less than 3000, it increased with increasing capillary number, while the Reynolds number is over 3000, the liquid film thickness is approximately a constant (Liu et al., 2022;Liu et al., 2021). Although there are several attempts to analytically scrutinize the impact of the slug characteristics on solid structures (Li et al., 2020;Wang et al., 2018), obtaining a reliable solution has not been achieved.
Except the aforementioned traditional steady slug flow, there exist a more violent scenario that due to suddenly opened valve the liquid slug accumulated behind the valve is forced by the high gas pressure to move along the voided pipe like a 'bullet' in a gun. The slug velocity can easily be over 30 m/s. Liquid slugs moving in voided lines are often observed in piping systems that carry highpressure steam such as thermal and nuclear power plant piping. Due to its dangerous and troublesome nature, the phenomenon has received much attention in the power industry (Bozkus & Wiggert, 1997;Fabre & Line, 1992). The formation of liquid slugs in voided lines can be due to different mechanisms and the collection of steam condensate water at lower segments of the pipeline is one of them (Bozkus et al., 2004). When the system is reactivated, under high steam pressure, initially stationary liquid slugs formed at low sections in a voided line will be forced to move along the pipe and possibly accelerate to high velocities. When the fast-traveling slug hits obstructions such as bends, orifice plates and tees, severe damage to the piping system may be caused by the resulting impact force. Therefore, liquid slugs traveling in voided pipelines represent a high potential risk that needs to be assessed for plant security, design of new facilities and revamping of existing ones.
For slug motion in a voided line and impact on an obstacle, previous experimental and mathematical models have been reviewed by Hou, Tijsseling, and Bozkus (2014). With regard to the mathematical models, they are mostly one-dimensional (1D) (Bozkus et al., 2004;Bozkus & Wiggert, 1997;Hou, Tijsseling, & Bozkus, 2014;Kayhan & Bozkus, 2011;Owen & Hussein, 1994;Tay & Thorpe, 2014;Tijsseling, Hou, and Bozkus, 2016). By carefully incorporating extra 'resistance' due to flow separation at the bend into the 1D models, the agreement between predicted pressures and measurements has been improved (Hou, Tijsseling, & Bozkus, 2014;Tijsseling, Hou, & Bozkus, 2016). However, in reality, the motion of an isolated slug in a pipeline and its passing through an obstacle is a 3D phenomenon, which cannot be well captured by a 1D model, in particular the deformation of the water-air interface, water slug aeration, velocity distribution inside the slug, and the flow contraction at an obstacle. For slug impact at a target, apart from the geometry of the target, the most important factor is the impinging slug itself (its density, velocity, length and front shape) and the driving pressure behind it. To accurately obtain this information, a two-or three-dimensional model is needed.
The first attempt was by Yang and Wiggert (1998) who developed a quasi-two-dimensional (2D) model which allowed air entrainment. The liquid slug was treated as a number of concentric cylinders sliding through each other. Due to wall shear effects, the outer cylinders moved slower than the inner ones, which mimicked the gas penetration into the slug at both its front and tail. Since the effect of gravity was neglected, the computed flow remained axisymmetric. The comparison between the numerical results and the experiments of Bozkus and Wiggert (1997) showed that this 2D model largely overestimated the measured impact pressures.
With the slug variables calculated by the analytical model derived by Bozkus and Wiggert (1997), such as velocity, length, and driving pressure, just prior to slug impact on the elbow, Hou, Tijsseling, and Bozkus (2014) numerically investigated the slug flow passing a 2D elbow by using the Smoothed Particle Hydrodynamics (SPH) method. That is, the slug motion in the pipe was predicted by the 1D model (Bozkus & Wiggert, 1997), and upon arrival at the elbow, the more advanced SPH method was utilized. As another novelty, by incorporating the mechanism of flow separation at the bend as captured by SPH (Hou, Kruisbrink, Pearce, et al., 2014), a new 1D slug impact model was developed. Hou, Tijsseling, and Bozkus (2014) validated the models against the experimental studies of Bozkus and Wiggert (1997) and found that the calculated and the measured results were in good agreement. Recently, based on the same procedure, Dincer et al. (2018) numerically simulated the slug impact process at the elbow by a new 2D computer SPH code. The slug variables prior to impact were calculated by the improved 1D model of Tijsseling, Hou, and Bozkus (2016) and the model was validated against the experiments of Bozkus et al. (2004) in a 10 cm diameter pipe, which was twice the diameter of the pipe in the previous study (Bozkus & Wiggert, 1997). Another SPH simulation of an accelerating liquid slug was reported by Korzilius et al. (2017). Different from the weakly compressible SPH model used in Hou, Tijsseling, and Bozkus (2014) and Dincer et al. (2018), the incompressible SPH was adopted for straightforward enforcement of the pressure-gradient condition and for obtaining solutions with smoother pressure fields. Furthermore, a novel hybrid model that enabled one to leave out the larger inner part of the slug was also developed, which saved a significant amount of computational time. The obtained pressure forces exerted by the slug on a bend were to a large extent consistent with experiments (Bozkus et al., 2004), and the predicted impact speeds agreed well with analytical data (Tijsseling, Hou, & Bozkus, 2016). Note that with different pipe fittings at the end of the pipe, the impact mechanisms will be different, but the movement of water slug is more or less the same, provided that the resistance and compression of the air ahead of the slug is negligible.
As a mature technology, Computational Fluid Dynamics (CFD) has become a robust and reliable engineering tool. In recent years, it has been applied for solving a wide range of unsteady hydraulic engineering problems, even though huge computational cost is necessary. The need for large computer resources is still a major limiting factor in the daily engineering practice of 3D timedependent numerical simulations (Martins et al., 2014). In particular, when water-air interface motion occurs, it is necessary to implement multiphase modeling, for which advances in CFD have provided the basis. The volume of fluid (VOF) model was developed by Hirt and Nichols (1981) for simulating two or more immiscible fluids with an identifiable interface between them. VOF method shares a single set of equations for both gas-liquid phases. It uses ancillary transport equations based on the volume fraction to distinguish between the phases at each cell. Compared with other multiphase flow models VOF model can effectively predict and identifying the flow pattern (Guerrero et al., 2017). VOF method can track the interface and the liquid holdup, and hence is capable of predicting the flow regime (Pineda-Perez et al., 2018). This model has been widely used for studying two-phase pressurized flows, such as the free surface-pressurized flow in a ceiling-sloping tailrace tunnel (Cheng et al., 2007), transient cavitating flow (Passandideh-Fard & Roohi, 2008), transient flows with entrapped air pockets (Martins et al., 2017;Zhou et al., 2011), and geyser formation in vertical shafts (Chan et al., 2018). Recently, Schmelter et al. (2021) used the VOF model to investigate the influence of the interface instability imposed by the inlet velocity perturbation, and analyzed the effect of inlet perturbations to slug growth and initiation in the pipe.
Compared with theoretical analyses by simple models and experiments at laboratory scale with limited instrumentation, numerical simulation with CFD can be an efficient alternative to gain detailed insight into flow and pressure fields, something which is highly needed for the present problem of slug motion in a voided line. Therefore, the main objective of the current work is to explore the use of a 3D CFD model describing unsteady twophase flow for a better understanding of the slug dynamics during its motion in a voided line with an end orifice. Published data from an experimental apparatus (Owen & Hussein, 1994) are used for comparison with the CFD results.
The rest of the paper is organized as follows. Section 2 presents the employed 3D mathematical models for fluids and moving water-air interfaces. The solution strategy and the employed numerical techniques are also described. To validate the CFD model, physical experiments available in the literature are used, of which the setups are detailed in Section 3. The obtained results are presented and discussed in Section 4 through comparison with laboratory experiments and previous theoretical models. Several key factors determining the validity of the assumptions in the 1D models are checked and analyzed. Section 5 lists the concluding remarks.

Navier-Stokes equations
The governing equations present the velocity and pressure characteristics of fluids by the conservation equations of mass and momentum. For incompressible fluid flow they are where t = time, x i , x j = spatial coordinate in i-and jdirection (i, j = 1, 2, 3), respectively, ν = kinematic viscosity, ρ = fluid density, u i , u j = velocity components, p = pressure, and g i = acceleration due to mass force in i-direction. Note that for the water air two-phase flow studied here, they share the same set of continuity and momentum equations.

Turbulence model
Due to the large scale of the simulation zone and the high Reynolds number, it is not possible yet to solve the 3D Navier-Stokes equations using direct numerical simulation. A common way is to use a turbulence model to close the system. Taking the standard turbulence model as an example, in the water-air two-phase flow, the turbulence model of mixed phase can be solved directly or the turbulence model of each phase can be solved separately.
In this study, a single-phase turbulence model is used to deal with the turbulence in multiphase flow. The standard two-equation k-model, which includes the transport k-equation and -equation, is given by where k = turbulence kinetic energy, = turbulence dissipation rate, µ = dynamic viscosity [negligible compared to eddy viscosity], µ t = turbulence eddy viscosity, σ k and σ ε = turbulence Prandtl numbers of k and , respectively, with σ k = 1.0 and σ ε = 1.3, constants C 1 = 1.44 and C 2 = 1.92, G k = turbulence kinetic energy generated by the average velocity gradient. Note that these default values are determined from a large number of physical experiments, and have been widely used in hydraulic transients in pipes with reliable results. The constants are the same for both water and air. The expression for µ t is where C μ = 0.09 is an empirical constant, and the expression for G k is Note that due to the extremely short duration of slug impact, turbulence and gravity may not play an important role in the impact process as assumed and validated in Hou, Tijsseling, and Bozkus (2014) and Dincer et al. (2018); however, they highly affect the acceleration of the isolated slug in the pipe.

VOF model
To track air-water interfaces, the VOF model based on multiphase flow theory is employed. It is an efficient numerical technique to track interfaces between two or more immiscible fluids. In the calculation, all fluids share the same turbulence model with the same values of the empirical constants. The interface is tracked by introducing the volume fraction of the i-th phase in each cell, α i , which satisfies where n = the number of fluids (gas and/or liquid). For the air-water flow considered herein, we have α w + α a = 1. When α w = 1 or α a = 1, a cell contains only water or air, otherwise it contains the air/water interface. The evolution of volume fraction for the water phase satisfies where u i is the velocity of both water and air, i.e. there is no slip at the interface of the two phases.
The fluid properties, such as ρ and µ, in each cell are adjusted according to the volume fractions, e.g. the mixed fluid density in each cell is

Boundary and initial conditions
The upstream boundary herein is a reservoir or tank that is modeled as an inlet with imposed pressure. The downstream boundary is an outlet consisting of an orifice open to atmosphere. With the change of wall humidity, the real wall boundary condition might be between the no-slip boundary and free slip boundary. However, as the wall is assumed to be fixed and dry here, the no-slip boundaries are selected for the orifice and the entire pipe wall. The initial liquid slug at rest (see Figure 1) is under given tank pressure separated by a sheet or a valve from stationary air under atmospheric pressure. The simulations with the models introduced above are carried out by Fluent in ANSYS fluent (version 18.0). Owen and Hussein (1994) for liquid slug motion and impact on orifice plate.

Test problem and numerical method
With given orifice geometry and pipe geometry, the slug velocity just before impact is the most important factor in the validation of the numerical model of liquid slug motion in a voided line with an end obstacle, because it determines the magnitude of the impact force (Hou, Tijsseling, & Bozkus, 2014). To our best knowledge, this velocity was measured only in the work of Owen and Hussein (1994), which therefore is used for validation herein.

Laboratory setup description
The test rig of Owen and Hussein (1994) is shown in Figure 1. The stainless-steel pipe length L p from valve to orifice was 13 m and the inside pipe diameter D was 50 mm. The downstream end of the pipe was outfitted with a circular orifice with a diameter of 25 mm. Air was used as the driving gas. As the volume of the air reservoir was large, the drop in pressure was less than 3% when the air was expelled. For the experiments, a quarter-turn butterfly valve was used between the air reservoir and the water-filled pipe to enable a rapid start of the event. The water was held between the closed butterfly valve and a thin polythene sheet, sandwiched between flanges, downstream of the valve. After sudden valve opening causing sheet rupture, the water slug was accelerated by the compressed air. To measure the arrival velocity of the slug at the orifice plate, two conductivity probes, 0.215 m apart, were used to detect the leading edge of the slug as it passed the probes whereby the time interval between the slug passing the two positions was recorded. To measure the impact at the orifice plate, a piezo-electric pressure transducer was located close to its upstream face at the top of the pipe (0.01 m in front of the orifice, see Figure 1). Since the end of the pipe was open to the atmosphere, the initial air pressure in the pipe was atmospheric. The initial air temperature in the reservoir was not measured but close to ambient. Due to the large tank pressure, one experimental run was completed within one second.

CFD model definition
The 3D geometry definition consists of specifying the shape of the physical boundaries of the two contained fluids air and water. After the geometry definition, CFD simulation requires the division of the fluid domain into a large number of smaller, non-overlapping subdomains resulting in a mesh of cells. The continuous equations that describe the unknown flow variables are numerically discretized in each of the defined cells. The final solution accuracy is strongly dependent on the number of cells in the mesh within the computational domain (Wang, 2004) and on the numerical time step. The calculation area of the CFD model is a circular pipe of 13 m in length and 0.05 m in diameter and an additional pipe segment with length 2D. The two pipe segments are connected by a plate with a centered orifice of 0.025 m in diameter and of 0.2D in thickness (see Figure 1). Instantaneous valve opening is assumed. The initialization of the fluid phases (liquid and gas) is carried out in the VOF model by defining the percentage of each phase contained in each cell.
The standard wall function method is used in this research. As the material of the pipe is stainless steel, according to Lv et al. (2002), the wall roughness height is in between 0.1 and 0.2 mm and it is taken as 0.1 mm herein. The inlet pressure (gauge) is set according to the test condition without the reported maximum 3% drop (Owen & Hussein, 1994).
In the 3D CFD model of this study, the finite-volume method (FVM) is employed to discretize the governing equations. The robust Pressure Implicit with Splitting of Operators (PISO) procedure is employed to accomplish the coupling of velocity field and pressure distribution. To precisely determine the location of the water-air interface, the geometric reconstruction method is employed. The numerical convergence can be assessed by progressively monitoring the global imbalances in conservation of mass and momentum through each iteration step. A converged solution, per time step, is achieved when the imbalances fall below a pre-set tolerance (10 −6 herein) (Martins et al., 2016;Martins et al., 2017).
Nine cases with different tank pressures (P 0 = 240. 0,309.3,386.4,452.0,482.8,544.3,552.0,586.8,and 656.2 kPa,gauge) at the downstream position of the valve (the tank with compressed air and valve resistance are not modeled) and a constant initial slug length (L s0 = 3.15 m) have been studied in this work.

3D-VOF mesh independence analysis
A larger number of cells generally leads to a more accurate solution, at the expense of computational resources. A systematic study has been carried out for obtaining the most efficient mesh, which meets two contradictory criteria: maximum accuracy and minimum computational effort. According to Wang (2004), to get reliable results, the ratio of the effective grid size in different directions should not be larger than two. Nine mesh sizes (uniform quadrilateral mesh) of the cross section of the pipeline are taken from 2 mm to 6 mm, resulting in a number of nodes between 148,500 and 1,340,000. Grid refinement is applied near the orifice. As the implicit algorithm was applied for time marching, numerical stability can be ensured even for a large flexible time step. However, For the case with a tank pressure of 240 kPa, the arrival velocity of the water slug at the orifice plate is shown in Figure 2. It is clear that the number of nodes has a large effect on the results, but that convergence can be achieved. As the grid number increases, the arrival velocity tends to increase. When the grid number is 900,000 or larger, the computed arrival velocity do not change much. Considering the accuracy versus computational time, the number of nodes used in all the simulations is 862,000. For one case, it takes roughly 48 h with 48 CPU cores.

Arrival velocity of the slug
The calculated arrival velocities of the slug at x = 12.99 m (just upstream of the orifice) are shown in Figure 3 in its dependence on the driving pressure, together with the measured velocities. To get the best agreement between the calculated and measured velocities for nine cases, the measuring conductivity probe is assumed to be at a vertical distance of 1.9 mm to the top of the pipe wall. For such a high-speed water slug, the Reynolds number is more than 10 million. The thickness of the corresponding boundary layer is so small that the measuring position was hardly in its vicinity. Some uncertainties that can lead to differences include tank pressure, pipe wall roughness height, shape of the leading water-air interface, cross-sectional velocity distribution, and measurement error. As noticed by Owen and Hussein (1994), the tank pressure after air release was not constant but dropped by a maximum of 3%. The used roughness height of 0.1 mm is the lower bound for stainless steel pipes (Lv et al., 2002). As examined below, the roughness height most likely affects the slug dynamics. Affected by gravity and air intrusion, the shape of the leading water-air interface can hardly remain flat. Due to wall shear stress and gravity, the velocity distribution in the cross section close to the leading water-air interface is extremely complex. The last two factors result in high sensitivity of the measured velocity to the length and circumferential location of the conductivity probe (see Figure 3).

Impact pressure
To measure the impact pressure at the orifice plate, a piezo-electric pressure transducer was located close to its upstream face at the top of the pipe (0.01 m in front of the orifice; see the reference point in Figure 1). Figure 4 shows the calculated pressure histories at the location of the transducer. An impulse-like response is observed for all the cases, and before the slug passes the measurement location, the pressure is close to the atmospheric pressure and hence equals to zero. In addition, the larger the tank pressure, the earlier the slug arrival at the orifice plate, and the larger the impact's peak pressure. These conclusions are according to expectations.
As shown in Figure 5, both the measured and calculated maximum impact pressure (gauge) increase linearly with the tank pressure, but the calculated peak pressures are smaller (54%) than those measured. The large deviation by a factor of two might be attributed to liquid elasticity (water-hammer effect), which − as experimentally observed (De Martino et al., 2008;Zhou et al., 2002) − may play a vital role in slug impact at a dead  end with a central orifice. While in the current CFD simulations, compressibility has not been considered neither in liquid nor gas. In fact, we had considered the air compressibility, and found that it had little effect on the movement of water slug due to the relatively large orifice size, but significantly slows down the computation. To save computational time, the current CFD simulations neglected the air compressibility. A second − probably more plausible − reason is that small amounts of air remain trapped and compressed at the upstream top corner of the orifice (Martin, 1996;Martin & Lee, 2000). Nevertheless, some useful results to explain the measured impact pressures are given next.
The impact pressures are displayed in Figure 6 as an amplification of the tank pressure. The amplification factor is more or less constant, no matter what the tank pressure is, which is consistent with the experimental data of Owen and Hussein (1994) who proposed the following relation for the force F on the orifice plate: where A and A 0 are the areas of the pipe and the orifice, respectively, U max is the slug velocity with which it hits the orifice plate, and L s0 is the initial slug length. The velocity U max depends on tank pressure, slug length, pipe length, slope and diameter, friction factor, amount of holdup, orifice size, etc. Note that the original correlation in Owen and Hussein (1994) is (A 0 /A) 2 , which is physically not correct, because for a given velocity, the smaller the orifice, the larger the impact force should be. The limit cases are a closed end (A 0 = 0) and an open end (A 0 = A). Equation (10) gives an infinitely large force for a closed end, which is not correct because gas compressibility and liquid inertia lead to an oscillating but finite pressure (Martin, 1976). When the end is open (A , p. 00 = A), the resistance of the orifice plate to the moving liquid slug can be neglected, and Eq. (10) is valid only for 'moderate' values of A 0 . It can be seen from Figure 7 that the present numerical data correlate well with Eq. (10) when the constant takes a value of about 4.6. To verify the CFD-VOF model for impact pressure calculation, a water jet impinging on a vertical wall is simulated and the results are shown in Appendix A.

Other flow characteristics
Based on the numerical simulations, some relevant hydraulic characteristics of the slug's motion in the voided line are studied, especially those that have hardly been measured in the experiments published so far, but that are of importance for the development of accurate 1D models. They are the changes of the slug itself (its velocity, length, shape, etc), the pressures at the slug tail and front. Hereafter, the slug front and slug tail are defined as the most upstream point of the leading and secondary water-air interfaces, respectively. The slug length is then the distance between the slug tail and the slug front; see the sketch shown in Figure 8.

Pressures at the slug tail and front
In the reported 1D models, the driving air pressure at the slug's tail is mostly assumed to be the tank pressure, although Bozkus and Wiggert (1997) modeled pressure waves in the driving air column. That is, the pressure loss (due to friction) and pressure gradient (due to acceleration) in the lengthening air column are neglected. If there is no liquid mass shedding from the slug, this assumption might be reasonable because the inertia of the air and the skin friction can be neglected relative to the large tank air pressure. However, as observed in experiments for pipe emptying (Laanearu et al., 2012), mass shedding occurs at the tail of the liquid column, due to which the driving air pressure at the slug's tail is different from the tank pressure.
The nine simulated driving pressures stay not equal to the prescribed tank pressure: after a constant period, they decrease nonlinearly in time (Figure 9 left) when the slug moves ( Figure 10 left) until its arrival at the downstream orifice, where S is the displacement of the slug's tail. For the studied nine cases, the slug arrives at the orifice plate at t max = 0.36 s (highest tank pressure) to t max = 0.61 s (lowest tank pressure). The normalized driving pressure v.s. the normalized moving time and the displacement of the slug's tail are shown in Figures 9 and 10 right, respectively. It is interestingly noted that the normalized driving pressure v.s. the normalized moving time follows the same quadratic relationship (Figure 9 right). So does the normalized driving pressures v.s. the displacement of the slug's tail (Figure 10 right). Although the larger the tank pressure is, the more the driving pressure at the slug tail reduces (see Figures 9 left and 10 left), the relative pressure reduction to the tank pressure is close to each other (see Figures 9 and 10 right). For example, for the case with a tank pressure of 240 kPa, the pressure decrease is 85 kPa (35% pressure reduction), whereas for the case of 656.2 kPa, it is about 260 kPa (40% pressure reduction). They imply that for a liquid slug with given length, the ratio of the driving pressure to the tank pressure evolves in a similar pattern no matter what the tank pressure is. Since the air density is so small that the pressure decease due to kinetic energy increase of the driving air is negligible. If the gas dynamics in the lengthening air column is further neglected, the large pressure decrease can only be attributed to the fact that part of the air pressure does work on the shed mass. In other words, the shed mass gives extra resistance to the blowing air. The same dimensionless driving pressure decrease pattern in time (Figure 9 right) has the same origin as that in distance (Figure 10 right).
The changes of the air pressure ahead of the slug are shown in Figure 11. In contrast to the pressure at the slug's tail, the air pressure at the slug's front increases nonlinearly in time (Figure 11 left) when the slug moves ( Figure 11 right) until its arrival at the downstream orifice. The normalized air pressure at the slug's front v.s. the normalized moving time and the displacement of the   slug's tail are shown in Figures 11 and 12 right, respectively. Similarly, it is noted that the normalized air pressure at the slug's front v.s. the normalized movement time follows the same pattern (Figure 11 right). So does the normalized air pressure v.s. the displacement of the slug's tail (Figure 12 right). Although the larger the tank pressure is, the more the air pressure at the slug's front increases (see Figures 11 left and 12 left), but the relative pressure increase to the tank pressure is close to each other (see Figures 11 and 12 right). For example, for the case with a tank pressure of 240 kPa, the pressure increase is 12.5 kPa (5.2% pressure increase), whereas for the case of 656.2 kPa, it is about 35 kPa (5.3% pressure increase). This implies that for a liquid slug with given length, the ratio of the air pressure at the slug's front to the tank pressure evolves in a similar pattern no matter what the tank pressure is.
Compared to the tank pressure, the air pressure buildup ahead of the slug is small. The air pressure at the slug front as function of the square of the cross-sectional average of the slug front velocity is shown in Figure 13 where a linear relationship is observed.
When the slug front arrives at the orifice, the ratio of the pressure at the slug's tail to the tank pressure is a constant for the studied nine cases. Similarly, the ratio of the air pressure buildup at the slug's front to the tank pressure is a constant too as shown in Figure 14. They are consistent with the observation shown in Figures 9 right, 10 right, 11 and 12 right.
At the time the slug front arrives at the orifice, the cross-sectionally averaged pressures normalized by the tank pressure along the pipeline are shown in Figure 15. The slug's tail is near S = 12 m, so that the large, nearly constant, negative pressure gradient between S = 12 m and S = 13 m is within the water slug. The upstream air pressure is nearly constant up to S = 3 m (largest tank pressure) and S = 6 m (lowest tank pressure). The decreasing air pressure in the middle is exactly where holdup (shed mass) is present. Apparently, the resistance of shed mass to the moving air above the water-air   interface is the main cause for decreasing pressure. The consequence is that the driving air pressure at the slug's tail is mostly assumed to be the same as the tank pressure in previous studies, which is clearly not the case. That is to say, to develop an accurate model for the concerned isolated slug problem, the decreasing driving-air pressure due to the interfacial drag between the low-velocity holdup and high-velocity air must be taken into account.

Flow pattern, holdup and air fraction
When the slug arrives at the orifice, the water-air pattern in the pipeline for nine cases is shown in Figure 16 together with the initial state. A similar holdup profile ('trapezoidal' distribution) is formed behind the moving slug, no matter what the tank pressure is. As explained above, the shed mass in the early stage of the slug motion  is blown downstream along the pipe and hence there is almost no water in the most upstream part (downstream of the valve, about 1 m). This is consistent with the experimental observation for pipe emptying with large air pressure (Laanearu et al., 2012). In addition, the shed mass accumulates at the bottom of the middle part of the pipe due to gravity and air blowing from behind and on top of it, resulting in increased height. A nose shape of the air behind the slug tail is seen in Figure 16, which is consistent with gravity current measurements on long air (Taylor) bubbles (Shosho & Ryan, 2001).
The length of the region of constant holdup cannot be determined from Figure 16, because it shows the central plane only. The thickness of the shed mass seems decreasing from 6 m to 11 m, but in fact, the holdup is more or less constant because part of it is distributed around the pipe cross-section. The shed mass needs time to accumulate at the pipe bottom (due to gravity). Figure 17 gives the proper picture, which shows the cross-sectional water fraction along the pipeline when the slug arrives at the orifice plate. In the part of the pipe downstream of the valve (about 1 m), the water fraction is almost zero. It increases from there and has a maximum at about 5 m. Then a more or less constant (about 0.25) level is maintained for all the considered cases. Therefore, a constant mass shedding rate may be assumed as done in previous 1D models.
When the slug arrives at the orifice, the front's water fraction is less than 0.9 (Figure 17), which might imply that there is no complete slug left because air penetration occurred. For an air fraction of 0.1, the corresponding thickness of the intruded air layer is 1.25 mm if pure air is assumed. However, due to complex water-air interactions in the wall boundary region, there is not a top layer with pure air, but a water-air mixture. An effective method to judge the slug flow and pseudo-slugs (Parsi et al., 2017) is to check whether there is a large pressure gradient between slug's front and slug's tail. For the studied isolated slug flow, the pressure at the slug's front is low, and it is high at the slug's tail, resulting in a large pressure gradient. This is not the case for pseudo-slugs, since in this instance the air can penetrate through the water slug (Parsi et al., 2017). During the entire movement of the water slug in the pipeline, by comparing Figure 9 and Figure 11, we can find that the driving air pressure scaled by tank pressure is from 1 to 0.6, while for the pressure scaled at the slug's front it is from 0 to 0.06. The minimum difference is more than a factor 10. Therefore, the intruded air is always separated by the liquid slug from the driving air at its tail, and the high driving pressure behind the slug cannot release via the intruded air from the slug's front. Figure 18 shows the water fraction along the pipeline at different times for the case with tank pressure of 656.2 kPa. The mass shedding starts at the early stage due to wall shear stress and gravity. Under the effect of Kelvin-Helmholtz instability, however, the early shed mass cannot stay stagnant, but is blown downstream along the pipe, resulting in an empty pipe section of increasing length (Figure 16). It is about 1 m for all the cases (Figure 17) when the slug arrives at the orifice plate.

Slug length variation during acceleration
The variation of the dimensionless slug length with time is shown in Figure 19. When the slug travels in the pipe, its length is nearly constant initially, but then suffers from Figure 16. The water slug in the pipe when it nearly hits the orifice plate, for different tank pressures; water is blue, air is red; vertical central plane is displayed. Figure 17. Water fraction distribution along the pipeline when the slug arrives at the orifice for different tank pressures. a rapid decrease due to mass shedding, the final result of which is seen in Figure 17. The final length of the slug is the same (about 1.06 m) for all tank pressures as shown in Figure 20. In other words, the total amount of holdup is independent of tank pressure, which is useful information for estimating the average mass shedding rate. It can be attributed to that the driving pressure is so large that the shedding mass is solely affected by the displacement of the slug.
To analyze the mass lost per traveled unit length, the variation of the slug length as a function of S is shown in Figure 21, where the decreasing pattern of slug length during traveling is seen, which is basically the same for all nine driving pressures. The constant mass-loss rate is of about 0.176 m/m, which in terms of mass is 0.346 kg/m. This is exactly what is assumed in 1D models (Bozkus & Wiggert, 1997;Tijsseling, Hou, & Bozkus, 2016). The linear fit in Figure 29 implies that L s approaches zero for a longer pipe length S, i.e. the slug would disintegrate. With the slug motion in the pipeline, its length continuously shortens, until the water slug disperses under the action of aeration, gravity and mass shedding. As no water mass is gained at the slug front, a steady state with constant slug length cannot be reached.

Velocity distribution within the slug
The cross-sectional average of the velocity in the slug is usually assumed to be uniform (Bozkus et al., 2004;Bozkus & Wiggert, 1997;   2014) or linearly decreasing from tail to front (Tijsseling, Hou, & Bozkus, 2016) in 1D models. The uniform velocity assumption fulfills mass conservation, while linearly decreasing velocity is based on linearly decreasing 'internal' holdup.
It is difficult to physically measure the velocity distribution in a deforming liquid slug, and CFD is, therefore, a justified alternative to obtain reasonable solutions by adopting proper setup and mesh size, The velocity distribution within a cross-section in the middle of the slug length at different times is depicted in Figure 22. In the initial stage, the cross-sectional velocity distribution is like that of developing entrance flow, which evolves to a laminar-like velocity profile. The latter is well known for accelerating flows (Koppel & Ainola, 2006). The initial symmetry of the distribution is clear. Although the effect of gravity is seen from the velocity distribution in When the slug arrives at the orifice plate for the case with tank pressure of 240 kPa, its shapes in the horizontal central plane (top view) and in the vertical central plane (side view) are shown in Figure 23. The slug's front is at S ≈ 12.96 m and the slug's tail at S ≈ 11.90 m, and hence the slug length is about 1.05 m. In the horizontal plane, the symmetry of the shed mass is as expected due to the shear force at the pipe wall. This observation is consistent with the assumption of Yang and Wiggert (1998) in their 2D model, where gravity was neglected. However, in the vertical plane, due to gravity, the shapes of the slug at both its front and tail are not axisymmetric. The distributions of horizontal velocity in different crosssections of the slug are presented in Figure 24. In the horizontal plane (Figure 24 left), the velocity distribution is symmetric, although at different cross sections it varies significantly. The gravity effect on the velocity distribution in Figure 24 (right) is more significant than that in Figure 22 (right). The front becomes wedge-shaped. The heterogenous velocity profiles will cause deformation of the slug's shape. For example, the velocity distribution is no longer symmetrical along the longitudinal direction of the water slug, resulting in its asymmetric deformation. Along the central axis (y = 0), the velocity of the slug decreases from its tail to its front.
The cross-sectional average velocities (used in 1D models) shown in Figure 25 are about constant along the whole pipeline as required by mass conservation. And the larger the tank pressure, the higher the flow velocity. The sudden velocity increase is due to the reduced cross-sectional area at the orifice plate.

Front and tail deformation
In 1D models, the slug front and tail are usually approximated by vertical planes where no air entrainment occurs, with Tijsseling, Hou, Bozkus, et al. (2016) as an exception. However, in real situations air intrusion at the top of the moving slug front may occur as observed in the experiments of rapid pipe filling (Hou, Tijsseling, Laanearu, et al., 2014). For the case with tank pressure of 240 kPa, the shape changes of the slug front and tail are shown in Figure 26. To clearly show the complete shape of the water slug, the middle part of the water slug is omitted. In the initial stage (0-0.2 s), the slug front can still be assumed to be vertical, although deformation starts to develop. Air intrusion indeed occurs at the sides, top and bottom of the slug front, leading to a bullet-shaped front. This bullet shape seems to be induced by the driving (Taylor) bubble of compressed air at the tail. So, the unsteady flow regime is annular instead of stratified. The slug shapes in the horizontal and vertical central planes are different because of gravity, which gives more water at the tail's bottom and more air at the front's top (Figure 26 right). When the water slug with changed front hits the orifice plate, the air around the humped front can be easily trapped at the orifice plate's top.

Slug motion
The simulated cross-sectionally averaged velocities of the mid slug in dimensionless form (normalized with √ P 0 /ρ w ) is shown in Figure 27, in which the same pattern is seen for all different tank pressures.
The cross-sectional acceleration of the slug in dimensionless form with P 0 /(ρ w L s0 ) against the dimensionless time with T max (the total time when the water slug arrive at the orifice) are shown in Figure 28. Under the combined effects of the increasing friction, the decrease of the driving pressure and the reducing slug mass, the slug acceleration decreases first and then increases again. In the first half of the pipeline, the slug acceleration decreases due to decreasing driving pressure at the slug tail, increasing skin friction and pressure buildup at the slug front. In the second half of the pipeline, the length (and mass) of the water slug decreases rapidly due to mass shedding. The driving pressure is still decreasing   (see Figure 9 right), the pressure buildup is increasing (Figure 11 right), and the shear stress on unit mass is increasing due to enlarging slug velocity (see Figure 27). However, the acceleration is more affected by the rapid mass reduction of the slug, which makes the acceleration increase rapidly.
The variation of the normalized slug velocities with the position of the slug's tail for different tank pressures is shown in Figure 29. With the slug movement, the slug velocities rapidly increase at the early stage (within the first 2.5 m), and then increase almost linearly until the slug arrives at the orifice. The link between Figures 27, 28 and 29 is: dU dx = 1 U dU dt .

Effect of wall roughness
Since the roughness of the pipe wall determines the friction force, it has a rather big effect on the slug movement  in the voided line. The roughness is usually specified by the manufacturer and can be affected by many factors, such as surface humidity, pipe age, and state of erosion.
To study the influence of roughness height on the slug motion, friction force, which is proportional to the square of the slug velocity, plays an important role. Therefore, when the largest tank pressure of 656.2 kPa is selected, the largest slug velocity and corresponding friction are obtained. In fact, other cases were studied too, and similar trends were obtained. Therefore, six different roughness heights from 0 to 1 mm are selected for the case with a tank pressure of 656.2 kPa. When the slug arrives at the orifice, the cross-sectional average velocities and pressures along the pipe are shown in Figures 30 and 31, respectively. Again, the sudden velocity increase is due to the reduced cross-sectional area at the orifice plate. As expected, the velocity and pressure decrease with increasing wall roughness. When the roughness height changes from 0 to 0.2 mm, the velocities and pressures decrease significantly, but when the roughness height is from 0.2 to 1.0 mm, the differences are not so significant. Note that for a roughness height of zero (hydraulically smooth pipe), skin friction still exists (dependent only on the Reynolds number). When the water slug reaches the orifice, the influence of the roughness height on the slug length, the average slug velocity and the driving air pressure are shown in Figure 32. For a roughness height less than 0.2 mm, the slug length and slug velocity decrease sharply. When the roughness height is larger than 0.2 mm, the decrease rate is relatively small. The driving air pressure decreases rapidly in the range from 0 to 0.4 mm, and when the roughness height is 0.4 to1.0 mm, the driving pressure decreases slightly. For a large roughness, more air will be initially attached to the wall surface, which will affect the aeration of the slug at the wall surface. For stainless steel, the value of the wall roughness height is in the range of 0.1-0.2 mm (Lv et al., 2002), which means that the influence of the roughness on the slug's motion cannot be ignored.

Conclusions
A 3D CFD model has been applied to simulate a water slug propagating in a horizontal voided pipe with an end orifice. The calculated arrival velocities of the slug at the orifice are in good agreement with the measurements of Owen and Hussein (1994), with an average deviation of about 6% when the measuring conductivity probe is assumed at a vertical distance of 1.9 mm below the top of the pipe wall. However, the predicted impact pressures are about half (54%) of the measured ones, and this large discrepancy might be attributed to the waterhammer effect or the air entrapped in front of the orifice, which cannot be accurately calculated with the current CFD model. In both the measurements (Owen & Hussein, 1994) and the numerical simulations, the ratio of the impact pressure to the tank pressure is constant, at least for the studied cases.
The driving air pressure at the slug's tail is not equal to the tank pressure as it decreases with the slug's advance. For all studied cases, a 35% to 40% pressure decrease is observed. When the slug has arrived at the orifice, the distribution of the cross-sectionally averaged air pressure along the pipe shows three stages (from upstream to downstream): a small decrease, a large decrease and a moderate decrease. The air fraction is almost 100% in the most upstream part (about 1 meter) of the 13 m long pipe, and this value decreases to 80% to 70% the next 4 m before an approximately constant mass shedding rate is achieved. This implies that the shed mass is not uniformly distributed, and the velocity of it is not zero.
The velocity in the liquid slug along its centerline decreases from tail to front, and the difference in velocity between tail and front increases with time due to the air intrusion at the slug's front during its propagation. The distribution of horizontal velocity in the cross-section at the slug's midpoint evolves from uniform to 'trapezoidal' and then to 'logarithmic' during the slug's movement. The cross-sectional averaged velocities are uniform along the pipeline, which fulfills the conservation-of-volume law for incompressible fluids.
The slug acceleration decreases first and then increases under the combined effects of reduced slug mass, nonlinear decrease of driving pressure, and increasing nonlinear friction.
The slug tail takes the form of a Taylor bubble, whereas the slug front is about vertical in the initial stage and then changes its shape due to air entrainment mainly but not exclusively at the top of the slug's front. The shapes of the slug front and tail are symmetrical in the horizontal plane, but in the vertical plane, due to gravity, the water fraction at the bottom is slightly higher than at the top. The length of the slug when it arrives at the orifice is independent of the tank pressure for all the simulated cases, which means that the slug length decreases at a constant rate per unit length.
The above interpretations are useful for a better understanding of the features of slug motion in a voided line, and helpful for the development of accurate and efficient 1D models in the future.