Numerical Simulation of Elongated Bubbles and Liquid Films in Horizontal Slug Flow

1 Introduction

Subsea multiphase transmission pipelines, frequently utilized in offshore oil and gas production and Floating Liquefied Natural Gas systems, predominantly feature horizontal or near-horizontal configurations [1,2]. In the co-transport of oil and gas within these complex systems, slug flow is a prevalent flow pattern [3,4]. It is fundamentally characterized by dynamic fluctuations, spatial intermittence, and complex phase distributions [5,6]. These transient fluctuations induce cyclic loading, posing threats to the structural integrity of the pipeline and the operational stability of downstream separation facilities [7,8]. Recent studies have underscored the importance of accurately modeling pressure loss mechanisms and near-wall flow characteristics in such pipeline systems [9,10]. Furthermore, the continuous push for high-fidelity computational predictions has led to extensive numerical and experimental validations of transient two-phase flow dynamics, particularly focusing on the precise capture of intricate interfacial structures and phase distribution in horizontal pipelines [11,12]. Investigating irregular flow behaviors, such as wavy two-phase flows in manifolds, has also proven critical for optimizing horizontal well deployments [13]. Additionally, the rheological properties of the transported medium, such as the flow behavior of waxy crude oil modified by bionic nanocomposite pour point depressants, play a vital role in determining overall pipeline flow characteristics [14].

Historically, early investigations into slug flow mechanisms relied on 1D theoretical formulations and macroscopic time-averaged variables [15], which laid the groundwork but are insufficient for capturing local transient phenomena. The initiation, development, and hydrodynamics of liquid slugs have been widely debated, governed primarily by two major theoretical frameworks: the interface instability theory and the liquid slug stability theory. Taitel and Dukler [16] established a foundational model for predicting flow regime transitions based on the Kelvin-Helmholtz (K-H) instability mechanism, proposing that local pressure variations driven by Bernoulli forces cause interfacial waves to grow and eventually bridge the pipe cross-section. Building upon this, Lin and Hanratty [17] and Wu and Ishii [18] utilized linear stability theory to predict slug initiation, concluding that the instability of finite-amplitude interfacial waves acts as the direct and necessary precursor to slugging. Furthermore, Barnea and Taitel [19] refined the K-H stability criteria for stratified flow by comparing viscous and inviscid approaches, highlighting the critical role of fluid viscosity in phase pattern transitions. Experimental observations by Kordyban [20] on the transition to slug flow in the presence of large waves further demonstrated the significant interfacial interaction occurring directly over wave crests.

Conversely, the liquid slug stability theory emphasizes mass conservation and dynamic balance at the slug boundaries. Jepson [21] modeled the transition to slug flow by focusing on the effect of pipe diameter on liquid flux and the physical slugging process. Ruder et al. [22] and Dukler et al. [23] thoroughly investigated the necessary conditions for the existence of stable slugs, establishing physical models to predict the minimum stable slug length required to sustain the flow. Providing a quantitative mechanism, Woods and Hanratty [24] related slug stability directly to the shedding rate, positing the conclusion that a liquid slug can only sustain itself if the rate of liquid entrainment at the front precisely matches the liquid shedding rate at its tail. The exact mechanisms of gas entrainment, however, remain deeply tied to 3D interfacial shear that 1D models cannot resolve.

While these 1D analytical approaches laid a macroscopic foundation, they inherently fail to capture the complex 3D topological evolution of elongated bubbles. As demonstrated by Lu [25], through various computational studies, 1D two-fluid models exhibit significant discrepancies compared to experimental data regarding slug frequency and spatial development distance. To resolve these fine-scale phase interface structures, 3D Computational Fluid Dynamics (CFD) has become an essential analytical tool [26]. The application of CFD has been significantly broadened in recent years, ranging from the three-dimensional analysis of severe liquid slug acceleration and impact within voided pipelines [27] to the precise optical monitoring and numerical evaluation of slug-flow velocities inside microchannels [28]. Additionally, hybrid approaches combining mechanistic CFD simulations with machine learning algorithms are now being deployed to predict slug formation in complex oil-gas-water three-phase pipelines [29]. The Volume of Fluid (VOF) method [30], known for its strict mass conservation, has been widely adopted in this domain. Ekambara et al. [31] utilized the VOF approach to simulate bubbly two-phase flow in horizontal pipes, resolving the internal phase distribution and concluding that a maximum void fraction consistently appears near the upper pipe wall. Similarly, Andrianto et al. [32] performed comprehensive CFD studies on gas-liquid slug flow, quantitatively predicting bubble lengths and liquid holdups that aligned closely with experimental measurements. Schepper et al. [33] applied CFD modeling to predict various gas-liquid flow regimes originally mapped by the empirical Baker chart, demonstrating the broad versatility of numerical methods.

