Effect of Changing Crude Oil Grade on Slug Characteristics and Flow Induced Mechanical Stresses in Pipes

: Slug multiphase ﬂow is known to be the most prevalent regime because of its extensive encounters associated with chaotic behaviour, complexity and instability that cause signiﬁcant ﬂuctuations in operating conditions and thus lead to undesirable effects. In this study, the effect of varying crude oil grades on slug characteristics is numerically investigated. A partitioned one-way coupling framework of ﬂuid–structure interaction (FSI) one-way coupling framework is adopted to investigate the inﬂuence of changing oil grades and slug characteristics on the maximum induced stresses in horizontal carbon steel pipe. It was found that increasing crude oil density causes frequent slugging and promotes the formation of liquid slugs further upstream near the inlet with high translational velocity and short wavelength. Therefore, the maximum induced stresses resulting from the interaction between slugs and the inner surface of pipes are strongly dependent on crude oil grade. In modelling extra heavy crude oil, a 40% increase in maximum induced stresses is recorded when the liquid superﬁcial velocity decreases from 1 to 0.86 m/s at a constant natural gas superﬁcial velocity.


Introduction
Transmission pipelines are widely used in the petrochemical, nuclear and petroleum industries as they are reliable, efficient and cost-effective methods for connecting supply and consumers [1]. Specifically, in the petroleum industry, hydrocarbon mixtures are produced from natural oil and gas reservoirs, and they are composed of two or more fluids of different phases. The gas-liquid two-phase flow is the most common type of multiphase flow in most oil and gas systems [2][3][4][5]. One of the main challenges in piping design is the occurrence of slug flow because its extensive encounter in most oil and gas production systems, chaotic behaviour and instability cause undesirable consequences [6,7]. Hydrodynamic slug flow could be developed from a stratified flow regime by the natural growth of small random disturbance in the gas-liquid flow interface or by liquid accumulation due to the instantaneous imbalance caused by hilly terrain pipes. Moreover, the high superficial velocity of the gas phase promotes the coalescence of kinematic waves which acts as an important mechanism for the formation of liquid slugs [8,9]. The occurrence of slug flow is not limited to a specific piping geometry or configuration and may occur in vertical, inclined and horizontal piping [10]. Slug flow is characterised as the alternating flow of liquid slugs and gas pockets. This alternating phenomenon exerts forces on pipe walls and causes stresses and vibration which are known to lead to serious pipeline damage.
Therefore, predicting slug behaviour and quantifying slug characteristics and the consequent stresses are essential steps in flow assurance and piping design. They should be predetermined in the design stage of any production system.
Researchers have investigated slug characteristics analytically, experimentally, and numerically. However, most existing research focused on air-water or gas-low viscous liquid two-phase flows and on determining slug characteristics by varying operating conditions. Nevertheless, computational fluid dynamics (CFD) modelling has become a viable option to predict the behaviour of slug flow and enable researchers to investigate different fluid and gas flow behaviours, including slug flow [11]. Laboratory-sized experiments are widely used for air-water two-phase flows in transparent piping to visualise the multiphase flow behaviour, including slug flow. However, natural gas-crude oil two-phase flows require a highly sophisticated setup and instrumentation for data acquisition to determine multiphase flow characteristics. Hence, CFD modelling techniques are widely used to predict multiphase slug flow behaviour and characteristics. Frank [12] performed the CFD modelling of slug flow by using the commercial CFX software to explore the slug characteristics whilst varying the piping length, inlet and outlet conditions and free surface flow agitation. Frank validated his model against the experimental work of Lex [13]. A reasonable agreement was achieved between the numerical cases of Frank [12] and the experimental work of Lex [13]. Woods and Hanratty [14] attempted to determine the effect of the liquid Froude number on the formation of liquid slugs in air-water flow through horizontal piping by identifying the liquid holdup at multiple locations. They stated that for air superficial velocity between 1 and 4 m/s, slugs were not formed before 40D (D represents the pipe diameter) from the piping inlet. Adibi et al. [15] investigated experimentally the effect of the liquid holdup and variation of gas and liquid superficial velocities on slug initiation. They used the air-water two-phase flow in a rectangular channel and concluded that for liquid holdup H L = 0.25, the slug initiation position occurred farther toward the downstream direction with the increase of the superficial velocities of air and water. For H L = 0.5, it exerted no significant impact on the slug initiation position. Recently, Xu et al. [16] experimentally investigated the formation of slug flow through particle image velocimetry in a fully turbulent air-water two-phase flow in horizontal piping at a constant gas and liquid superficial velocity. They analysed the slug flow characteristics and developed a correlation to predict slug translational velocity on the basis of the maximum local velocity in the slug tail. Doyle et al. [17] investigated the effect of flow rate and physical geometry in CO 2 -water two-phase flow in a developed flow regime. They found that at a relatively low gas superficial velocity, bubble flow is highly likely to occur relative to slug flow. Mohmmed et al. [18] numerically investigated the slug characteristics in air-water two-phase flow in horizontal piping by using STAR-CCM+ software and validated the numerical results against those in the literature [13]. Mohmmed [18] also utilised his validated numerical model to simulate the two-phase flow of air-oil with specific grade to explore the effect of changing the fluids' physical properties on slug characteristics. On the basis of the limited data set utilised in their study, they concluded that liquid physical properties have a significant effect on slug characteristics. Mohmmed [19] performed experimental and numerical investigation to study the effect of varying water and air velocity on the slug behaviour and characteristics in Plexiglass pipe. Mohmmed used STAR-CCM+ software in his numerical investigation. Sam et al. [20] numerically investigated the pressure gradient, slug liquid holdup and slug frequency of oil-gas two-phase flow in horizontal piping. They developed a correlation to predict slug liquid holdup and slug frequency as a function of gas superficial velocity. However, they reported a limitation in which the proposed correlations are only applicable to fluids with similar properties.
Pao et al. [21] numerically investigated the slug liquid holdup and transition velocity for oil-gas oil vapor in horizontal piping and investigated their effect on varying gas superficial velocities. Losi et al. [22] statistically analysed the results from their experimental investigation on highly viscous oil-air slug flow in horizontal piping. They concluded that liquid physical properties greatly influence slug characteristics, and that slug length decreases as liquid viscosity increases. Baba et al. [23,24] experimentally investigated the effect of varying liquid viscosities on slug length and slug translational velocity by utilising two gamma densitometers with fast sampling frequencies. They utilised dimensional anal-ysis and established an empirical correlation to predict slug length and slug translational velocity. Shadloo et al. [25] investigated the effect of changing the operating conditions and physical properties of two-phase flow on pressure drop through the implementation of a multilayer perceptron artificial neural network model. However, their findings are used to validate the results of the developed algorithm. Abdulkadir et al. [26] constructed an experimental facility and utilised nonintrusive instrumentation technique for the measurement of the Taylor bubbles, liquid slugs' translational velocity, slug frequency and liquid slug lengths of an air-silicon oil two-phase flow in a pipe flow rig. They concluded that slug translational velocity is significantly dependent on the superficial velocity of the mixture. As the gas superficial velocity increases, the Taylor bubbles and liquid slug lengths increase. As the liquid superficial velocity increases, the slug frequency increases.
Several researchers have investigated the cascading effect of slug characteristics on the deformation of the mechanical properties of piping by utilising conventional analytical and experimental methods [27][28][29] in which slugs are geometrically idealised as a moving load in pipes with a fixed length, velocity and frequency. However, the slug characteristics were predefined as an input to the model. Chica et al. [30] implemented the two-way coupling of fluid-structure interaction (FSI) to investigate the resultant effect from slug flow into the structural domain by varying the volume fraction and mixture velocity in a subsea jumper for gas-low viscous liquid. In addition, they developed a structured methodology to determine the likelihood of failure and predict the reduction of the designed life of a piping system due to slug flow behaviour. An and Su [31] and Wang et al. [32] analysed the effect of hydrodynamic slug forces on structural behaviour by changing the volumetric gas-liquid fractions and flow rate. They found that slug velocity exerts a significant effect on the vibration response of piping systems. Mohmmed et al. [33] developed a one-way coupled FSI model to establish a prediction model that correlates slug frequency and induced stresses by using air-water mixture flowing in a horizontal Plexiglass piping. Furthermore, Mohmmed et al. [19] also extended his air-water slug flow investigation and studied the effect of slug liquid length and slug translational velocity on the mechanical stresses in a Plexiglas pipe. Liu et al. [1] analytically explored the dynamic responses of gas-liquid two-phase slug flow in a simply supported piping system by varying the piping materials and translational velocities. They concluded that for a high gas superficial velocity, the dynamic response manifests as a periodic-like motion.
Recently, several researchers utilised various intrusive and nonintrusive methods to investigate the flow pattern and identify the two-phase flow characteristics [1][2][3][4][5]. Almutairi et al. [34] used electrical capacitance tomography (ECT) to measure liquid holdup and utilised the high-speed camera images (HSCI) to investigate multiphase flow patterns. Xu et al. [35] combined the typical electrical capacitance tomography (ECT) with the convolutional neural network (CNN) algorithm and predicted the two-phase flow patterns and 'gas-liquid' flow rates. Roshani et al. [36] measured the void fraction and identified the annular regime by utilising gamma radiation attenuation and artificial intelligence. Furthermore, Roshani et al. [37] experimentally determined the flow pattern and predicted the oil-water-gas volume fraction by using sodium iodide crystal detectors, X-ray tube combined with group method of data handling (GMDH) neural network. Fang et al. [38] used ultrasonic phased array technology to identify the 'air-water' vertical flow regimes in multiphase flow.
In sum, most previous studies investigated the effect of varying the operating conditions of slug flow on the induced stresses in pipes whilst keeping the liquid physical properties unchanged. As mentioned previously, slug flow occurs frequently in oil and gas systems and is known to cause serious damage to pipelines. Therefore, the effect of changing the crude oil grade on slug characteristics and the resulting induced stresses should be investigated in detail. In this study, we aimed to investigate numerically the effect of varying crude oil grades on slug characteristics and the structural responses through the implementation of a partitioned one-way coupled FSI technique of a natural gas-crude oil two-phase flow in a horizontal carbon steel (ASTM SA106 Grade B) piping. The developed CFD and structural models were validated against the experimental results for air-water two-phase flow of Mohmmed [19].

