Abstract

Wax deposition in field-scale crude oil pipelines poses a significant challenge to the oil and gas industry, leading to reduced flow rates, increased pressure drops, and potential blockages. Understanding the mechanisms governing wax deposition is crucial for developing effective mitigation strategies. This study investigates the impact of multiphase flow conditions, including water-in-oil emulsion, wax precipitation kinetics, shear dispersion, and molecular diffusion, on wax deposition in field-scale crude oil pipelines. A numerical model is developed that employs second-order semi-implicit temporal discretization schemes, such as Crank–Nicolson and Adams–Bashforth methods, in conjunction with a bivariate spectral collocation scheme using Chebyshev–Gauss–Lobatto grid points. The impact of various flow parameters, including Reynolds number (Re), mass Grashof number (Gr), Schmidt number (Sc), and Weber number (We), on the flow variables, wall shear stress, and heat and mass fluxes are investigated. The numerical simulations demonstrate that flow parameters significantly influence the flow behavior, wall shear stress, wall heat flux, and wall mass flux in waxy crude oil pipelines. Specifically, the aggregation of wax crystals in the pipeline decreases by at most 2.5% with increasing Reynolds number from 2.2361 to 3.1361. Conversely, it increases by at most 3.4% with increasing mass Grashof number from 5 to 11 and by at most 4.8% with increasing Weber number from 1.0 to 2.5. Furthermore, the Nusselt number increases from 1.9907 to 4.9834 with increasing Reynolds number from 2.2361 to 5.2361 and from 1.9907 to 2.0225 with increasing mass Grashof number from 5 to 20. It also increases from 1.9907 to 2.0434 with increasing Weber number from 1.0 to 2.5. The insights gained from this study can be applied to optimize pipeline design, operational parameters, and wax deposition mitigation strategies, leading to enhanced pipeline performance and reduced operational costs. The numerical model developed in this work serves as a valuable tool for simulating and predicting wax deposition behavior under various operating conditions.

1. Introduction

Wax deposition in pipelines, caused by the solidification and adhesion of paraffin molecules at temperatures below the wax appearance temperature, leads to flow restrictions, reduced throughput, and increased operational costs. Numerical studies of wax deposition in large-scale crude oil pipelines are crucial for understanding and mitigating these challenges [1]. These studies employ computational models and simulations to analyze the complex interplay between fluid flow, heat transfer, and mass transfer during wax deposition, helping pipeline operators to optimize operations, design effective wax control strategies, and plan maintenance schedules to prevent excessive wax deposition and its associated problems. Additionally, these studies facilitate the development of cost-effective measures for wax removal and pipeline cleaning, and they help in tackling complex engineering challenges, as evidenced by the work of Ali et al. [2] on the Carreau fluid model in roll coating, Ali et al. [3] on Levenberg–Marquardt artificial neural networks for reverse roll coating, and Ali et al. [4] on perturbation-based solutions for non-Newtonian fluids in roll coating.

The bivariate spectral collocation method (BSCM) is a powerful numerical approach for solving partial differential equations (PDEs) involving two independent variables [5]. It combines the advantages of spectral methods, known for their high accuracy and convergence rates, with the ability to handle complex geometries and boundary conditions. In BSCM, the solution to the PDE is approximated using a set of globally defined basis functions, typically Lagrange interpolating polynomials, over a rectangular or nonrectangular domain. The collocation points, where the solution is evaluated, are chosen carefully to ensure optimal accuracy. BSCM has been successfully applied to a wide range of problems in fluid dynamics, heat transfer, and other scientific disciplines.

In recent years, various experimental and theoretical studies have been conducted to model the formation and deposition of solid wax crystals on the inner surface of crude oil pipelines. These studies have received immense attention due to their importance in understanding and mitigating wax deposition challenges. For instance, Kim et al. [6] modeled the consolidation of wax deposition for a progressive cavity pump using computational fluid dynamics. Mrinal et al. [7] investigated a transient 3D computational fluid dynamics model of a progressive cavity pump. Waheed and Megahed [8] studied the heat transfer mechanism of the non-Newtonian micropolar slip fluid flow over a stretching sheet in the presence of the melting heat transfer to heat generation or absorption in the slip flow regime. Additionally, Singh et al. [9] conducted a numerical study on the formation and aging of the wax deposit by performing experiments using laboratory flow loops and considering pipelines with externally cooled walls.