Despite these advancements, accurately capturing localized transient phenomena remains challenging. As investigated by Ramdin and Henkes [34,35], who established 2D and 3D models to simulate Benjamin bubbles in a horizontal tube and Taylor bubbles in a vertical tube, 2D axisymmetric simplifications were found to deviate significantly from experimental results. They concluded that fully 3D models are imperative when surface tension effects are dominant and the asymmetric 3D deformation of the bubble nose becomes critical. Additionally, fundamental literature has evaluated the kinematic behaviors of elongated bubbles, emphasizing the dependence of bubble propagation on gravity currents in horizontal configurations, as detailed by Zukoski, Benjamin, and Taitel and Barnea [36,37,38].

Vallée et al. [39] conducted experimental investigations and CFD simulations of horizontal stratified two-phase flow, noting that traditional VOF formulations combined with standard RANS turbulence models often fail to resolve fine-scale interfacial fragmentations. In multiphase flows, standard LES sub-grid scale (SGS) models often struggle near the wall and phase interfaces due to excessive turbulent viscosity generation. To address this limitation, the dynamic Smagorinsky-Lilly model calculates the sub-grid coefficients dynamically based on local flow structures, preventing the over-damping of interfacial fluctuations. Ujang et al. [40] observed the detailed initiation and evolution of slugs, pointing out that the highly unsteady velocity profiles within the underlying liquid film and the intense wedge effects are key unmapped factors driving the transient dynamics.

To address the limitations of existing methods, this study introduces a CLSVOF-LES framework. While previous VOF/RANS studies have mapped macroscopic slug behaviors, they often fail to capture micro-scale localized stagnation and interface fragmentation. This paper provides physical insights into the gas wedge mechanism and its direct correlation to localized pressure surges. The practical engineering implications of these findings lie in providing transient structural loading data for subsea pipeline fatigue analysis.

2 Experimental Setup and Numerical Framework

2.1 High-Frequency Interface Detection Physical Platform

To capture the translational motion of the elongated bubble and the transient evolution of the underlying liquid film in a horizontal pipe, this study relies on a high-frequency multiphase flow physical detection platform to obtain baseline data, as illustrated in Fig. 1. The working fluids utilized were air (ρ = 1.225 kg/m3, μ = 1.789 × 10−5 Pa·s) and water (ρ = 998.2 kg/m3, μ = 1.003 × 10−3 Pa·s, σ = 0.072 N/m). The main test section consists of a transparent horizontal polymethyl methacrylate (PMMA) pipe with a total length of 13 m and an inner diameter of 50 mm, which ensures fully developed flow fields and meets the observation requirements for the kinetic energy dissipation of elongated bubbles.

Along the pipeline, 20 groups of non-intrusive double parallel conductance probe arrays are arranged at key nodes. While it is recognized that conductance probes cannot perfectly resolve the precise shape of the elongated bubble close to its highly curved nose [34], they are robust for extracting macroscopic velocity and lengths. The distance between adjacent probes is 0.1 m near the inlet and 0.5 m in the fully developed sections. Compared to traditional macroscopic liquid holdup measurements, this array, combined with a cross-correlation signal processing algorithm, can extract the transit time difference of the head and tail of the elongated bubble at a sampling frequency of 2000 Hz. Specifically, the LabVIEW-based cross-correlation algorithm utilized a sliding window size of 1000 data points (0.5 s) with a 50% overlap percentage to calculate the maximum correlation coefficient R X Y ( τ ) . The time delay Δ t corresponding to the peak R X Y was used to calculate the translational velocity v t = Δ L / Δ t . The measurement uncertainty of the conductance probe system was quantified at ±4.5%, established through static liquid level calibration procedures (using a reference probe for dimensionless voltage conversion).