Modelling
The numerical model of the fluid and solid domains was developed using Altair AcuSolve and Altair OptiStruct software, respectively. The partitioned one-way coupled FSI technique was adopted by mapping the fluid domain and resultant pressure into the boundary conditions of the structural domain. Both models were validated against the experimental work and numerical modelling results of Mohmmed [19]. Then, the modelling was extended to investigate the effect of varying crude oil grades on slug characteristics and the induced stresses on the gas-crude oil two-phase flow in the horizontal carbon steel piping.

Mathematical Formulation of CFD Domain
For the CFD model in this study, the back-and-forth error compensation (BFECC) algorithm with the level set (LS) method was adopted in the physical modelling in Altair AcuSolve software. The integration of BFECC and LS improves the accuracy of a semi-Lagrangian integration for improved computational performance [39]. Various numerical studies have reported that the implementation of the BFECC algorithm results in secondorder accuracy in the space and time domains and generates a relatively realistic animation of the two-phase flow [40][41][42][43]. The CFD model numerically solves incompressible Navier-Stokes equations as follows: where → V and p, respectively, represent the velocity and pressure profiles, σ and ρ, respectively, denote the fluid stress tensor and density and → f is the summation of the bulk and surface forces per unit volume. In addition, ∇. → V is defined as zero to physically represent the time rate of changing the volume of the flowing fluid element. In the case of fluids with a constant density despite the applied pressure, they adhere to the following equation: where µ is the fluid absolute viscosity and D( → V) represents the strain tensor rate of the fluid. In this work, both fluids satisfied the aforementioned equations. For the multiphase flow simulation, the numerical model was supported with a combination of predefined boundary and initial conditions at the fluid domain and the interface which are intended to solve Equations (1) and (2). In this numerical model, the no-slip boundary condition at the interface was considered to ensure that the defined pressure and velocity profiles were continuous along the interface. The LS method was represented by the Lipschitz scalar function ∅, which is defined as a signed distance function that measures geometrical properties and tracks the transit topological changes in the interface of fluid regions Ω i and Ω j which share the same boundary ∂Ω with an interface denoted by Γ ij . The signed distance function is defined as The LS function reflects the shortest distance from the interface to the grid cell with a positive value for the liquid phase (Ω i ), a negative value for the gas phase (Ω j ) and a zero level implicitly representing the interface (Γ ij ). As the fluid regions Ω i and Ω j evolve with time and satisfy Navier-Stokes Equation (2), the LS function is governed by the following expression: However, an auxiliary equation is required to solve Equation (5) at each time step because of the large variation and the difficulty of preserving the intrinsic property of |∇∅| = 1 as a result of the velocity field in the second term (V.∇∅). Moreover, the curvature of the interface can develop irregularities on the LS function and thereby lead to dissipative errors. Therefore, the BFECC algorithm was introduced to compute the error compensation (e) from the LS function (∅ n ) at time t n by evolving the reversible differential Equation (5) forward in time to measure ∅ n+1 at time t n+1 and backward in time to measure ∅ n 1 at time t n . As the same advection scheme was adopted for the forward and backward operations, the errors were expected to be similar in both steps. Therefore, the numerical error was estimated as follows: This error was subtracted before the final forward advection step, and then, ∅ n+1 was measured. If the forward and backward advection operations returned the exact value, then ∅ n = ∅ n 1 .