Other significant contributions to the understanding of wax deposition in crude oil pipelines include: Stubsjoen’s exploration of both numerical and analytical modeling of paraffin wax in crude oil pipelines, as reported in Stubsjøen’s [10] study; Fusi’s numerical study of waxy crude oil flow in a laboratory test loop, as documented in Fusi’s [11] study; Banki et al.’s investigation of the numerical modeling of wax deposition in oil pipelines for the laminar flow regime based on the enthalpy–porosity approach, as presented in Banki et al.’s [12] study; Zhang et al.’s study on wax deposition and development of a model to predict the temperature profile and location of wax deposits based on the coupling process involving heat and fluid, as detailed in Zhang et al.’s [13] study; Ying et al.’s experiment to analyze the heat transfer of oil phase change in the case of overhead pipeline shutdown, as described in Ying et al.’s [14] study; and Magnini and Matar’s investigation of the deposition of wax in crude oil pipelines through interface-resolved numerical simulations, as documented in Magnini and Matar’s [1] study. These studies have collectively provided valuable insights into the complex phenomenon of wax deposition in crude oil pipelines.

Previous research on wax deposition in crude oil pipelines has largely focused on single-phase flow scenarios, limiting the applicability of existing wax deposition models to real-world multiphase flow scenarios encountered in the petroleum industry. These multiphase flows, including oil–water two-phase and water–gas–oil three-phase flows, are becoming increasingly prevalent due to the rising water content in oil-bearing rocks and the presence of water-in-oil emulsions during crude oil extraction. To address these research gaps, this study focuses on developing numerical solutions for waxy crude oil flow in pipelines with heat and mass transfer, aiming to extend the applicability of wax deposition models to real-world oil pipelines by incorporating the complexities of multiphase flow scenarios.

The novelty of this research lies in the comprehensive numerical investigation of wax deposition in oil–water two-phase flow in the presence of water-in-oil emulsions, incorporating critical factors such as the porosity of the deposited wax layer, the influence of surface tension-induced forces, the effects of internal heat absorption or generation, the kinetics of wax precipitation, the impact of energy flux due to viscous dissipation, and mass flux due to molecular diffusion and shear dispersion. This comprehensive approach provides a deeper understanding of the complex wax deposition phenomenon in a more realistic multiphase flow environment. The remaining sections of the paper are structured in the following manner: Section 2 presents the formulation of the mathematical model for the flow of waxy crude oil in pipeline systems; Section 3 presents the numerical techniques used to solve the model equations; and Section 4 presents the study results and a comprehensive discussion of the results. The results are validated in Section 5. Finally, the summary and conclusions drawn from this study are outlined in Section 6.

2. Mathematical Formulation

In this study, we investigate the 2D unsteady flow of waxy crude oil within a model of a pipeline with a semi-infinite length, circular cross-section, and an inner radius denoted as , as depicted in Figure 1. The pipeline is inclined at an angle denoted as with respect to the horizontal. We utilize a cylindrical coordinate system denoted as , where represents the radial distance measured from the pipeline’s central axis, represents the circumferential direction, and indicates the axial direction. At the initial time , waxy crude oil at a uniform temperature is introduced at the pipeline inlet. The inner pipeline surface, assumed to be smooth, impermeable, and rigid, is maintained at a constant temperature . The temperature can either be higher or lower than .

The interface between the discrete solid phase and the continuous fluid phase is sharp. The outward-drawn unit normal vector to this interface is given by . The fluid phase comprises three pseudo-components: oil, wax, and emulsions.

We employ the pseudo-single phase approach, which treats the three-phase mixture of water, oil, and gel as a single fluid system. The physical properties of this pseudo-fluid are determined by averaging the relevant physical properties of water, oil, and gel. This averaging is weighted according to their respective volume fractions, as given in Yang et al.’s [15], Zheng et al.’s [16], Al-Ahmad et al.’s [17], and Ochieng et al.’s [18] studies:where the subscript denotes “mixture fluid” and denotes “fluid phase.”

The following assumptions are adopted in this study: no gas is present in the pipeline, the flow is axisymmetric, molecular diffusion and shear dispersion are the sole mechanisms of wax deposition, the fluid particles exhibit no slip at the interface between the fluid and solid phases, and the thermophysical properties are assumed to be constant with the exception of the temperature—and concentration-dependent density variation included in the body force term. Therefore, the Boussinesq approximation is employed to model the flow within the boundary layer.

Utilizing the assumptions outlined above, we obtain the following dimensionless equations [18]:

Equation of continuity:

Equations of conservation of linear momentum:

Equation of energy:

Equation of wax concentration:

Equation of wax precipitation kinetics:

Equation of oil volume fraction:

Equation of deposit growth:

Equation of deposit aging:

The following dimensionless numbers and parameters are employed in the model Equations (13)–(21):