Furthermore, to continuously monitor the transient hydrodynamic forces, high-frequency dynamic pressure transducers are flush-mounted along the bottom of the test section. Keller pressure sensors (range 0–2 bar, accuracy ±0.1 FS) were installed precisely at z = 12.6 m (L/D = 252) from the inlet to measure local absolute pressure fluctuations. Additionally, Rosemount differential pressure sensors were deployed between z = 7.6 m (L/D = 152) and z = 12.1 m (L/D = 242). Synchronized with the conductance probes at a 2000 Hz sampling rate, these sensors capture the instantaneous axial pressure gradients and the localized pressure surges induced by the violent impact of the liquid slug front and the high-speed gas wedge.

Figure 1: Experimental system.

2.2 CLSVOF Coupled Model

Because the tail of the elongated bubble is accompanied by gas-liquid two-phase shear, wake shedding, and gas entrainment phenomena, the traditional Volume of Fluid (VOF) method is prone to generating parasitic currents and blurring the interface when dealing with large-curvature interface deformations. Therefore, this study introduces the Coupled Level-Set and Volume of Fluid (CLSVOF) method [41] to construct a 3D transient multiphase flow model. This model combines the advantages of the VOF method [30], which satisfies mass conservation, and the Level-Set method [42], which calculates the interface normal vector and curvature.

The advection equations for the phase volume fraction α and the Level-Set signed distance function φ are expressed as follows: ∂α∂t+∇·(αu)=0(1) ∂φ∂t+u·∇φ=0(2) where u represents the velocity vector field of the fluid mixture.

Within the CLSVOF framework, the normal vector n and curvature κ of the interface are directly resolved from the smooth Level-Set function φ , thereby improving the computational accuracy at the gas-liquid boundary of the elongated bubble:

The Navier-Stokes equation governing the momentum transfer of the flow field adopts a single-fluid formulation. The surface tension is transformed into a volumetric force source term and incorporated into the momentum equation via the Continuum Surface Force (CSF) model: ∂(ρu)∂t+∇·(ρuu)=−∇P+∇·μ(∇u+∇uT)+ρg+σκδ(φ)∇φ(4) where the mixture density ρ and mixture dynamic viscosity μ undergo a smooth transition based on the reconstructed Heaviside function. To prevent artificial smearing, the numerical interface thickness was controlled, with the Heaviside transition function spanning exactly three grid cells across the phase boundary.

The computational domain consists of a 13.0 m length pipe (260D) to match the experimental spatial development length. The inlet boundary condition was specified using constant superficial velocities for gas and liquid phases, initialized dynamically as a stratified flow based on a theoretical liquid holdup ratio of 0.5 to trigger natural slugging. It should be noted that this initial liquid holdup ratio was selected based on the theoretical equilibrium level for the tested flow conditions; however, the natural slugging process is governed by the K-H instability mechanism and is not uniquely determined by this initial value. A sensitivity analysis confirmed that varying the initial holdup ratio between 0.4 and 0.6 resulted in less than 3% variation in the fully developed slug characteristics at the measurement sections, indicating that the downstream results are robust to this initialization parameter. This initialization strategy allows the CLSVOF-LES model to explicitly capture the onset of Kelvin-Helmholtz (K-H) instability, as interfacial waves naturally grow and bridge the pipe cross-section to form fully developed slugs. The outlet boundary was set to a fixed atmospheric pressure with a zero-gradient condition for velocity. No periodic boundary conditions were applied.

The numerical simulations were performed using the open-source CFD toolbox OpenFOAM. The governing equations were discretized using the finite volume method. For temporal discretization, a first-order implicit Euler scheme was applied to ensure stability for transient multiphase flows [43]. Spatial discretization of the convective terms in the momentum equation was handled using a second-order upwind scheme, while the diffusion terms were discretized using a second-order central differencing scheme [44]. The advection of the volume fraction and level-set functions was captured using the Multidimensional Universal Limiter with Explicit Solution combined with a van Leer flux limiter to maintain a sharp and bounded interface [45]. The pressure-velocity coupling was resolved using the PIMPLE algorithm, which is a hybrid of the PISO [46] and SIMPLE algorithms, allowing for robust coupling in transient calculations. To properly capture the highly transient interfacial dynamics and ensure numerical stability, the selection of the time step size was guided by the Courant–Friedrichs–Lewy (CFL) condition [47].