CFD Simulation
As mentioned previously, the fluid domain was simulated in Altair AcuSolve software. The horizontal piping had an internal diameter (D) of 74 mm and a length of 8 m at room temperature and atmospheric pressure condition ( Figure 1). The boundary conditions imposed on the domain were as follows: (1) a smooth piping wall with no-slip condition is applied; (2) inflow velocity for both phases is given as an input to the inlet surface; (3) atmospheric outflow pressure is applied at the outlet surface. Initially, the piping was filled with stratified two-phase flow with 50% liquid holdup. A user-defined function was then incorporated at the inlet surface boundary to impose perturbations. The volume fraction of the liquid phase was set as a function of time: The pipe and fluid domain mesh are shown in Figure 2. A fine mesh was placed at the fluid-internal pipe surface interface so that the fluid behaviour and the slug effect on the pipe surface could be accurately observed and determined. The maximum absolute cell length was 2 mm with an aspect ratio of 1.63 and maximum skewness angle of 37.4 • . Mesh independence study was carried out by varying the number of cells. The number of cells was varied to optimise the number of cells and to ensure that mesh independence was reached. As shown in Figure 3, increasing the number of elements to more than 388,504 resulted in a 0.5% change in the slug frequency values. Therefore, we set the optimum number of elements in this study to 388,504.