The following thermodynamic model, presented in Cragoe’s [19] study and Al-Ahmad et al.’s [17] study, is adopted in this study:where denotes the shift factor while denotes the molecular weight of waxy crude oil. Note that the American Petroleum Institute (API) gravity is a scale used to grade crude oils, calibrated in degrees API (API). Thus, crude oils are classified into three categories based on their API gravity: heavy crude oils are those whose API gravity is less than 22.1°API, intermediate crude oils are those whose API gravity ranges between 22.1°API and 31.5°API (inclusive), and light crude oils are those whose API gravity is greater than 31.5°API. In this study, a heavy waxy crude oil with an API gravity of 18°API is considered.

The corresponding boundary and initial conditions for the flow are formulated in the following manner:

Here, the parameter , , and . Note that when , there is no gel layer present.

This study focuses on the following engineering parameters of interest: the skin friction coefficient (), the local Nusselt number (), and the local Sherwood number (). These parameters are expressed in dimensionless form as follows:

The skin friction coefficient, the local Nusselt number, and the local Sherwood number characterize the shear stress, the rate of heat transfer, and the rate of mass transfer at the wall of the crude oil pipeline, respectively.

3. Numerical Solution

The governing equations given by Equations (13)–(21), subject to the boundary and initial conditions given in Equation (33), are solved numerically using a combination of the second-order semi-implicit finite difference method [20] and the BSCM [5]. The model equations are discretized in two separate stages: temporal discretization and spatial discretization, as shown below.

3.1. Temporal Discretization

The nonlinear terms are discretized in time using the Adams–Bashforth second-order method, which is an explicit numerical scheme with second-order accuracy. The linear terms involving spatial derivatives are discretized in time using the Crank–Nicolson method, which is an implicit numerical scheme with second-order accuracy. This decoupling of the governing PDEs helps save on computer memory. Finally, the forward Euler method is used to discretize the time derivatives. Therefore, the time-discretized numerical schemes are as follows (for ):

The arbitrary functions and are given by the following equations:

The intermediate functions , , , , , , , , and are expanded using first-order Taylor series about the point to get the following equations:

The time-discretized boundary and initial conditions are given as follows:

The linear iterative schemes given by Equations (37)–(54) and the corresponding boundary and initial conditions given by Equations (55)–(59) are discrete in time but continuous in space.

3.2. Spatial Discretization

The linear iterative schemes given by Equations (37)–(54) are discretized in space using the BSCM based on the Chebyshev–Gauss–Lobatto grid points [5]. The domain is transformed into the new domain using the following linear transformation:

Similarly, the domain is transformed into the new domain using the following linear transformation:

Here, , and , where is a sufficiently large finite number to approximate the asymptotic behavior at infinity. Lagrange fundamentals (or Lagrange coefficients) are chosen as the basis functions. We use the bivariate Lagrange interpolating polynomials to approximate the unknown functions , , , , , and as follows:

The Lagrange cardinal polynomials are defined by the following equation:

The above interpolations utilize symmetrically distributed Chebyshev–Gauss–Lobatto grid points defined on the domain by the following equation:for , where and denote the number of collocation (or grid) points in and direction, respectively. The Chebyshev–Gauss–Lobatto grid points are indexed from right to left of the domains in and since , , , and . Hence, we take and as the computational grids.

The Chebyshev–Gauss–Lobatto grid points are selected to discretize the continuous spatial derivatives (in both and ) and convert them into a discrete matrix form at the collocation points using the standard Chebyshev derivative matrices and [21]. This is illustrated as follows:where

The partial derivatives of the other dependent variables, i.e., , , , , and , with respect to and are similarly transformed to discrete matrix form. Substituting the respective discrete derivative matrices into the temporal schemes above yields matrix systems of the form:where is the coefficient matrix, is the unknown column vector, and is the solution matrix. The corresponding boundary conditions are imposed on the main diagonal of the subblock matrices of .

The system represented by Equation (79) is solved iteratively starting from appropriate initial guesses. The iteration process is repeated for until the prescribed absolute error tolerance is achieved. MATLAB® software is employed for the computer simulations.

4. Results and Discussion

The flow variables investigated in this study include axial velocity, radial velocity, fluid temperature, wax concentration, wax precipitation kinetics, and oil volume fraction. Various flow parameters were varied, including Reynolds number (), thermal Grashof number (), mass Grashof number (), Eckert number (), Schmidt number (), and Weber number (), at the final time step. These parameters were input into a computational program for independent variation. The results of the parametric study are presented in graphs and tables and subsequently discussed.