2.3 Large Eddy Simulation (LES) Framework

The rapid thinning of the liquid film beneath the elongated bubble and the strong shear vortex shedding at the bubble tail exhibit highly unsteady characteristics. Traditional Reynolds-averaged Navier-Stokes (RANS) models generally average out these critical small-scale turbulent dissipation processes. Consequently, this model innovatively adopts Large Eddy Simulation (LES) to resolve the localized transient flow field. Through spatial filtering operations, the motion of large-scale eddies is directly solved by the filtered Navier-Stokes equations, while the effect of small-scale eddies on the large-scale ones is closed using a Sub-grid Scale (SGS) stress model. In this study, the dynamic Smagorinsky-Lilly SGS model is selected. The sub-grid stress tensor τ i j S G S introduced in the filtered momentum equation is defined as: τijSGS=−2μsgsSij¯+13τkkSGSδij(5) where the sub-grid viscosity μ s g s = ρ ( C s Δ ) 2 S ¯ . The dynamic constant C s adjusts dynamically in time and space according to the local flow field, with recorded fluctuation values ranging between 0.08 and 0.16 during intense interfacial shearing, and Δ represents the filter scale.

2.4 Mesh Generation and Independence Test

To ensure the fidelity of the LES framework in resolving the near-wall turbulent boundary layer and the interfacial dynamics, a structured grid resolution strategy was implemented, as depicted in Fig. 2. A static, locally refined hexahedral mesh was utilized, concentrating grid elements in the phase transition zones and utilizing boundary layer encryption to avoid the computational overhead of Adaptive Mesh Refinement (AMR) while maintaining interface sharpness. The first near-wall grid cell height was controlled to ensure that the dimensionless wall distance, y + , remains predominantly less than 1.0 across all simulated operating conditions.

Figure 2: Mesh generation.

To quantitatively verify mesh independence, three resolutions were tested: Coarse (1.2 million cells), Medium (2.5 million cells), and Fine (4.1 million cells). The percentage difference in bubble translational velocity between the Medium and Fine meshes was 1.2%, while the deviation in pressure gradient was 1.8%. Consequently, the Medium mesh was selected to balance accuracy and computational cost. In addition to spatial grid resolution, a time-step independence test was conducted, which confirmed that dynamically adapting the time step to maintain a maximum Courant number of C 0 < 0.25 provides an optimal balance between temporal accuracy and computational efficiency.

3 Results and Discussion

To characterize the underlying physics, the flow conditions are evaluated using fundamental dimensionless numbers. The mixture Reynolds number ( Re m = ρ m v m D / μ m ) ranges from 45,000 to 120,000, confirming highly turbulent flow. The mixture Froude number ( F r m = v m / g D ) spans from 1.4 to 5.7, indicating a transition from gravity-dominated to inertia-dominated regimes. The Eötvös number ( E o = Δ ρ g D 2 / σ ≈ 340 ) confirms that surface tension effects are secondary to buoyancy and inertia on the macroscopic scale. Nevertheless, although surface tension is globally subdominant at this high Eötvös number, the CLSVOF method remains essential because the accurate computation of interface curvature and normal vectors is critical for resolving localized interfacial phenomena—particularly the fine-scale fragmentation at the bubble tail and the sharp curvature at the bubble nose—where local length scales are significantly smaller than the pipe diameter and the effective local Eötvös number is correspondingly reduced.

3.1 Bubble Shape and Velocity Distribution

The cross-sectional contours of the volume fraction and velocity fields provide insights into the localized structural topologies and macroscopic kinematic behaviors of elongated bubbles. At relatively low superficial gas velocities, the gas-liquid interface of the underlying liquid film remains hydrodynamically stable. Under these conditions, the liquid level exhibits spatial uniformity across consecutive cross-sections from the trailing edge to the leading edge of the bubble, as depicted in Fig. 3.

Figure 3: Cross-sectional contours of volume fraction and velocity field (V_sg = 1 m/s, V_sl = 1 m/s).

Conversely, an increase in the superficial gas velocity induces axial stretching of the elongated bubble, accompanied by the onset of interfacial instabilities and visible fluctuations along the underlying liquid film. Furthermore, the elevation of the liquid film near the bubble nose gradually thickens in the direction of flow, reshaping the frontal region of the gas cavity into a distinct, wedge-shaped profile, as illustrated in Fig. 4. As documented in the literature, at low Froude numbers where gravity dominates, the nose is typically located at the top of the pipe. However, under the high Froude number conditions ( F r m > 3 ) simulated here, inertia forces the nose further down from the top wall, a phenomenon captured by the present model.