Mathematical Formulation of Solid Domain
An explicit dynamic solver was employed for the finite element (FE) model, which utilises the widely used Newmark integration scheme for the structural response analysis to determine the induced stresses caused by the slug flow. The explicit finite element equation of motion is expressed as where M, C and K are constant values that represent the mass matrix, damping matrix and stiffness matrix, respectively; q R u represents the coordinates used in the mechanical system configuration; F is the load vector at the t th time step; .. q .
, q and q are the acceleration, velocity and displacement, respectively.
The solution was derived through a one-step time integration by solving the equation of motion in the Newmark scheme. The results were formulated in terms of time history. The vector iterations were calculated through the following system: ..
, q 0 and q 0 are the acceleration, velocity and initial displacement, respectively. The resultant load at time t was determined by The displacement q t+∆t at time t + ∆t was calculated using the following expression: The dynamic stress σ t+∆t at time t + ∆t can be predicted from the displacement at time t + ∆t through the following expression: The Young's modulus, Poisson's ratio and stain displacement matrix are denoted by E, ν and B, respectively.

Solid Model
The developed model was validated against the experimental measurement of Mohmmed [19]. The developed solid model's material properties which are similar to the Plexiglass piping specifications of [19], are listed in Table 1.  Figure 4 shows that fixed support boundary conditions were adopted at the inlet and outlet surfaces to simulate the operating conditions of the experimental work of Mohmmed [19]. The resultant pressure profile from the fluid model domain was extracted from the CFD simulation and mapped in the solid structure. A mesh independence study of the solid model was conducted. The number of elements was varied, as shown in Figure 5. Increasing the number of elements to more than 86,880 resulted in a 0.21% change in the value of the maximum induced stresses predicted by the present model. Therefore, the optimum number of elements of 86,880 was used in the present FSI model. The developed model showed a 10% variation against the experimental measurement of Mohmmed [19]; hence, the results of both works could be considered to have a reasonable agreement. The developed model could thus be used to predict the induced stresses in carbon steel ASTM SA106 Grade B pipes, which are widely used in oil and gas plants.