4.1. Effects of Varying Reynolds Number

It is observed in Figure 2 that an increase in Reynolds number leads to an increase in the radial velocity profiles of waxy crude oil in the pipeline under laminar flow conditions. The Reynolds number represents the ratio of inertial forces to viscous forces acting on a fluid element. The observed trend can be attributed to the enhanced shear stress acting on the fluid. As the Reynolds number increases, the shear stress, which represents the frictional force between adjacent fluid layers, becomes strong enough to induce mixing between the layers. This mixing, known as shear dispersion, facilitates the distribution of wax particles more evenly across the pipe radius. The increased shear stress enables the waxy molecules to overcome the cohesive forces that tend to aggregate them, allowing them to slide past each other more easily. Consequently, the tendency of wax particles to deposit on the pipe walls decreases, and the radial velocity profiles become more uniform across the radial direction.

It is observed in Figure 3 that an increase in Reynolds number results in a decrease in the axial velocity profiles of waxy crude oil in the pipeline under laminar flow conditions. This trend can be attributed to the enhanced drag force acting on the fluid. As the Reynolds number increases, the drag force, which opposes the movement of fluid particles, also intensifies. The increased drag force causes the fluid particles to decelerate in the axial direction, leading to a reduction in the average axial velocity. Additionally, with an increase in Reynolds number, the axial velocity profile becomes more nonuniform. The slower moving fluid particles near the oil–gel interface exert a retarding effect on the faster moving fluid particles at the pipe centerline, contributing to the decline in average axial velocity.

It is observed in Figure 4 that an increase in Reynolds number leads to an increase in the temperature profiles of waxy crude oil in the pipeline under laminar flow conditions. As the Reynolds number increases, the shear stress, which represents the frictional force between adjacent fluid layers, also intensifies. This increased shear stress causes the waxy crude oil molecules to rub against each other with greater force, generating heat due to internal friction. The generated heat is then transferred to the surrounding fluid, resulting in an overall increase in temperature. Furthermore, with an increase in Reynolds number, the thickness of the thermal boundary layer, which is a thin layer of fluid adjacent to the oil-deposit interface where temperature gradients are significant, decreases. This reduction in boundary layer thickness can be attributed to the increased shear stress, which disrupts the stagnant layer of fluid near the pipe wall and promotes more efficient heat transfer. The thinner boundary layer allows heat to dissipate more effectively from the oil-deposit interface into the bulk fluid, contributing to the overall increase in temperature profiles.

It is observed in Figure 5 that an increase in Reynolds number leads to an increase in the total concentration profiles of wax in crude oil pipelines under laminar flow conditions. This trend arises because higher Reynolds numbers signify stronger fluid motion driven by inertial forces relative to viscous forces. As the Reynolds number increases, the waxy crude oil encounters greater fluid movement in the radial direction, promoting mixing within the pipeline. This increased flow facilitates the distribution and dispersion of waxy components more evenly throughout the fluid, resulting in higher total concentration profiles. In essence, the enhanced fluid motion counteracts the tendency of waxy components to settle or adhere to the oil–gel interface, leading to a more uniform distribution of waxy crude oil within the pipeline.

It is observed in Figure 6 that an increase in Reynolds number leads to a decrease in the aggregation degree profiles of wax crystals in crude oil pipelines under laminar flow conditions. This trend is attributed to the stronger fluid motion driven by inertial forces relative to viscous forces at higher Reynolds numbers. When the Reynolds number increases, the fluid’s enhanced radial velocity and turbulence promote better mixing and dispersion of waxy components. This vigorous flow hinders the waxy particles from sticking together or forming aggregates, resulting in lower aggregation degree profiles. Essentially, the enhanced fluid motion disrupts the tendency of waxy components to clump or adhere to each other, thereby reducing aggregation within the pipeline.

It is observed in Figure 7 that an increase in Reynolds number leads to a decrease in the volume fraction occupied by crude oil in the pipeline under laminar flow conditions. As the Reynolds number increases, shear stress between fluid layers also increases. This increased shear stress causes the waxy crude oil molecules to deform and flow more readily, promoting the dispersion of wax crystals within the fluid. The increased shear stress breaks up the agglomerates of wax crystals, causing them to become more uniformly dispersed within the fluid, resulting in an overall increase in wax concentration. This increase in wax concentration leads to a reduction in the volume fraction of crude oil. Moreover, the increased shear stress also minimizes wax deposition on the oil–gel interface. As the waxy crystals are more effectively dispersed, they are less prone to precipitation and adhesion to the pipe wall, reducing the accumulation of wax deposits. This reduced deposition enables more wax crystals to remain in suspension within the bulk fluid, further contributing to the reduction in the volume fraction of crude oil. This phenomenon bears practical implications for the efficient transport and processing of crude oil in pipelines, as it directly impacts the quality and composition of the transported fluid.