Figure 4: Cross-sectional contours of volume fraction and velocity field (V_sg = 2 m/s, V_sl = 1 m/s).

This kinematic evolution is governed by interfacial momentum transfer; the high-speed gas phase exerts an interfacial drag that accelerates the slower underlying liquid layer. As the superficial gas velocity increases and the bubble elongates, the liquid film is afforded a developmental length to be propelled by this interfacial drag. Concurrently, the elongated bubble undergoes deceleration, which leads to the synchronization of the gas and liquid phase velocities at the wedge-shaped nose, as shown in Fig. 4. Additionally, the internal velocity distribution within the gas cavity exhibits spatial heterogeneity. The high-velocity gas core is predominantly confined to the upper cross-section of the pipe, resulting in a velocity gradient where the flow near the upper wall is faster than that in the lower vicinity adjacent to the liquid film. For instance, under Vsg = 2 m/s, Vsl = 1 m/s, the gas core velocity reaches a maximum of approximately 3.8 m/s, while the underlying liquid film abruptly decelerates to below 0.5 m/s, yielding a distinct velocity slip of over 3.3 m/s directly across the interface.

These micro-scale interfacial momentum transfer characteristics and continuous advection mechanisms dictate the macroscopic kinematic evolution of the two-phase flow. To quantitatively substantiate the accuracy of the proposed 3D framework, the numerical results were compared against the physical platform’s probe data. By applying the cross-correlation algorithm to the voltage signals from the parallel conductance probes located at L/D = 258 (Z = 12.9 m), the experimental translational velocity ( v t ) was calculated. The experimental time-averaged length ( L b ) of the elongated bubbles was then mathematically inferred by multiplying v t by the transit duration captured by the probes.

To move beyond qualitative contour observations and quantify the phase velocity slip, the vertical velocity distribution profile at the cross-section of a fully developed elongated bubble was extracted, as shown in Fig. 5. Under the baseline condition (Vsg = 1 m/s, Vsl = 1 m/s), the velocity profile exhibits a discontinuity at the phase boundary. In the upper region of the pipe, the gas core velocity maintains a relatively uniform profile due to the low dynamic viscosity of the gas phase. Conversely, an intense velocity gradient is observed exactly at the gas-liquid interface. The velocity drops across this interface, and the underlying liquid film decelerates towards the bottom pipe wall to satisfy the no-slip boundary condition.

Figure 5: Vertical velocity distribution profile across the cross-section of an elongated bubble.

When the superficial gas velocity is elevated to Vsg = 2 m/s, this vertical velocity differential widens. The upper gas core accelerates, while the lower liquid film velocity increases at a slower rate. This quantified upper-lower velocity discrepancy illustrates the interfacial shear stress. It is this velocity slip that provides the continuous drag required to propel the liquid film forward and trigger the turbulent hydraulic jump at the bubble tail. Although the global mixture Reynolds numbers (Table 1) indicate a turbulent flow regime, the velocity profiles in Fig. 5 deviate from classical fully developed turbulent power-law shapes. This deviation is primarily attributed to phase stratification and the dampening effect of the gas-liquid interface, which suppresses local turbulent fluctuations and imposes a steep velocity gradient that governs the flow structure in this cross-section.

As detailed in Table 1, the simulation demonstrates exceptional fidelity in predicting the phase velocity slip. Under a constant superficial liquid velocity (Vsl = 1.0 m/s), as the superficial gas velocity (Vsg) increases from 0.8 m/s to 2.0 m/s, both experimental and numerical results reveal a substantial acceleration and longitudinal stretching of the bubble cavity. The relative errors for both v t and L b are strictly confined within 6.0%. The unit cell length ( L u , comprising the elongated bubble and liquid slug) and the ratio ( L b / L u ) were also predicted consistently, indicating the model is not unduly sensitive to inlet boundary initialization once the flow reaches the fully developed probe sections. The simulated unit cell length ( L u ) has been explicitly included in the table for comprehensive comparison.

3.2 Pressure Distribution in the Liquid Film