Results and Discussion
This section is divided into two parts, namely, CFD modelling and solid modelling. As mentioned previously, Altair AcuSolve and Altair OptiStruct are respectively used in the CFD modelling for the fluid domain and the solid modelling for the one-way coupled FSI.

Numerical Model Validation
As mentioned, this research investigates varying crude oil grade effect on slug characteristics and induced stresses. Due to the fact, and to the authors best knowledge, there is no published research that investigates the effect of crude oil slug flow on induced stresses in horizontal pipes. Hence the results of 'water-air' slug flow published by [19] was used to validate this model which could be considered as reasonably acceptable validation; taking into consideration that the software solvers (CFD and FE) used in this research solve the same governing equations for both crude oil-natural gas and 'water-air' multiphase flow. Past research also performed their validation and verification models of crude oil against 'water-air' multiphase published results [18,30].

CFD Fluid Domain Model Validation
The developed CFD model was further validated quantitatively by comparing the results of the predicted slug characteristics (slug initiation position, slug translational velocity, slug liquid length and slug frequency) for the cases shown in Table 2 with the measured results by Mohmmed [19]. Selected cases (1, 2 and 4) were used to validate qualitatively the present model with the slug contours of the STAR-CCM+ model and slug images from Mohmmed [19]. The results of the qualitative analysis of slug formation and growth and the liquid length of the present model were compared with the slug contours in the experimental work and numerical model of Mohmmed [19] (Figure 6). The captured slugs present different hydrodynamic behaviours and liquid lengths, which change as a function of gas and liquid superficial velocities. The slug flow for j L = 1 m/s and j G = 2.44 m/s shows the shortest liquid slug which increases in length as the water superficial velocity decreases, as in j L = 0.86 m/s and j L = 0.7 m/s. As shown in Figure 6, the qualitative results of the developed CFD model are in good agreement with those of the experimental work and CFD model of Mohmmed [19].  Figure 7 shows the quantitative validation of slug initiation position, slug translational velocity, slug liquid length and slug frequency against those reported by Mohmmed [19]. In general, the slug characteristics are strongly dependent on the water and air superficial velocities. As the water superficial velocity increases, the slug formation tends to occur further downstream from the piping inlet. Moreover, a reduction in the slug liquid length is observed as the water superficial velocity increases. The comparison of the results of the present model with those of Mohmmed [19] shows that the highest deviations recorded for the slug initiation position, slug translational velocity, slug liquid length and slug frequency are 11.8%, 11.5%, 10.8% and 11.8%, respectively. Mohmmed [19] validated his experimental work by using the STAR-CCM+ numerical model [19] and found that the highest deviation recorded for slug length was 17.6% when j G = 2.44 m/s and j L = 0.86 m/s. Therefore, relative to the experimental work of Mohmmed [19], the present model achieves a reasonable prediction of slug characteristics.