4.2. Effects of Varying Mass Grashof Number

It is observed in Figure 8 that an increase in the mass Grashof number leads to a decrease in the radial velocity profiles of waxy crude oil in the pipeline. Mass Grashof number indicates the ratio of the buoyancy forces acting on the wax crystals to the viscous hydrodynamic forces. The observed trend arises from the enhanced buoyancy forces acting on the wax crystals associated with higher mass Grashof number. This increased buoyancy forces cause the wax crystals to rise toward the oil–gel interface, forming a layer of wax-enriched fluid with higher viscosity compared to the surrounding bulk oil, leading to a reduction in its flow rate. Therefore, the velocity boundary layers adjacent to the oil–gel interface thicken. The thickened boundary layer impedes radial fluid motion, consequently reducing the radial velocities across the pipe’s cross-sectional area.

It is observed in Figure 9 that an increase in the mass Grashof number leads to an increase in the axial velocity profiles of waxy crude oil in the pipeline. As the mass Grashof number increases, the species buoyancy forces become more dominant relative to viscous hydrodynamic forces. Consequently, the buoyancy of the waxy crude oil increases, making it more likely to move upward within the pipeline. The upward movement of waxy crude oil induces a flow of oil in the axial direction of the pipeline. This axial flow thins the hydrodynamic boundary layer, the layer of oil adjacent to the pipe wall that is slowed down due to friction. The thinned hydrodynamic boundary layer facilitates the unrestricted flow of waxy crude oil in the axial direction, resulting in an increase in the axial velocity profile.

It is observed in Figure 10 that an increase in the mass Grashof number leads to a decrease in the temperature profiles of waxy crude oil in the pipeline. With increasing mass Grashof number, buoyant forces acting on the fluid become more dominant than viscous forces. This means that the warmer, less dense oil at the pipe centerline rises more readily, while the cooler, denser oil near the oil–gel interface sinks more readily. This circulation of oil promotes a more uniform distribution of heat throughout the pipe, consequently leading to a reduced overall temperature profile.

It is observed in Figure 11 that an increase in the mass Grashof number leads to a decrease in the total concentration profiles of wax molecules in the crude oil pipeline. The observed trend arises from the increased dominance of species buoyancy forces over viscous hydrodynamic forces. With increasing mass Grashof number, buoyant forces propel wax molecules away from the oil–gel interface and toward the bulk of the fluid. The enhanced natural convection currents within the oil promote a more comprehensive mixing of wax molecules across the pipeline’s cross-section. Consequently, the concentration of wax molecules near the pipe centerline diminishes, resulting in a decline in the overall concentration profile.

It is observed in Figure 12 that an increase in the mass Grashof number leads to an increase in the aggregation degree profiles of wax crystals on the pipe wall in crude pipelines. This is because as the mass Grashof number increases, the buoyancy forces acting on the fluid also intensify. These intensified buoyancy forces drive a more vigorous convective flow, which promotes mixing and shear stress within the fluid. This increased turbulence in the flow creates eddies and vortices that promote the dispersion and collision of wax crystals, increasing the likelihood of their aggregation. Aggregation occurs when wax crystals collide and adhere to each other, forming larger clusters. These larger clusters are more likely to deposit on the pipe wall due to their increased inertia and reduced ability to remain suspended in the flow.

It is observed in Figure 13 that an increase in the mass Grashof number causes an increase in the volume fraction occupied by crude oil in the pipeline. The observed trend is because as the mass Grashof number increases, the buoyancy forces acting on the wax crystals become stronger, causing them to disperse more effectively within the crude oil. This dispersion reduces the tendency of wax crystals to accumulate and form deposits on the pipe wall, thereby increasing the volume fraction of crude oil flowing through the pipeline.

4.3. Effects of Varying Eckert Number

It is observed in Figure 14 that an increase in the Eckert number leads to an increase in the temperature profiles of waxy crude oil within the pipeline. This phenomenon is attributed to the enhanced viscous heating effect. The Eckert number represents the ratio of viscous dissipation to thermal conduction. As the Eckert number increases, the viscous shear stresses acting on the oil also increase. This increased shear stress generates more heat through viscous dissipation, which in turn increases the oil temperature. The increased temperature enhances Brownian motion of the wax crystals, causing them to collide more frequently and break down into smaller particles. These smaller particles are less prone to deposition and contribute to a more uniform temperature distribution across the pipeline.

