Three-Dimensional CFD Modeling on the Thermal Characteristics of Buried Oil Pipeline Involving the Heat Transfer of Wax Layer

Abstract

It is common to have a wax layer on the inner wall of waxy crude oil pipeline. However, the study of the wax layer on the heat transfer of buried pipeline systems is inadequate due to its instability of composition and properties; it may lead to the inaccurate prediction of the pipeline temperature field. Based on the finite element simulation technology, a three-dimensional heat transfer model of buried crude oil pipeline involving wax layer was proposed and solved. The thermal effect of the wax layer on pipeline system was analyzed quantitatively. Numerical results show that the average deviation of soil temperature near the pipeline reach 1.42 K when there is a 4 mm wax layer. Among different thermal conductivity models of heterogeneous materials, the EMT model plays best in predicting the conductivity of a waxy layer. By setting different working conditions, the influence mechanisms of several thermal influencing factors are discussed. The results show that the thermal influence range of heated pipe is positively associated with oil temperature and velocity. The core thermal response zone is about 12 m along the X-axis. Beyond 8 m depth from ground surface, the temperature fluctuation of soil is almost unaffected by the oil pipeline.

1. Introduction

Crude oil still plays an important role in today’s energy structure, and the amount of crude oil production and throughput is considerably large. Waxy crude oil is widely distributed in the oil fields worldwide. Due to its high viscosity and high thermal sensitivity [1], waxy crude oil is usually heated during transportation to improve flow performance [2], and is usually buried by soil between oil stations. The large temperature variation of waxy crude oil, and a plenty of influencing factors present challenges to the flow safety and accurate simulation of the transportation process. Simulation and prediction of the buried oil pipeline are of great significance to the safety, efficiency, and energy-saving of crude oil transportation.

In the field of buried waxy crude oil pipeline, systematic research begins with the establishment of a mathematical model and the solution of a numerical model. Barletta et al. [3] investigated the steady-periodic heat transfer between offshore buried pipelines for the transport of hydrocarbons and their environment. The unsteady two-dimensional conduction problem was written in a dimensionless form, and it was transformed into a steady two-dimensional problem and solved numerically by means of the finite-element software package. Yu et al. [4] established a physical model of a buried oil pipeline under normal operation, and an approach combining unstructured-finite-volume and finite difference methods was applied to solve the governing equations, but the fluid is treated as a boundary condition rather than material when discrete computing domain with unstructured mesh. Yu et al. [5] then proposed the pod-Galerkin reduced-order model to realize the efficient solution of the pipeline heat transfer model. Zhou et al. [6] further improved the solution method by adopting PSO-DE algorithm. Scholars have made progress in thermal characteristics analysis of buried crude oil pipelines. Xing et al. [7] proposed a mathematical model for the preheating process of waxy crude oil pipelines, solved the governing equations by combining the finite volume method with the finite difference method, and discussed factors such as the medium, the plan, the temperature, and the flow rate for preheating. In the last decade, studies with the cooling process which combine the non-Newtonian properties of waxy crude oil have been carried out. Cheng et al. [8] established the temperature drop model of the oil pipeline after shutdown, and emphatically analyzed the movement of the phase transition interface of wax deposition from pipe wall to pipe center. Z. Dai [9] used an enthalpy–porous medium method to establish the thermal model of crude oil in shutdown, and concentrated on the solidification law and allowable shutdown time of crude oil. H. Dong et al. [10] studied the thermal system coupling by the crude oil, pipeline, and soil during the pipeline restart process, emphatically analyzed the influence of thermal characteristics of pipeline restart under different working conditions, and divided the pipeline into three areas with different dominant mechanisms after restart. Zhao et al. [11] used the momentum source term method and additional specific heat capacity method to investigate the thermal characteristics of wax crude oil pipeline during its shutdown and studied the quantitative effects of crude oil viscosity, specific heat capacity, pipeline diameter, and insulating thickness on the thermal characteristics of waxy crude oil pipeline during the shutdown. With the improvement of computer performance and the popularization of commercial simulation software, the research on CFD simulation of buried crude oil pipeline has also made progress. Wang et al. [12] established a two-dimensional steady-state model of buried oil pipeline, and used the software Fluent 6.0 to analyze the morphology of soil temperature field, but has not yet gelled into a satisfying result. Desamala et al. [13] investigated the transition boundaries of different flow patterns for moderately viscous oil–water two-phase flow through a horizontal pipeline, and adopted Volume of Fluid approach including effect of surface tension to predict the flow pattern. However, only the 2d-simulation of the organic glass tube with an internal diameter of 0.025 m was carried out. Combined statistical experimental design with a developed computational fluid dynamics model, Rukthong [14] explored the heat conduction and the thermal convection behavior in the thin-wall pipeline with laminar unsteady state flow, and located an origin point where wax precipitate in the pipeline by using response surface methodology. However, in this model, there was no other heat transfer material existed outside the pipe and was simply set as the ambient temperature. Ji [15] simulated the temperature drop processes at stable operation and shutdown of buried crude oil pipeline by software Fluent 16.0, and then proposed the method of dividing the temperature-drop stage during the shutdown of the pipeline by the relation of the storage modulus, loss modulus, and loss angle with temperature. However, the heat transfer material in the model is not fully considered and the wax layer is neglected.