FSI Model Validation
The maximum induced stresses due to the effect of the slug flow characteristics on the piping were used to validate the FSI model in the current work. Cases 1, 2 and 4 from the CFD model were utilised for the FSI model validation, and the results were compared against the findings of Mohmmed [19]. Figure 8 illustrates the maximum induced stresses plotted against the water superficial velocity and the quantitative comparison. Through the comparison of the results from the present model with those by Mohmmed [19], the deviations for Cases 1, 2 and 4 are recorded as 9.5%, 8.7% and 8.8%, respectively. The deviations indicate the reasonable agreement of the results compared.

Effect of Varying Crude Oil Grades on Slug Characteristics
After performing the validation, the model was modified to include the effect of varying the crude oil grade (density and viscosity) on the slug characteristics in the carbon steel pipe. Table 3 lists the properties of the natural gas-crude oil two-phase flow and the modelled cases' superficial velocities. In practical applications, the oil superficial velocity lies in the range of 0.1-1 m/s whilst the natural gas superficial velocity lies between 0.5 and 4 m/s [44]; these values explain the selection of the cases under investigation herein. Table 3. Physical properties of two-phase flow. 1 Natural gas phase density (ρ L ) and viscosity (µ G ) is 0.9 kg/m 3 and 0.02 cP, respectively. Figure 9 illustrates the slug initiation position for the modelled cases. Under the constant superficial velocities of gas and crude oil, an increase in crude oil density causes the slug initiation position to decrease. Moreover, the slug initiation position for light crude oil increases by 25% when the liquid velocity increases at a constant gas velocity j G = 2.44 (m/s). As the density of crude oil increases, the location of the slug initiation position moves closer to the inlet and decreases by 84% for very heavy crude oil, 48% for the heavy crude oil and 16% for medium crude oil relative to light crude oil.  Figure 9 also shows that the variation in slug initiation position is minimal when the crude oil grade increases and is highly noticeable for light crude oil. Meanwhile, decreasing the gas velocity to 0.7 m/s at 0.93 m/s liquid velocity results in a noticeable increase in slug initiation position for light-grade crude oil; this increase is minimal for very heavy crude oil. The contours shown in Figure 10 indicate that for light crude oil (low-viscosity crude oil), the slug initiation occurs further toward the downstream distance. As the crude oil grade increases at a constant gas-liquid velocity, the slug formation occurs close to the pipe inlet.

Slug Liquid Length
Slug liquid length is one of the main characteristics of slug flow, and an incorrect prediction may cause flooding in the downstream equipment and thereby result in production deferment. Figure 12 illustrates the effect of changing the oil grades on the slug liquid length. For all investigated cases, increasing the crude oil density decreases the slug length. Moreover, decreasing the liquid velocity at a constant gas velocity causes an increase in slug length. By contrast, up to 30% decrease in slug length is recorded when the gas velocity decreases at liquid velocity of 0.93 m/s as compared with Case 1.

Slug Translational Velocity
The effect of varying crude oil grades on slug translational velocity is presented in Figure 13. The slug translational velocity for high-density crude oil is higher than that for low-density oil. This difference can be explained by considering the mass per unit volume of liquid slug. As shown in Figure 12, higher density crude oil has shorter slug length. Although the slug length decreases, the decrease is larger than the increase in oil density. Hence the mass of the slug for higher density oil is less than the mass of lower oil density which causes the slug translation velocity to increase.  The slug frequency increases when the crude oil density increases. Moreover, as the density of crude oil increases, the difference in slug frequency values decreases at different fluid and gas velocities. For instance, the slug frequency for very heavy crude oil fluctuates between 1.89 and 2 Hz for all examined superficial velocities.
For light crude oil, the slug frequency is highly influenced by the superficial velocity of the gas and liquid phases, and the variation is more noticeable than that of heavy crude oil. When the gas and liquid superficial velocities are 0.7 and 0.93 m/s, respectively, the slug frequency for light crude oil is 1.17 Hz. As the gas and liquid superficial velocities increase to 2.44 and 1.0 m/s, respectively, the slug frequency noticeably decreases to 0.6 Hz. By maintaining the gas phase superficial velocity at 2.44 m/s and decreasing the crude oil superficial velocity to 0.86 m/s, the slug frequency further decreases to 0.53 Hz.