The spatial distribution of the pressure gradient within both the elongated bubble and the underlying liquid film reveals phase-dependent hydrodynamic characteristics. An analysis of the cross-sectional pressure contours, presented in Fig. 6, indicates that, spanning from the tail to the nose of the bubble, the local pressure within the gas core is maintained at a higher level than that in the adjacent liquid film. This transverse pressure gradient is dictated by the dynamic pressure stemming from the gas-phase kinetic energy, coupled with the geometric curvature of the phase interface. It should be noted that the pressure contours presented in this study capture different aspects of the pressure field. Fig. 6 displays the absolute pressure field to illustrate the macroscopic pressure drop toward the atmospheric outlet. Conversely, to better visualize the localized hydrodynamic fluctuations and interfacial forces, subsequent contours are plotted using relative gauge pressure referenced to the outlet boundary. In these relative pressure contours, localized negative values indicate low-pressure wake regions induced by high-speed gas shearing and vortex shedding.

Figure 6: Cross-sectional contours of volume fraction and pressure field (V_sg = 1 m/s, V_sl = 1 m/s).

Furthermore, a continuous pressure drop is observed in both phases along the axial flow direction. Analysis of the axial pressure distribution along the pipe reveals a characteristic sawtooth profile: a linear attenuation trend occurs within the bulk liquid slugs due to frictional resistance, interspersed with abrupt pressure drops directly across the bubble interfaces. This axial pressure attenuation is attributed to viscous energy dissipation mechanisms, which are induced collectively by frictional resistance at the pipe wall and interfacial shear stresses generated between the high-velocity gas core and the slower-moving liquid film.

To quantitatively substantiate the cross-sectional contours, the instantaneous axial pressure profile along the bottom wall (10.0–11.5 m) was extracted, as plotted in Fig. 7. This profile explicitly visualizes the transient hydrodynamic characteristics of horizontal slug flow. Within the bulk liquid slug body, the pressure attenuates with a relatively linear trend dominated by wall friction. However, at the gas-liquid interface near the bubble nose, the abrupt displacement of the liquid film generates a sharp, asymmetric stagnation spike. In the numerical analysis, this stagnation point is dynamically identified as the spatial coordinate where the axial velocity of the liquid phase, relative to the moving bubble nose, approaches zero. This location corresponds to the peak of the spatial pressure gradient along the bottom wall. Notably, elevating the superficial gas velocity to Vsg = 2 m/s amplifies this stagnation impact to approximately 2150 Pa, corroborating the peak amplitude data recorded by the physical transducers in Table 2. Downstream of the nose, the pressure within the gas cavity drops to a flatter, lower-magnitude state, before transitioning into wide-band pressure fluctuations driven by the turbulent hydraulic jump at the bubble tail.

Figure 7: Quantitative axial distribution of relative dynamic pressure along the pipe bottom.

Under elevated superficial gas velocities, as demonstrated in Fig. 8, the unsteady hydrodynamic forces and interfacial interactions exerted on the liquid film become intensified. Driven by the increased gas momentum, the leading edge of the elongated bubble extends and morphs into a wedge-shaped configuration. This high-speed gas wedge imposes a localized impelling and compressive effect on the liquid film situated downstream.

Figure 8: Cross-sectional contours of volume fraction and pressure field (V_sg = 2 m/s, V_sl = 1 m/s).

Consequently, this fluid displacement triggers a localized hydrodynamic stagnation effect, resulting in a pressure surge at the frontal region of the liquid film. Quantitative extraction of the velocity magnitude distribution near the bubble nose shows that the liquid decelerates from the mixture velocity ( v m ) to near zero in the reference frame of the interface. Analysis of the local pressure coefficient reveals that this stagnation spike aligns with the theoretical stagnation pressure estimate ( Δ P ≈ 0.5 ρ l ( v t − v l ) 2 ), corroborating the physical mechanism.