4.4. Effects of Varying Weber Number

It is observed in Figures 15 and 16 that an increase in the Weber number causes a decrease in both radial and axial velocity profiles of waxy crude oil within the pipeline. The Weber number represents the ratio of inertial forces to surface tension forces within the fluid. With an increase in the Weber number, inertial forces, which are responsible for dispersing fluid particles, become relatively more dominant in comparison to surface tension forces, which act to draw fluid particles together. Consequently, the fluid’s resistance to deformation or breakup due to surface tension weakens as inertial forces take precedence. This weakening leads to a decrease in the fluid’s capacity to sustain radial and axial velocities, resulting in a decrease in the velocity profiles.

It is observed in Figure 17 that an increase in the Weber number leads to a decrease in the temperature profiles of waxy crude oil within the pipeline. This decrease is attributed to the enhanced mixing and heat transfer between the oil and the surrounding environment. With an increase in the Weber number, inertial forces, which are responsible for fluid mixing and turbulence, become relatively more dominant compared to surface tension forces, which tend to dampen mixing. This increased dominance of inertial forces leads to enhanced turbulence and mixing within the oil, promoting heat transfer from the oil’s core to its outer regions. Consequently, the temperature at the pipeline’s centerline, which represents the core of the flow, decreases. The increased turbulence also results in the formation of a thinner thermal boundary layer near the pipe wall. The thermal boundary layer is a region of the fluid where temperature gradients are significant. A thinner boundary layer indicates more efficient heat transfer between the oil and the surrounding environment, further contributing to the decrease in temperature of the bulk oil.

It is observed in Figure 18 that an increase in Weber number causes a decrease in the total concentration profiles of waxy crude oil within the pipeline. The observed trend is attributed to enhanced mixing and shear-induced wax dispersion. With an increase in the Weber number, inertial forces, which are responsible for fluid mixing and turbulence, become relatively more dominant compared to surface tension forces, which tend to promote wax particle aggregation and deposition. The increased turbulence promotes the mixing of the oil phases, distributing waxy particles more uniformly throughout the pipeline cross-section. This reduces the concentration of wax particles in the near-wall region, where deposition is most likely to occur.

It is observed in Figure 19 that an increase in the Weber number causes an increase in the aggregation degree profiles of wax crystals within crude oil pipeline. This increase is attributed to the weakening of surface tension forces relative to inertial forces as the Weber number increases. This weakening allows wax crystals to overcome surface tension and collide with each other more frequently. Consequently, wax crystals are more likely to aggregate or stick together and form larger clusters, leading to higher aggregation degree profiles near the oil–gel interface.

It is evident in Figure 20 that an increase in the Weber number causes an increase in the volume fraction occupied by crude oil within the pipeline. The observed trend is because higher Weber numbers represent stronger inertial forces relative to surface tension. As a result, the fluid is less constrained by surface tension and has a greater tendency to flow as a continuous, bulk fluid rather than forming stable droplets. This leads to a higher volume fraction occupied by crude oil within the pipeline because the fluid is less likely to adhere to the walls and is more inclined to occupy a larger portion of the pipeline’s internal volume.

4.5. Skin Friction Coefficient and Rates of Heat and Mass Transfer

The local skin friction coefficient, the local Nusselt number, and the local Sherwood number are computed. The parameters , , , , , and are varied on the local coefficient of skin friction (), local Nusselt number (), and local Sherwood number (), and their numerical values are presented in Table 1.