Effect of Varying Crude Oil Grades on Maximum Induced Stresses
Slug flow occurrence is known to cause serious damage to pipelines. Hence, the effect of slug characteristics on induced stresses in carbon steel pipes is investigated in this work. Figure 15 shows that varying crude oil grades cause significant variation in maximum induced stresses values. The flow of light crude oil records the lowest induced stresses. However, as the crude oil density increases, the induced stresses increase.  Figure 15 also shows that the flow of light and medium crude oil in carbon steel pipe records a marginal increase in stresses when the liquid and gas velocities change. By contrast, when heavy and extra heavy crude oil flow in pipes, the margin of increase in induced stresses due to changes in gas and liquid velocities is significant.
Moreover, for very heavy crude oil, the maximum induced stresses decrease by up to 33% when the liquid velocity increases at a constant gas velocity. Meanwhile, decreasing the gas velocity at 0.93 m/s liquid velocity records the lowest induced stresses for all crude oil grades.
The slug initiation position is a function of the operating conditions and physical properties of flowing media. Slug formation in piping systems subsequently exerts hydrodynamic loads on the piping's inner surface. The perturbation of the slug that moves with a velocity higher than the gas and liquid superficial velocities causes high momentum, which is then transferred to the piping wall where it causes induced stresses. Figure 16 illustrates the relation between the slug initiation position and the induced stresses for Case 1. As the slug formation occurs at 64D (D represents the pipe diameter) for the light crude oil, the maximum induced stress observed is 48 MPa at the location of slug formation. When the crude oil density increases, the liquid slug translational velocity increases (Figure 13), thereby causing the slug to move with high momentum. This condition leads to high induced stresses, as seen in the medium crude, heavy and very heavy crude oil grades with maximum recorded induced stresses of 54, 189 and 216 MPa, respectively, as shown in Figure 16. Slug liquid length is one of the main variables that characterise slug flow and could be associated with maximum induced stresses when the liquid impinges on the internal wall of the pipe. Thus, the relationship between slug length and maximum induced stresses for various crude oil grades should be understood. Figure 17 shows that for all cases (Cases 1, 2 and 3), when the slug length increases, the recorded maximum induced stresses decrease. For heavy and extra heavy crude oil at a constant gas velocity, increasing the liquid velocity causes a decrease in slug length and consequently leads to high stresses. A possible reason is that the stresses are related to both the magnitude of force and the area. Hence, the induced stresses caused by the slug flow is affected by the fluid density and mass of the slug and the size of the crude oil slug impinged area at the pipe surface. However, decreasing the gas velocity to 0.7 m/s at 0.93 m/s liquid velocity records the lowest induced stresses for all investigated cases. It is important to mention that this research mainly investigates varying the liquid velocity at constant gas velocity.
The results of the induced stresses indicate the need to carefully perform pipeline designs that handle heavy and extra heavy crude oil because the maximum induced stresses may change significantly after a slight change in liquid superficial velocities. Moreover, when the slug moves at a velocity that is higher than the gas and liquid superficial velocities, the slug moves with higher momentum, resulting in higher maximum induced stresses.

Conclusions
The effect of varying crude oil densities on slug characteristics and consequently on induced stresses was numerically investigated using gas-crude oil two-phase flow in a horizontal carbon steel pipe.
The variation of crude oil properties exerts a significant effect on the slug characteristics of the various cases investigated.
For light crude oil, a 14% decrease in slug translational velocity, 16% decrease in slug length and 13% increase in slug frequency are recorded when the liquid superficial velocity increases at a constant gas velocity. Moreover, increasing the crude oil density results in an increase in maximum induced stresses. Meanwhile, varying liquid and gas velocities causes a small variation in maximum induced stresses when light and medium crude oil are modelled. When heavy and extra heavy crude oil are modelled, the variation in induced stresses is significant. Hence, when flow operating conditions change, the design of pipelines that handle heavy and extra heavy crude oil should be carefully performed because small variations in liquid flow velocity may cause significant changes in maximum slug-induced stresses.

Data Availability Statement:
The data presented in this study are available on request from the corresponding authors. The data are not publicly available due to privacy.