To bridge these micro-scale hydrodynamic mechanisms with macroscopic pipeline performance, the dynamic pressure transducers physically flush-mounted at L/D = 252 were utilized. A crucial distinction must be made between the spatial pressure contours and the time-series pressure data. The spatial contours represent instantaneous snapshots typically captured within the gas core or underlying liquid film, where the relative dynamic pressure remains in the order of hundreds of Pascals. However, in slug flow, the maximum mechanical stress does not occur within the film, but rather at the highly transient impact of the liquid slug front. The time-series data in Fig. 9 specifically captures these transient stagnation pressure spikes. When the high-velocity liquid slug body arrives and abruptly displaces the gas, the localized kinetic energy is converted into a stagnation pressure surge, reaching amplitudes ( Δ P p e a k ) exceeding 2.0 kPa, matching the predictions in Table 2. As illustrated in Table 2, the time-averaged experimental pressure drop (measured by the differential sensors between 152D and 242D) increases with elevated superficial gas velocities, a trend replicated by the CLSVOF-LES simulation.

Figure 9: Comparison of simulation and experiment about pressure over time (V_sg = 2 m/s, V_sl = 1 m/s).

More critically, the experimental dynamic pressure sensors physically detected the localized anomalous pressure surges every time a slug nose passed. To visually demonstrate this transient phenomenon, Fig. 9 compares the simulated and experimental pressure time-series recorded at the downstream sensor location (L/D = 252) under the high-velocity condition (Vsg = 2 m/s, Vsl = 1 m/s). As illustrated, both the numerical model and the physical sensors capture distinct, high-frequency pressure spikes corresponding strictly to the arrival of the slug fronts. The peak dynamic pressure amplitude ( Δ P peak ) of these surges reaches 2.04 kPa.

It should be noted that the baseline shift observed in the experimental data is a physical artifact caused by the backpressure from the downstream gas-liquid cyclone separator, whereas the numerical domain outlets to the atmosphere. Despite this baseline offset, the transient surge amplitudes ( Δ P ) agree well. This physical reading corroborates the localized hydrodynamic stagnation effect identified in the simulated pressure contours. The maximum relative error of 6.34% between the experimental and simulated Δ P peak indicates that the LES framework captures the high-frequency micro-shock waves impinging on the pipeline wall.

As observed in Fig. 9, the experimental and simulated unit cell periods exhibit variability and some discrepancies. This divergence fundamentally stems from the distinct mechanisms governing flow evolution in physical and numerical environments. In the physical experiment, natural slugging is driven by inherent stochastic perturbations, leading to a broader distribution of slug periods. Conversely, the numerical simulation is initialized with an idealized, uniform stratified flow. Consequently, slugging in the CFD model develops strictly through Kelvin-Helmholtz instabilities without external random noise, resulting in more regular and concentrated periodic cycles. Due to this variability, a statistical analysis was performed as shown in Fig. 10. The macroscopic parameters presented in Table 1, specifically the unit cell length ( L u ) and bubble length ( L b )—represent the ensemble-averaged mean values extracted from a statistical sample of consecutive unit cells, rather than values selected from a single arbitrary cycle.

Figure 10: Statistical distribution of experimental and simulated unit cell periods at V_sg = 2 m/s, V_sl = 1 m/s.

3.3 Effect of Gas Velocity on Bubble Morphology

The geometric evolution and spatial topology of the elongated bubble are analyzed by considering the combined mixing Froude number, which encapsulates the interplay of both gas and liquid velocities. As illustrated in Fig. 11, under a constant superficial liquid velocity, an increase in gas kinetic energy induces a longitudinal stretching of the macroscopic bubble structure. Specifically, the morphological transformation at the bubble nose is notable; the augmented gas momentum forces the leading edge to elongate axially, thereby reducing its downward penetration angle into the liquid phase. This adaptation reflects the capability of the high-speed gas core to overcome the opposing liquid inertia, facilitating a transition over the underlying liquid film.

Figure 11: Flow parameter contour plot of longitudinal section (V_sg = 0.8, 1, 2 m/s, V_sl = 1 m/s).

Furthermore, at the trailing edge, elevated gas velocities destabilize the gas-liquid interface, amplifying the intensity of the hydraulic jump. The hydraulic jump interface transitions from a relatively coherent state to an irregular profile, marking the onset of interfacial fragmentation. The interfacial shear stress driven by the accelerating gas flow triggers a micro-entrainment mechanism. Consequently, the tail of the elongated bubble is sheared into a prolonged, stepped configuration that penetrates the adjacent liquid slug. This stepped tail functions as an aeration source, entraining dispersed micro-bubbles into the frontal region of the liquid slug body, which alters the local void fraction distribution and the multiphase mixture dynamics within the pipeline.