From the table, the following observations are noted:(i)An increase in the Reynolds number () increases the skin friction coefficient but leads to a decrease in the Nusselt number and Sherwood number. This phenomenon is a consequence of the relationship between wall shear stress and the velocity gradient. As rises, the velocity profile increases, resulting in higher wall shear stress and, consequently, an elevated skin friction coefficient. Additionally, the increase in causes the thermal boundary layer to thicken, leading to a decreased Nusselt number. The Sherwood number, on the other hand, diminishes because higher Reynolds numbers cause a thickening of the concentration boundary layer, reducing the rate of species transport within it.(ii)An increase in the thermal Grashof number () increases the skin friction coefficient and Sherwood number but has no effect on the Nusselt number. The rise in the skin friction coefficient is attributed to the thermal Grashof number causing an increase in velocity, which subsequently thins the velocity boundary layer. Conversely, the elevation in the Sherwood number is a consequence of the thermal Grashof number leading to a reduction in the thickness of the concentration boundary layer, thereby increasing the rate of species transport within this layer. The Nusselt number, however, remains unaffected.(iii)An increase in the mass Grashof number () results in an increase in the skin friction coefficient, Nusselt number, and Sherwood number. This observed pattern is because the higher values of the mass Grashof number lead to the thinning of the velocity, thermal, and concentration boundary layers. This, in turn, results in a higher rate of transportation within these boundary layers, accounting for the increased values of the skin friction coefficient, Nusselt number, and Sherwood number.(iv)An increase in the Eckert number () causes an increase in the skin friction coefficient and Sherwood number but decreases the Nusselt number. This pattern arises from the fact that higher values of the Eckert number result in an increased fluid velocity, subsequently leading to higher wall shear stress. As the Eckert number rises, it also contributes to the thickening of the thermal boundary layer, which reduces the heat transfer rate at the pipeline wall, consequently lowering the Nusselt number. However, this thickening of the thermal boundary layer enhances the rate of species transport, leading to an increase in the Sherwood number.(v)An increase in the Schmidt number () causes a decrease in the skin friction coefficient and Nusselt number but increases the Sherwood number. This phenomenon is attributed to the influence of Schmidt number on fluid behavior. As increases, it decelerates the axial velocity of fluid particles, causing the velocity boundary layer to thicken, which in turn reduces the motion of fluid particles and leads to a decrease in the skin friction coefficient. Additionally, the thermal boundary layer thickness increases with higher , resulting in reduced heat transfer at the pipeline wall and a lower Nusselt number. In contrast, the concentration boundary layer thickness decreases with increasing , leading to an enhanced rate of species transportation and an increase in the Sherwood number.(vi)An increase in the Weber number () decreases the skin friction coefficient but it increases the Nusselt number and Sherwood number. This trend can be attributed to the influence of the Weber number on the behavior of the fluid. An increase in results in the deceleration of fluid particle velocities, leading to the thickening of the velocity boundary layer, which subsequently reduces particle motion and causes a decrease in the skin friction coefficient. Moreover, an increase in causes the thermal boundary layer thickness to decrease, resulting in a higher rate of heat transfer. The Sherwood number experiences an increase because higher values of the Weber number lead to a reduction in concentration boundary layer thickness, enhancing the rate of species transport.

5. Validation

The findings from this research are validated against experimental data from Ying et al. [14]. In particular, the fluid temperature profile is compared in both studies, as shown in Figures 21 and 22. It is observed that the temperature profiles follow a similar trend as time increases. This validation confirms the model’s accuracy in predicting wax deposition under laminar flow conditions.

6. Summary and Conclusions

This study presents a novel numerical investigation of wax deposition from multiphase flow in field-scale crude oil pipelines. Waxy crude oil and the solid wax deposit are treated as two immiscible phases separated by a smooth, continuous interface. Two deposition mechanisms, such as molecular diffusion and shear dispersion, are considered. The model equations, in the form of coupled nonlinear PDEs governing the flow, are discretized in time using a second-order semi-implicit finite difference scheme and in space using the BSCM. The study analyzes the impact of varying the flow parameters on the flow variables. The following conclusions are drawn from this study:(i)The aggregation of wax crystals in the pipeline decreases by at most 2.5% with increasing Reynolds number from 2.2361 to 3.1361. However, it increases by at most 3.4% with increasing mass Grashof number from 5 to 11. Additionally, it increases by at most 4.8% with increasing Weber number from 1.0 to 2.5.(ii)The total concentration of waxy components in the pipeline increases by at most 10.0% with increasing Reynolds number from 2.2361 to 3.1361. However, it decreases by at most 8.0% with increasing mass Grashof number from 5 to 11. Additionally, it decreases by at most 20% with increasing Weber number from 1.0 to 2.5.(iii)The fraction of the volume occupied by waxy crude oil in the pipeline decreases by up to 10.0% with increasing Reynolds number from 2.2361 to 3.1361. However, it increases by up to 9.0% with increasing mass Grashof number from 5 to 11, and by up to 25% with increasing Weber number from 1.0 to 2.5.(iv)The temperature of waxy crude oil in the pipeline increases by at most 0.001% with increasing Reynolds number from 2.2361 to 3.1361, and by at most 0.002% with increasing Eckert number from 1.2 to 1.5. However, it decreases by at most 0.0034% with increasing mass Grashof number from 5 to 11. Additionally, it decreases by at most 1.3% with increasing Weber number from 1.0 to 2.5.(v)The radial velocity of waxy crude oil in the pipeline increases by at most 1.0% with increasing Reynolds number from 2.2361 to 3.1361. However, it decreases by at most 1.8% with increasing mass Grashof number from 5 to 11. Additionally, it decreases by at most 1.3% with increasing Weber number from 1.0 to 2.5.(vi)The axial/streamwise velocity of waxy crude oil in the pipeline decreases by at most 0.7% with increasing Reynolds number from 2.2361 to 3.1361 and by at most 1.2% with increasing Weber number from 1.0 to 2.5. However, it increases by at most 3.0% with increasing mass Grashof number from 5 to 11.(vii)Skin friction coefficient increases from 0.1230 to 0.4022 with increasing Reynolds number from 2.2361 to 5.2361 and from 0.1230 to 0.1277 with increasing mass Grashof number from 5 to 20. It also increases from 0.1230 to 0.1334 with increasing Weber number from 1.0 to 2.5.(viii)The Nusselt number increases from 1.9907 to 4.9834 with increasing Reynolds number from 2.2361 to 5.2361 and from 1.9907 to 2.0225 with increasing mass Grashof number from 5 to 20. It also increases from 1.9907 to 2.0434 with increasing Weber number from 1.0 to 2.5.(ix)The Sherwood number increases from 1.3916 to 7.2234 with increasing mass Grashof number from 5 to 20 and from 1.3916 to 1.6002 with increasing Weber number from 1.0 to 2.5.