Studies on the physical properties and deposition laws of the wax layer have been carried out. Some scholars describe the wax deposits by experiment. Most of these experiments used instruments to measure pipe samples with wax deposits. Isaksen et al. [16] used ultrasonic positioning technology to measure the distribution of wax deposits online in the testing pipeline. Rainer et al. [17] used laser measuring devices to measure the thickness and roughness of the wax deposits of single-phase crude oil under different temperatures and flow conditions. The approach of experimental testing is unrepresentative, because the quality of different crude oil varies greatly, and the operating state of a certain crude oil pipeline is not constant. Using mathematical models to describe the wax deposits is more popular, including thermodynamic models and kinetic models. Most thermodynamic models adopt the method of phase equilibrium or phase transformation to predict the wax deposition point [18] and amount of sediment [19]. In this field, kinetics models of wax deposition have been established including the Burger model [20], the Hsu model [21], the Hernandez model [22], and the Huang model [23]. However, these models were almost focused on sediment thickness or effective volume in the pipe. In recent years, some prediction models based on intelligent algorithms were used to describe the wax deposition prediction [24]. These models are essentially based on the statistics of wax deposits in particular oil pipelines. Cheng et al. [25] and G. Yu et al. [26] used empirical values as the thermal conductivity of wax layer in the calculation of crude oil pipeline. Jin et al. [27] studied the effective thermal conductivity of the wax layer by using the assumption of dualistic heterogeneous mixture hypothesis, and pointed out that the thermal conductivity of wax deposits is highly positively correlated with the wax mass fraction. However, the influence of wax layer on the thermal calculation of buried crude oil pipeline transportation is not clear in the existing studies.

To sum up, the current studies on the simulation of buried waxy crude oil pipeline have the following shortcomings. Firstly, the existing methods lack the complete three-dimensional simulation of buried pipeline system coupling heat and mass transfer. Most studies, especially those that take analytical solutions, only build two-dimensional models, and decoupled the true three-dimensional system into two parts: the oil flow and the thermal system cross-section. Secondly, there is lack of quantitative influence of wax layer on heat transfer of hot oil pipeline. In many cases, the wax layer has been ignored or the thermal conductivity of the wax layer has not been rigorously demonstrated. These simplifications will affect the thermal calculation results and lead to inaccurate temperature field prediction.

The purpose of this work is to realize the three-dimensional simulation of buried crude oil pipeline involving the heat transfer of wax layer. In the process of simulation, we quantify the unclear influencing factors such as wax layer to optimize the accuracy of the simulation as much as possible. A three-dimensional coupled heat transfer model of oil flow-wax layer-steel pipe-anticorrosive coating-insulating layer-soil was proposed. The necessity and quantitative deviation of wax layer on heat transfer of buried crude oil pipeline system were proved. Based on the thermal conductivity model of heterogeneous materials, the effective thermal conductivity of wax layer was predicted, and the optimal model was evaluated and recommended. The response law of thermal characteristics of pipeline system changing with thermal influencing factors was revealed.

The structure of this paper is organized as follows. Firstly, the background and pipeline data of the study were stated. Then, the mathematical model, the governing equations, the numerical scheme and boundary condition of the coupled heat transfer pipeline system were described. Finally, the simulation and discussion of the waxy crude oil transportation process under different working conditions were carried out.

2. Background

2.1. Characteristics of the Pipeline

In this work, we investigated the pipeline section of the Huage pipeline in Qinghai oilfield, China, which transported waxy crude oil (QH oil) with solidifying point of 305 K. The total length of the pipeline is 39,672 m; and the outer diameter of the steel pipe is 355.6 mm. The buried depth of the pipeline is 1.3 m on average.

The external anticorrosive coating is applied widely to the protection on the outside of the pipes in soil. To reduce heat dissipation, the oil pipeline is covered by the insulating layer such as polyurethane foam. The properties of all the materials used in the modeling were listed in

Table 1. Properties of the materials used in the modeling.

	Parameters	Thermal Conductivity W/(m·K)	Density kg/m3	Specific Heat J/(kg·K)
Material
Steel pipe		48	7850	483
Anticorrosive coating		0.17	1050	1680
Insulating layer		0.032	48	1900
Wax layer		0.1804	830	1900
Soil		1.2	1400	2100
Air		0.0242	1.225	1006.4

, including the thermal conductivities, the density, and the specific heat.

Figure 1 shows the geometry of buried oil pipeline system in this work. The entire system is a rectangular prism containing oil fluid, several circular materials, and the soil. The cross-section of the system is called the thermal affected zone [28]. According to the thermal test results from existing reports [29], a sufficiently large heat affected zone was adopted in this system to capture thermal characteristics effectively. The length of $H_d$ and $L$ are both 10 m.

2.2. Description of the Crude Oil

The QH oil has the characteristic of high wax content. Wax precipitation test showed that the mass fraction of wax deposits was up to 18%. According to the experiment, the physical property parameters of QH oil change obviously with temperature. The density, thermal conductivity, and volumetric expansion coefficient of crude oil can be calculated as follows:

$$\rho =893.4-0.193\times (T-273.15)$$

(1)

$$\lambda =0.154-0.0000829\times (T-273.15)$$

(2)

$$\beta =0.8255/(1118.65-0.8255\times (T-273.15))$$

(3)