3.4 Effect of Liquid Velocity on Bubble Morphology

Conversely, the superficial liquid velocity imposes a compressive and hydrodynamic dampening effect on the morphological evolution of the elongated bubble. Maintaining a constant superficial gas velocity, an increment in the superficial liquid velocity bolsters the liquid phase inertia. This augmented liquid momentum accelerates the liquid bridging process—a key mechanism in slug formation—which effectively compresses the elongated bubble, resulting in a reduction of its overall axial length, as evidenced by Fig. 12.

Figure 12: Flow parameter contour plot of longitudinal section (V_sg = 1 m/s, V_sl = 0.8, 1, 2 m/s).

Beyond this macroscopic compression, the elevated liquid velocity reshapes the localized hydrodynamics at the bubble tail. A higher liquid flow rate inherently thickens the liquid film beneath the gas cavity. This thickened film acts as a hydrodynamic buffer, mitigating the intensity of the hydraulic jump at the trailing edge. As a result, the vertical amplitude of the hydraulic jump is attenuated.

To validate the morphological dampening effect observed in the CFD contours, the instantaneous liquid film thickness was extracted along the axial span of the elongated bubble, as illustrated in Fig. 13. The curves capture the chaotic, irregular interfacial wave signatures inherent to the resolved VOF interface. Under a lower liquid velocity (Vsl = 0.8 m/s), the film is severely compressed to approximately 10 mm. This thin layer is vulnerable to shear, triggering finite-amplitude Kelvin-Helmholtz instability waves that eventually culminate in chaotic fragmentation at the trailing-edge hydraulic jump. Conversely, elevating the liquid velocity to Vsl = 2 m/s structurally augments the underlying film thickness to nearly 18 mm. This quantitatively confirms the formation of a critical “hydrodynamic buffer” layer. As evidenced by the significantly attenuated fluctuation amplitude at the tail region of the profile, this thicker fluid bed successfully dampens the interfacial shear, stabilizing the local hydrodynamics and repressing chaotic gas entrainment.

Figure 13: Quantitative extraction of liquid film thickness profiles along the elongated bubble.

Moreover, this dampening of interfacial fluctuations induces a morphological shift: rather than undergoing vertical fragmentation, the bubble tail elongates into a stable, stepped structure. This extended tail reconfigures the subsequent gas entrainment process, transitioning the gas-liquid mixing mechanism from chaotic interfacial shear to a more stabilized, continuous stratified aeration into the liquid slug body.

3.5 Limitations of the Numerical Model

While the CLSVOF-LES framework captures complex transient dynamics, certain limitations exist. First, the computational cost associated with resolving LES sub-grid scales and sharp 3D interfaces is substantial, rendering it currently impractical for simulating multi-kilometer pipeline networks. Second, the capture of fine interfacial fragmentation is inherently sensitive to grid resolution; droplet shedding smaller than the minimum grid scale cannot be explicitly resolved. Consequently, there are limitations in fully quantifying the secondary breakup of micro-bubbles entrained deep within the liquid slug.

4 Conclusions

In this study, a 3D transient multiphase flow model utilizing the CLSVOF method and a LES framework was developed and experimentally validated to investigate the kinematics of elongated bubbles and liquid film dynamics in horizontal slug flow. The principal conclusions are as follows:

• (1)The proposed CLSVOF-LES model resolves the 3D interfacial topologies and phase velocity slip, predicting the translational velocity and length of elongated bubbles with a relative error of less than 6.0%. As gas velocity increases, the bubble elongates and develops a distinct wedge-shaped nose, driving momentum transfer that synchronizes the phase velocities at the bubble front.
• (2) The high-speed gas wedge exerts fluid displacement on the downstream liquid film, triggering a localized hydrodynamic stagnation effect. This localized interaction induces high-frequency pressure surges at the liquid film front, which serve as the primary catalyst for transient mechanical stresses and pressure pulsations within the pipeline.
• (3) Elevating the superficial liquid velocity thickens the underlying liquid film, acting as a critical hydrodynamic buffer that attenuates the hydraulic jump at the trailing edge. This dampening effect suppresses chaotic interfacial fragmentation, reshaping the bubble tail into a stable stepped configuration and transitioning gas entrainment into a continuous, stratified aeration process within the liquid slug body.