The insights gained from this study provide valuable guidance for optimizing pipeline operations, designing effective wax control strategies, and enhancing pipeline integrity management in field-scale crude oil transportation systems. Pipeline operators can utilize the model to identify critical flow parameters influencing wax deposition and optimize these parameters to minimize wax accumulation. While the current study focuses on laminar flow, future research should extend this study to incorporate the effects of turbulence and droplet interactions, as these factors could also play a vital role in wax deposition.

Nomenclature

Symbols

Cp:Specific heat at constant pressure, (J/(kg·K))
Cinterface:Wax concentration at oil–gel interface, (kg/m3)
Dp:Shear dispersion coefficient, (m)
Dd:Molecular diffusion coefficient of wax, (m2/s)
:Dimensionless diameter of water droplet
g:Gravitational acceleration, (m/s2)
k:Thermal conductivity, (W/m·K)
K1:Smooth positive function of fluid temperature, (1/s)
K2:Smooth positive function of fluid temperature, (1/Pa)
L:A sufficiently large finite number to approximate the asymptotic behavior at infinity, (dimensionless)
P:Pressure, (dimensionless)
Q:Volume flow rate, (m3/s)
R:Radius of clean pipe, (m)
Reff:Effective radius for oil flow, (m)
Tinterface:Fluid temperature at oil–gel interface, (K)
(, ):Dimensionless cylindrical coordinate variables
:Dimensionless time
(u,v):Dimensionless velocity components
x:Weight fraction of wax crystals in the gel layer
Da:Darcy number
Ec:Eckert number
GrT:Thermal Grashof number
GrC:Mass Grashof number
Pr:Prandtl number
Pe:Peclet number
Re:Reynolds number
St:Stanton number
Sc:Schmidt number
We:Weber number
Cf:Skin friction coefficient
Nuz:Nusselt number
Shz:Sherwood number.

Greek Symbols

αavg:Average aspect ratio of the wax crystals
αm:Aggregation degree of wax, (dimensionless)
βC:Concentration volume expansion coefficient, (m3/kg)
βT:Thermal volume expansion coefficient, (1/K)
:Dimensionless deposit thickness
μ:Coefficient of dynamic viscosity, (Ns/m2)
ϕ:Dimensionless total concentration of wax
ϕi:Proportion of volume occupied by the th phase
ρ:Fluid density, (kg/m3)
σ:Surface tension coefficient, (N/m)
Θ:Dimensionless temperature of waxy crude oil
ϕ:Angle of elevation of pipe from the horizontal.

Subscripts

f:Fluid phase
i:(Oil, water, gel)
water:Water droplets
oil:Crude oil
gel:Deposit or gel-like layer
mix:Mixture fluid
wall:Condition at the pipe inner wall
interface:Condition at the oil–gel interface
∞:Condition at the free stream.

Data Availability

The findings of the present study are supported by data (in the form of MATLAB® source code) which can be availed by the corresponding author upon reasonable request. The source code was generated by MATLAB® software version 9.10.0.1602886 (R2021a).

Conflicts of Interest

The authors wish to declare no conflicts of interest regarding the publication of this research article.

Acknowledgments

Francis Oketch Ochieng is grateful to his supervisors, Prof. Mathew N. Kinyanjui, Prof. Jeconia O. Abonyo, and Dr. Phineas R. Kiogora, for their insightful comments and guidance throughout the development of this research article.