where $\rho$ is the density of crude oil at different temperatures, kg/m³; $\lambda$ is the thermal conductivity of crude oil, W/(m∙K), $\beta$ is the volume expansion coefficient of crude oil.

The viscosity—temperature correlation equation was fitted by following formula.

$$\lg \eta =\{\begin{array}{l}2.035-0.021\times (T-273.15),311K<T\le 316K\\ 1.597-0.011\times (T-273.15),316K<T\le 343K\end{array}$$

(4)

where $\eta$ is the dynamic viscosity of crude oil, mPa∙s; $T$ is the temperature of crude oil, K.

2.3. Problem Definitions and Conditions

In order to improve the accuracy of the simulation, the parameters of the materials used in the model should be as close to the reality as possible, such as anticorrosive coating, thermal insulating layer, and wax layer. Among them, the parameters of wax layer is not clear. It is because of the difficulties in observation since the wax layer is located on the inner wall of pipe and the specificity of composition in different kinds of crude oil. The thermal conductivity of the wax layer could be predicted based on crude oil’s properties and state. The temperature and velocity of crude oil at the inlet of the pipeline were selected as the thermal influencing factors of the fluid. The influence of the external environment on the heat transfer of the buried pipeline system should be based on the climatic conditions where the pipeline is located. The climate along the investigated pipeline varies significantly in different seasons, with an average annual temperature of 283 K and average wind speed of 21 m/s. The air temperature and the convective heat transfer intensity in different seasons were selected as the thermal influencing factors of the external environment.

All these analyses were to understand the thermal characteristics of the buried waxy crude oil pipeline through numerical simulation coupling flow and heat transfer. The response law of thermal characteristics of pipeline system with the change of thermal influencing factors was also desired. To ensure the accuracy of the simulation, not only the materials involved in heat transfer and the corresponding thermal conductivities must be selected and calculated precisely, but also the numerical schemes and boundary conditions must be set properly.

3. Mathematical Model

3.1. Assumption and Governing Equation

In the process of establishing and solving the model, some reasonable assumptions were adopted. We consider the oil flow pattern to be fully developed when it flows through pipelines because the length of the pipe in this work is 200 m and is large enough. In reality, the soil around the pipe is a mixture of complex components. To define the thermal parameters of soil, it is assumed that the soil around the pipeline is uniform and isotropic [30]. Studies show that the thickness of the wax layer on the inner wall of the crude oil pipeline is unevenly distributed [31]. To facilitate the establishment of the model, the wax layer on the inner wall is assumed to be of the same thickness along the entire pipe.

The oil flow in the pipeline satisfies the conservation equations of mass, momentum, and energy. Heat conduction and natural convection are considered in three dimensions. Phase transition is not considered because the operating temperature of crude oil is much higher than the wax precipitation point. The governing equations of the fluid in the pipe are as follows.

(1) Mass conservation equation of flow

$$\frac{\partial {\rho }_0}{\partial \tau }+\frac{\partial }{\partial x}({\rho }_{0}u)+\frac{\partial }{\partial y}({\rho }_{0}v)+\frac{\partial }{\partial z}({\rho }_{0}w)=0$$

(5)

(2) Momentum conservation equation of flow

$$\frac{\partial u}{\partial \tau }({\rho }_{0}u)+\frac{\partial u}{\partial x}({\rho }_{0}uv)=-\frac{\partial p}{\partial x}+{\rho }_0{g}_{\mathrm{x}}+F_{\mathrm{x}}$$

(6)

$$\frac{\partial v}{\partial \tau }({\rho }_{0}v)+\frac{\partial }{\partial y}({\rho }_{0}vw)=-\frac{\partial p}{\partial y}+{\rho }_0{g}_{\mathrm{y}}+F_{\mathrm{y}}$$

(7)

$$\frac{\partial w}{\partial \tau }({\rho }_{0}w)+\frac{\partial }{\partial z}({\rho }_{0}wx)=-\frac{\partial p}{\partial z}+{\rho }_0{g}_{\mathrm{z}}+F_{\mathrm{z}}$$

(8)

where $F_{\mathrm{x}}$, $F_{\mathrm{y}}$, $F_{\mathrm{z}}$ are the external component of volume force in x, y, z direction, N/kg; $g_{\mathrm{x}}$, $g_{\mathrm{y}}$, $g_{\mathrm{z}}$ are the component of gravity volume force in x, y, z direction, N/kg.

(3) Energy conservation equation of flow

$$\frac{\partial t}{\partial \tau }+\frac{\partial }{\partial x}({\rho }_{0}ut)+\frac{\partial }{\partial y}({\rho }_{0}vt)+\frac{\partial }{\partial z}({\rho }_{0}wt)=\frac{{\lambda }_0}{c_{\mathrm{p}}}\left(\frac{{\partial }^{2}t}{\partial x^2}+\frac{{\partial }^{2}t}{\partial y^2}+\frac{{\partial }^{2}t}{\partial z^2}\right)$$

(9)

where ${\lambda }_0$ is the thermal conductivity of crude oil, W/(m∙K); $c_{\mathrm{p}}$ is the specific heat capacity of crude oil, J/(kg·K); $t$ is the temperature of the liquid, K.

In this buried crude oil pipeline system, mathematical models of heat transfer between various solid materials, including wax layer, steel pipe, anticorrosive coating, insulating layer and soil, were established.

(1) Heat transfer equation for each layer of pipeline

$${\rho }_n{c}_n\frac{\partial T_n}{\partial \tau }=\frac{\partial }{\partial x}({\lambda }_n\frac{\partial T_n}{\partial x})+\frac{\partial }{\partial y}({\lambda }_n\frac{\partial T_n}{\partial y})+\frac{\partial }{\partial z}({\lambda }_n\frac{\partial T_n}{\partial z})$$

(10)

where ${\lambda }_n$ is the thermal conductivity of the number n pipe layer material, W/(m∙K); ${\rho }_n$ is the density of the number n pipe layer material, kg/m³; $c_n$ is the heat capacity of the number n pipe layer material, J/(kg∙K); $T_n$ is the temperature of the number n pipe layer, K; the number n (n = 1, 2, 3, 4) denote, respectively, wax deposit, pipeline wall, anticorrosive coating, and insulating layer.

(2) Heat transfer equation for soil

$${\rho }_{\mathrm{S}}{c}_{\mathrm{S}}\frac{\partial T_{\mathrm{S}}}{\partial \tau }=\frac{\partial }{\partial x}({\lambda }_{\mathrm{S}}\frac{\partial T_{\mathrm{S}}}{\partial x})+\frac{\partial }{\partial y}({\lambda }_{\mathrm{S}}\frac{\partial T_{\mathrm{S}}}{\partial y})+\frac{\partial }{\partial z}({\lambda }_{\mathrm{S}}\frac{\partial T_{\mathrm{S}}}{\partial z})$$

(11)

where ${\lambda }_{\mathrm{S}}$ is the thermal conductivity of soil, W/(m∙K); ${\rho }_{\mathrm{S}}$ is the density of the soil, kg/m³; $c_{\mathrm{S}}$ is the heat capacity of the soil, J/(kg∙K); $T_{\mathrm{S}}$ is the temperature of the soil outside the pipe, K.

3.2. Boundary Conditions

According to the characteristics of the buried crude oil pipeline system, the boundary conditions of the numerical models are defined as follows.

(1) Boundary conditions of the solid domain

The upper boundary of the model is the interface between soil and air, where the convective heat transfer coefficient and fluid temperature are specified.

$${\frac{\partial T_{\mathrm{S}}}{\partial y}|}_{y=+H_1}=\frac{a_{\mathrm{a}}}{{\lambda }_{\mathrm{S}}}(T_{\mathrm{S}}-T_{\mathrm{a}})$$

(12)

where $T_{\mathrm{S}}$ is the temperature of the soil outside the pipe, K; $T_{\mathrm{a}}$ is the temperature of the air, K; ${\lambda }_{\mathrm{S}}$ is the thermal conductivity of soil, W/(m∙K); $a_{\mathrm{a}}$ is the heat transfer coefficient of air to the surface of the soil, W/(m² K), $H_1$ is the distance from the origin of the coordinates to the ground surface, m.

The lower boundary of the model is the soil thermostatic layer, and belongs to the first boundary condition, where the temperature values are specified.

$${T_{\mathrm{S}}|}_{y=-H_2}=const$$

(13)

where $H_2$ is the distance from the origin of the coordinates to the soil thermostatic layer, m.

The left and right boundaries of the model are the limits of the heat affected zone, where the heat flux is specified to be zero.

$${\frac{\partial T_{\mathrm{S}}}{\partial x}|}_{x=-L}=0$$

(14)

$${\frac{\partial T_{\mathrm{S}}}{\partial x}|}_{x=+L}=0$$

(15)

where $L$ is the half-width of the heat affected zone.

(2) Boundary conditions of the fluid domain

Since the inlet of the pipe is set as the velocity inlet and a given temperature value, the boundary conditions are constrained by the following equation.

The velocity inlet boundary is used at the pipe inlet, where the velocity and temperature are given.

$${V|}_{z=0}=V_0$$

(16)

$${T|}_{z=0}=T_0$$

(17)

The pressure outlet boundary is used at the pipe outlet, where the pressure is given.

$${P|}_{z=z\max }=P_0$$

(18)

where $V_0$, $T_0$, $P_0$ are the given initial values.

3.3. Numerical Scheme

Commercial software Fluent, ANSYS Corporation, Pittsburgh, PA, USA, was used to solve the model in this work, and the software version was 19.0. Based on the pressure base solver, the PISO algorithm was used to solve the pressure–velocity coupling. According to the Reynolds number, the computation of turbulent flow in the oil pipe was studied. The standard k-ε model was used to solve the turbulent oil pipe flow, due to its computational speed and stability. The governing equations, species transport equation, k equation, and ε equation were all discretized by a second-order upwind scheme. The PRESTO method was adopted as pressure interpolation mode. The viscosity, density, and thermal conductivity of crude oil were defined using the DEFINE_PROPERTY macro, which was loaded as an interpreted UDF.

All soil boundaries in the computational domain were set as Wall boundary conditions. Taking into account the natural convection patterns between air and soil on the ground, the convective heat transfer coefficient at the upper boundary was set. The lower boundary was the constant temperature boundary because of the soil’s thermostatic layer. The boundary on both sides was adiabatic. The inlet boundary of the pipeline was the velocity boundary and the outlet was the pressure boundary. The calculation is performed on a GenTai HPC Tower with Professional Version win10, Intel Xeon (R) CPU Silver 4114 @ 2.20 GHz, 20 processors, and 96 GB RAM. The convergent criteria for all calculations are set as the residual in the control volume for each equation is smaller than 10⁻⁶.

3.4. Computational Domain and Mesh

The entire computational domain is a 20 m × 10 m × 200 m rectangular prism. In the X, Y, and Z directions, the coordinate ranges are (−10 m, 10 m), (−4.5097 m, 5.492 m), and (0 m, 200 m), respectively. In the system, the thickness of wax layer, anticorrosive coating and insulating layer is 4 mm, 6 mm, and 40 mm, respectively.

Figure 2 shows the mesh of the simulation model in this work. ICEM software was used to generate hexahedral structured mesh, and O-type mesh was generated for the fluid domain. To simulate the temperature field more accurately, considering the inhomogeneity of temperature gradient distribution, a dense grid near the pipe and a sparse grid far from the pipe were adopted. Ensure the accuracy of computing, using the local mesh refinement approach at the region of pipe wall, anticorrosive coating, and insulating layer. Moreover, to extract the calculation results, several cross-sections along the axial direction of the pipeline were established in the computational domain.

To demonstrate the heat transfer law between flow and solid materials further, several auxiliary paths were set on the cross-section with axial coordinate 100. The details of all auxiliary paths were shown in

Table 2. Coordinates and physical meaning of auxiliary paths.

Number	Start Point	Endpoint	Length	Physical Meaning
Path 1	(0, −4.5, 100)	(0, 5.5, 100)	222 m	The centerline of the pipeline cross-section
Path 2	(−10, 5, 100)	(10, 5, 100)	20 m	The depth of 0.5 m below the ground surface
Path 3	(−10, 1.5, 100)	(10, 1.5, 100)	20 m	The depth of 4 m below the ground surface
Path 4	(−10, −0.5, 100)	(10, −0.5, 100)	20 m	The depth of 6 m below the ground surface
Path 5	(−10, −1.5, 100)	(10, −1.5, 100)	20 m	The depth of 7 m below the ground surface
Path 6	(−10, −2.5, 100)	(−10, −2.5, 100)	20 m	The depth of 8 m below the ground surface
Path 7	(−10, −3.5, 100)	(10, −3.5, 100)	20 m	The depth of 9 m below the ground surface
Path 8	(−10, 4.19, 100)	(−0.4, 4.19, 100)	9.6 m	Test line for grid independent verification

. These auxiliary paths came in two types, one type was perpendicular to the X-axis, and the other was perpendicular to the Y-axis. Path1 is the centerline of the cross-section at z = 100 m. It passes through coordinates (0 m, 0 m), and belongs to the first type. The paths of the second type are the horizontal line at different depths on the cross-section, including Path 2 to Path 7.

Mesh independence verification was carried out by computations of four groups of grids with different numbers, and the results shown in Figure 3. The numerical results of soil temperature along Path 8 were taken as the verification benchmark. When the mesh number reached 1000 thousand, the numerical results along the test line had no distinguished change. Therefore, considering both accuracy and efficiency of calculation, the grid with 1,003,968 meshes and 1,101,190 nodes was used for numerical calculation.

4. Numerical Results and Discussion

All simulations were carried out on FLUENT 19.0. These models successfully predicted the flow in the pipe and the temperature field in the whole solution domain.

4.1. Thermal Conductivity of the Wax Layer

Waxy crude oil pipelines usually form the wax layer on the inner wall after running for a long time. The wax layer is inhomogeneous and is easily affected by the state of the crude oil in transit. Therefore, adopting empirical values in calculation may produce deviation. To calculate the effective heat transfer coefficient of the wax layer accurately, the thermal conductivity models of heterogeneous materials [32,33] were adopted, including the Series model, the Parallel model, the ME model, and the EMT model. All these models were described in

Table 3. Description of the effective thermal conductivity models.

Name	Expression	Fundamental Assumption
Series model	$k_{\mathrm{e}}={\left(\frac{v_1}{k_1}+\frac{v_2}{k_2}\right)}^{-1}$	Assuming that the two components are distributed in different horizontal layers when the heat flow is in the vertical direction.
Parallel model	$k_{\mathrm{e}}=v_1{k}_1+v_2{k}_2$	Assuming that the two components are distributed in different vertical layers when the heat flow is in the vertical direction.
ME model	$k_1<k_2,k_{\mathrm{e}}=\frac{k_2{v}_2+k_1{v}_1\frac{3k_2}{2k_2+k_2}}{v_2+v_1\frac{3k_2}{2k_2+k_1}}$	Assuming that one substance is distributed in different groups as small spheres which are far apart. The thermal conductivity of the continuous phase in the mixture is lower than that of the dispersed phase.
EMT model	$v_1\frac{k_1-k_{\mathrm{e}}}{k_1+2k_{\mathrm{e}}}+v_2\frac{k_2-k_{\mathrm{e}}}{k_2+2k_{\mathrm{e}}}=0$	The composition of the mixture is assumed to be completely random, and there is no distinction between continuous and dispersed phases.

.

In the equation above, $k_{\mathrm{e}}$ is the effective thermal conductivity of the mixture, W/(m∙K); $k_1$ and $k_2$ are the thermal conductivity of the two components, W/(m∙K); $v_1$ and $v_2$ are the volume fractions of the two components, respectively, %. Based on the known thermal conductivity and volume fraction of each component, the effective thermal conductivity of the mixture can be calculated.

The wax layer is the compound of wax crystals and liquid oil. The amount of wax crystals in the wax layer mixture is related to thermal history, shear rate, and composition of the crude oil [34]. According to the physical properties of QH oil, the mass fraction of wax crystal was set to 50%. At 323 K, the thermal conductivity of wax crystal is 0.25 W/(m∙K), and that of QH oil is 0.125 W/(m∙K). The effective thermal conductivity of the wax layer calculated by the four methods above is 0.1667, 0.1875, 0.1786, and 0.1804 W/(m∙K), respectively. Figure 4 shows the effective thermal conductivity of the wax layer obtained by the above calculation method.

As can be seen from the Figure 4, the effective thermal conductivity increases with the mass fraction of wax crystal, and the calculation result of the Parallel method is the largest while the Series method is the smallest, and the calculation result of the ME and the EMT method is moderate. When the mass fraction of wax crystal is in the range of 30% to 65%, the effective thermal conductivity of the wax layer is in the range of 0.155–0.195 W/(m·K).

In this work, combined with the calculation method of effective thermal conductivity above, we compared the simulation results with the experimental measurements. At the pipeline site, five temperature detecting sensors are set at the starting point (A to E). At the end of the pipeline, a similar location of sensors was adopted (F to J). The installation positions of the sensor are shown in Figure 5. The distance from point A to the ground surface is 1.0 m.

The working conditions at field were as follows: the oil velocity was 0.7842 m/s, the ambient temperature was 263.15 K, and the wax layer thickness was 4 mm, the actual oil temperature at inlet and outlet of the pipe were 323.15 K and 314.85 K, respectively. It is worth noting that the total length of the pipe in the field is 39,672 m, while the length of the numerical model is 200 m, thus the outlet section in the model cannot represent the real outlet of the field pipe. To solve this problem, in the numerical cases, the actual inlet oil condition (323.15 K oil temperature) was used at the inlet of the model to obtain the temperature field in the plane where points A–E located. We then used the actual outlet oil condition (314.85 K oil temperature) at the inlet of the model to obtain the temperature field in the plane where points F–J located. In this way, the deviations at all detecting points could be obtained to verify the accuracy.

By using different effective thermal conductivity models, four groups of thermal conductivity were calculated and were adopted in the numerical cases. We made an auxiliary line from point C to point E and derive all the temperature values along the auxiliary line in the four groups of cases. The variation curves of temperature from point C to E with different wax layer thermal conductivity are shown in Figure 6. It can be seen that the calculated temperature values show a trend of gradual decrease, with minor deviations of the four models. The absolute deviation of the calculated values at each position from point C to E is less than 1 K. The temperature distributions of the curves using the ME method and EMT method are the closest to the measured value, indicating that the ME method and the EMT method are more suitable for estimating the effective thermal conductivity of waxy layer. Although the two models have similar computational accuracy, the EMT model is better because it has no restrictions on the continuous phase and dispersed phase, and the assumption of the EMT model is more suitable for the wax layer with flaky or needle-like structure. Therefore, the EMT model was used to calculate the effective thermal conductivity of the wax layer in the subsequent simulation cases.

We then verify the numerical results with the measured values of all temperature detection points, including A–E and F–J. The details information of all temperature detection points are shown in

Table 4. The details of the detecting points.

Points	A	B	C	D	E	F	G	H	I	J
Measured (K)	290.3	293.8	292.2	289.4	287.2	286.2	288.4	287.3	286.1	285.2
Numerical (K)	289.858	294.125	292.031	289.081	286.359	285.697	288.881	287.763	285.629	284.361
RD (K)	0.442	0.325	0.169	0.319	0.841	0.503	0.481	0.463	0.471	0.839
ARD (%)	0.15	0.11	0.05	0.11	0.29	0.18	0.17	0.16	0.16	0.29

. Most calculations have quite a good agreement with the measured values. Though the maximum difference between the numerical and measured values is as large as 0.84 K, the average difference is only 0.49 K. We introduce a relative deviation (RD) defined as the measured value divided by the absolute difference between the calculated value and the measured value. With this definition, we can calculate that the average relative error (ARD) and maximum relative error (MRD) of all temperature detection points are 0.16% and 0.29%, respectively. From an engineering application point of view, the agreement is quite good.

4.2. Effect of Wax Layer with Specific Thickness

The thickness of the wax layer on the inner wall changes with the state of crude oil. We study the effect of wax layer with specific thickness on heat transfer of buried crude oil pipeline system. The cases including the wax layer with 0 mm and 4 mm thickness were simulated, respectively. In these cases, the inlet temperature and velocity of oil flow was 323.15 K and 0.7842 m/s. The Reynolds number is 20,733, and the convective heat transfer coefficient between air and ground surface was 25 W/(m∙K). The temperature field of pipeline cross section in the case with 0 mm and 4 mm wax layer are shown in Figure 7.

As can be seen from Figure 7, in the cross section, the contour varied greatly in the case of the model with a wax layer or not. According to the isotherms of the two cases that the area of the isotherm envelope was reduced in the case with 4 mm wax layer. Specifically, the lowest position of the 286.8 K isotherm was about 3.4 m when the 4 mm wax layer is taken into account, and was about 3.6 m without the wax layer. Similarly, compared with Figure 7a, the distribution range of 284.6 K isotherm in Figure 7b also reduced significantly. The intersection points between 284.6 K isotherm and the right half of the X-axis decreased from 0.7 m to 0.3 m. The results above proved that the influence of the wax layer on the heat transfer is not to be ignored. Due to a greater thermal resistance when the wax layer is formed on the inner wall of the pipe, the heat transfer from the oil to the soil is reduced. Therefore, to obtain an accurate temperature field, it is more reasonable to consider the wax layer in the thermal calculation of crude oil pipeline system.

We further discuss the quantitative influence of wax layer on soil temperature field. Figure 8 shows the temperature distribution along the centerline (auxiliary Path 1) on the pipeline cross section which considers the wax layer with 0 mm and 4 mm thickness. The range of the x-axis is (−4.5 m, 5.5 m) since the data points were along the path from the lower boundary to the upper boundary of the model. The temperature curve along the centerline was divided into five parts for discussion. In part I, the temperature remains unchanged from the lower boundary of the soil to the position of ordinate about 1 m, indicating that the heated pipeline has no obvious thermal influence on the soil in this part. In part II, the soil temperature increased gradually, indicating that the temperature of the soil is positively correlated with the distance between soil and pipe within the range of thermal influence. In part III, the temperature rose sharply due to a large temperature gradient caused by the high thermal resistance of the insulating layer. In part IV, the temperature reached the maximum value and remained constant, which corresponds to the temperature of the oil flow in the pipeline. In part V, temperature decreases gradually, indicating that the soil temperature is closer to the ambient temperature as it gets closer to the ground surface.

By comparing the temperature curves with 0 mm and 4 mm thickness wax layers, it can be concluded that the increase in wax layer thickness is negatively correlated with soil temperature outside the pipeline, and this trend is mainly manifested in part II and part V. Within a 1 m distance from the pipe, the average temperature deviation reached 1.42 K, and the maximum deviation of the two curves reached 2.90 K at x = 3.75, in part II.

4.3. Effect of the Oil Flow at Inlet

The heated crude oil dissipates heat during pipeline transportation, and the oil flow acts as the heat source and affect the formation of temperature field. The influence of oil flow on heat transfer of buried pipeline system was studied quantitatively by controlling the state of oil flow at the inlet. We first conducted a comparative test of oil flow with the same flow rate but different temperatures. In the three simulation cases adopted, the oil flow temperature at inlet was set at 323 K, 332 K, and 341 K, respectively. The Reynolds numbers were 21,151, 26,952, and 33,226, respectively. The rest calculation conditions were the same, including oil velocity (0.8 m/s), the convective heat transfer coefficient at the soil surface (25 W/(m·K)), and the soil thermostatic layer (283 K). Figure 9 shows the temperature fields on the cross sections at different oil flow temperatures.

It can be concluded from the contours that the shape of the isotherm near the pipeline changes obviously when the temperature of oil flow increases, but is not obvious in the parts that are far from the pipe. Here we introduce the definition of core thermal response zone. It refers to the area where the temperature gradient exceeds 1.5 K/m of the soil temperature field on the cross section. According to the isotherm distribution, the temperature values within the coordinate range of (−10 m, −6 m) and (6 m, 10 m) change a little. In the two regions, the temperature profile was close to the soil in the natural state. Therefore, the range of the core thermal response zone on the X-axis can be set within (−6 m, +6 m). Based on this conclusion, it is feasible to reduce the size of the X-axis appropriately when the computation needs to be simplified by reducing the computational domain.

Figure 10 shows the temperature distribution along the centerline on the cross section with different oil temperatures at inlet. As can be seen, the higher the oil flow temperature, the greater the soil temperature near the pipeline. However, only the soil temperature within 2 m of the pipeline has a significant positive correlation with oil temperature, the temperatures of the three curves are virtually unanimous beyond 6 m from the pipe.

Figure 11 shows the temperature distribution at different depths in the case of 341 K oil flow inlet temperature. As can be seen, with the increase in depth, soil temperature fluctuation becomes smaller. For the lines at the depth of 6 m and 7 m, the temperature fluctuates a lot, and the maximum amplitudes are 2.2 K and 1.5 K, respectively. While the temperature curves show slight fluctuations at 9 m depth and 8 m depth, where amplitude does not exceed 1 K and the peaks are consistent. In other words, soil temperatures at the depth of 8 m and 9 m from the ground surface are very close regardless of oil flow temperature changes.

We further discuss the influence of velocity on heat transfer of buried pipeline system. The velocity is an important parameter determined by the throughput of the pipeline, and the change of velocity will affect the intensity of convective heat transfer. Thus, we carried out a comparison of two simulation cases, in which the velocity inlet was set as 0.9 m/s and 1.2 m/s, respectively. In addition, the Reynolds numbers were 23,893 and 31,860, respectively. Contours of two cases at different velocity were shown in Figure 12.

As can be seen from the contours, the range of the core thermal response zone increases significantly with the increase in velocity. Specifically, when the velocity increases from 0.9 m/s to 1.2 m/s, the maximum value of the 283.1 K isotherm in the X direction expands from 3.5 m to 5.2 m. Therefore, it is only appropriate to reduce the computational domain in the condition of a small velocity of oil flow.

This work uses predefined points along the auxiliary lines to analyze temperature fields in detail, including the centerline and horizontal line at different buried depths on the cross section. Figure 13 shows the temperature distribution along the centerline on the cross section in the cases of different oil flow rates. As can be seen, the greater the velocity, the greater the soil temperature near the pipeline. It is worth noting that the increase in velocity has a significant thermal effect on soil that is far from the pipeline, indicating that the enhanced flow can increase the scope of the core thermal response zone.

In Figure 14, the temperature curves along different burial depths were made in the simulation cases with 0.9 m/s and 1.2 m/s velocity. As can be seen, at the same velocity, the peak of the temperature curve is negatively correlated with buried depth; while the peak of the temperature curve is positively correlated with the flow at the same buried depth. The peaks of the two temperature curves at 8 m depth are equal, indicating that the thermal influence of oil flow is weak when it is far away from the pipeline.

4.4. Effect of the External Environment

The change of external environment will affect the formation of temperature field of buried pipeline system. We study the effect of the external environment in order to make the reasonable transportation plan of the crude oil. Since the soil temperature fields tend to change slowly, the alternation of four seasons in a year was used to represent the fluctuation of the external environment, and suggestions were provided for hot oil pipeline transportation. In the simulation scheme, spring and autumn were combined, and the boundary conditions were set to match the characteristics of different seasons. The settings of numerical cases for different seasons are listed in

Table 5. The settings of numerical cases for different seasons.

	Parameters	Air Temperature K	Convective Heat Transfer Coefficient W/(m∙K)	Soil Thermostatic Layer K	Reynolds Number
Seasons
Spring/Autumn		276	20	283	25,220
Summer		293	30	288	25,225
Winter		253	30	281	25,223

, including the temperature of air, the convective heat transfer coefficient at the upper boundary, the temperature of soil thermostatic layer at the lower boundary, and the Reynolds number.

Figure 15 shows the contours of the pipeline cross section in different seasons. It can be seen that the isotherms in different seasons have significant changes due to the act of temperature gradient. The core thermal response zone of the heat pipe presents the largest scope in summer and the smallest scope in winter. It suggests the external environment has the greatest influence on the temperature field of buried pipeline system in winter.

Through the temperature distribution along the centerline on the cross section in the cases of different seasons, we further analyzed the influence of external environment on the distribution of temperature field. As can be seen from the Figure 16, the soil temperature at the same position is much higher in summer than in other seasons. Except in winter, the temperature distribution along the centerline below the pipeline shows a trend of unilateral increase. However, in winter, the soil temperature curve beneath the pipe varies complicated, and its lowest value appears at about 0.9 m below the oil pipe. Therefore, attention should be paid to the endothermic effect of external environment on the temperature field of buried pipeline system in winter, and the temperature of oil flow should be appropriately raised.

Figure 17 shows soil temperature distribution at a depth of 0.5 m in different seasons. In the main part of the curves, the values are close to air temperature and only rise in the range of −2 m to 2 m along the X-axis. It indicates that the temperature at the depth of 0.5 m is mainly determined by the environment. With the help of the temperature curve, we can estimate the soil temperature near the pipeline. For instance, in summer, the soil temperature at 0.5 m depth which is directly above the pipe exceeds the air temperature by about 3 K.

Figure 18 shows soil temperature distribution at a depth of 8 m and 9 m in different seasons. It can be seen from the Figure 18 that in different seasons, although the curves of the two depths have different trend and deviation ways, the maximum deviation of the temperature curves at 8 m depth and 9 m depth is less than 2.5 K. In other words, no matter the external environment changes, the soil temperature distributions at 8 m and 9 m depth are similar.

5. Conclusions

The CFD modeling successfully predicted the temperature field of buried crude oil pipeline. Our numerical results well agree with the real measured data. The heat transfer effect caused by the wax layer and the thermal response law of influence factors can be summarized as follows.

(1)

The numerical results proved that the effect of the wax layer on heat transfer cannot be ignored. Thermal resistance caused by wax layer will change the temperature field distribution of buried pipeline system. If the wax layer is ignored, the average deviation of soil temperature is 1.42 K and the maximum deviation is 2.90 K within 1 m away from the pipe. The thermal conductivity of waxy layer can be predicted by using heterogeneous material thermal conductivity model, and the EMT model has the best accuracy and applicability. When the mass fraction of wax crystal is in the range of 30% to 65%, the effective thermal conductivity of the wax layer is in the range of 0.155–0.195 W/(m∙K).

(2)

According to the comparison of simulation results under different working conditions, the thermal influence of heat pipe is positively correlated with oil flow temperature and oil flow rate. The core thermal response zone is defined as the area where soil temperature varies significantly with oil flow. The range of this region on the X-axis is suggested to be within the (−6 m, 6 m) coordinates. When the flow rate increased from 0.9 m/s to 1.2 m/s, the core thermal response zone advanced about 3.4 m in the X-axis direction. Based on the simulation results of different seasons, the soil temperature at the depth of 0.5 m is almost equal to the ambient temperature; and the soil temperature beyond 8 m depth from the ground fluctuates little so that it can be viewed as constant.

Our methods and results provide useful reference for long-distance oil pipeline transportation, and further research based on this work is promising. The results of this work provide initial conditions and a theoretical basis for the study of complex conditions of the buried crude oil pipeline, such as pipeline shutdown and restart.