Study on Asymmetry Concentration of Mixed Oil in Products Pipeline

: Long-distance pipelines transporting multiple product oils such as gasoline, diesel and jet fuel, are important facilities for transporting fossil energy. One major concern in operation is the energy consumption of the pipeline. Energy consumption should be made optimized tracking batches of oils and cutting mixed oil, which requires an accurate prediction of concentration curve. In engineering, the concentration curve is usually assumed to be symmetric, but it is actually asymmetric, which may lead to estimation errors. Thus, the asymmetric concentration of mixed oil should be studied. The formation mechanism of the asymmetry of concentration curve has not been clearly clariﬁed. A new method is proposed to measure the asymmetry of the concentration curve. Quantitative analysis is carried out for each factor on the asymmetry distribution of concentration curve. Based on the convection–di ﬀ usion equation, a modiﬁed oil-mixing model considering near wall adsorption e ﬀ ect is established. The model shows a good agreement with the Jablonski empirical formula. The error, compared with the experimental results, is less than 5%. The main ﬁndings are: (1) deviation volume has a negative correlation with pipe diameter and mean velocity; (2) adsorption coe ﬃ cient has a greater impact on the length ratio of front and tail oil than di ﬀ usion coe ﬃ cient; (3) the inﬂuence of all factors considered on the total length of mixed oil, front oil, tail oil and trail oil are basically the same; (4) if the limit of adsorption concentration in adsorption layer is 1, the reasonable value of adsorption coe ﬃ cient a and b should be around 0.4. The results reveal the mechanism of asymmetric concentration of product oils and can provide practical suggestions to deal with the mixed oil.


Introduction
An efficient way to transport different types of product oils over a long distance is to use one pipeline to transport the different oils sequentially because the construction of multiple pipelines would require much larger investments and energy consumption. Consecutive oils in pipeline form one interfacial section containing mixed oil [1]. A long-distance pipeline has several mixed oil sections. Oil transportation enterprises require the mixed oil sections to be as few and short as possible, to ensure a high-quality product. This requires scheduling, batch interface tracking and mixed oil cutting in engineering. The complexity of these operations increase with the increasing types of product oil and the branch points [2]. Due to the high efficiency of computer-aided programming of this complex problem, the requirements for numerical simulation of mixed oil flow in products pipelines have gradually increased [3]. As the positions of the front oil and tailing oil cannot be accurately obtained during the operation of the products' pipeline, the mixing situation of the oils is difficult to obtain quickly and accurately. This results in the predicted value of the oil quality differing greatly from the actual value, which deeply affects the automated operation of the product pipeline [4].
Due to the non-uniform distribution of flow velocity on the pipe cross-section and the adsorption of oil molecules on the pipe wall, the oil-mixing concentration curve is usually asymmetrical, which greatly affects the accuracy of mixed oil cutting. The asymmetry of the oil concentration was observed by early scholars in experiments, but it can hardly be predicted accurately [5]. With the rapid developments of computers, an increasing number of scholars choose to study the mixed oil problems through numerical simulations to obtain the asymmetric phenomena of the mixed oil concentration [6]. However, the formation mechanism of the asymmetric concentration has not been revealed. With the increasing requirements of product pipeline scheduling [7] and mixed oil cutting [8], it is very important in practical operation to predict the mixed oil concentration curve quickly and accurately [9]. The fast prediction depends on an efficient mathematical model reflecting the asymmetry of oil concentration. If oil flow is treated as a one-dimensional flow, the variation in oil concentration can only be considered along the axis of the pipeline. It cannot describe the asymmetry of oil concentration in radius direction. Therefore, mathematical models of asymmetric distribution of mixed oil should be two-dimensional or three-dimensional.
Studies using two-dimensional model are summarized as follows. Lu [10] proposed a time-splitting method to decompose the two-dimensional convection-diffusion equation into a convection equation and a diffusion equation. The convection equation was solved by the characteristic line method, while the diffusion equation was solved by the finite difference method. This treatment is easy to converge. Deng and Pu [11] proposed that oil mixing was the result of the comprehensive action of radial diffusion and axial extension. The simulation results showed that the tail of mixed oil was longer than the head. The concentration of mixed oil was asymmetrically distributed. Engineering factors that affect the tailing of mixed oil were listed but the mechanism was not analyzed. Zheng and Jiang [12] and Ma and Bai [13] also found asymmetric distribution of mixed oil concentration using two-dimensional cylindrical coordinates and an adaptive grid. Wang [14] used a commercial software to establish a two-dimensional flow and mass transfer coupling model for alternating oil transportation. He analyzed the mixed oil process of gasoline as the front oil and diesel as the tailing oil in the batch transportation using numerical simulation, and obtained the relationship of the mixed oil head shape with time.
Some other scholars tried three-dimensional models. Zhao [15] established a mixed oil model based on computational fluid software by using the Reynolds time-averaged method. Simulation results show the asymmetry of the mixed oil concentration curve. In his opinion, the asymmetry is caused by flushing action of the tailing oil on the front oil, the non-uniform velocity distribution, and the initial mixed oil. Zhang [16] simulated the oil mixture and reached the conclusion that the oil mixture concentration curve is asymmetric. Wang et al. [17] considered the influence of turbulence and obtained the asymmetric concentration curve. They concluded that the cutting effect of mixed oil was better if the front oil was gasoline rather than diesel. Dai et al. [18] considered the influence of gravity on oil mixing. They found that the length of the mixed oil will increase with increasing inclination angle when diesel is in front of gasoline. When gasoline is in front of diesel, the length of mixing oil shortens with increasing inclination angle. He et al. [19] simulated the turbulent flow of gasoline and diesel in the inclined pipeline and obtained similar conclusions with Dai et al. [20].
The three-dimensional models can describe the flow details more precisely. However, they cannot be used for long-distance pipelines due to their much longer computational time compared to the two-dimensional models. Moreover, the asymmetry is mainly in the radial direction. Thus, the two-dimensional model, balancing precision and computational speed is a better choice. It is used in this paper to study the asymmetry of mixed oil concentration. From the literature review, we realize Energies 2020, 13, 6398 3 of 24 that researchers mentioned the phenomenon of the asymmetry of mixed oil concentration but did not analyze its mechanism. We are trying, in this paper, to study the asymmetry mechanism by proposing a method to measure the asymmetry and quantify the relevant influencing factors. It is very useful for guaranteeing the precise operation of multiproduct pipelines.

Flow-Control Equation
The pipe flow of oils is actually three-dimensional (3D). However, oil concentration varies mainly in the radial and axial directions, and therefore the transport of the oil concentration can be considered as a two-dimensional (2D) convection-diffusion problem. The basic form of the equation is where, u is the flow velocity, m/s; c is the oil concentration; D is the turbulent diffusion coefficient, m 2 /s; r is the radial distance, m; x is the axial distance, m; t is the time, s. This equation shows the general law of convection-diffusion of oil concentrations in radial and axial directions. The left-hand side of the equation contains the partial derivative of concentration with respect to time and convection terms, while the right-hand side of the equation contains the diffusion term. In the two-dimensional model, the values of oil concentration, flow velocity and diffusion coefficient are different at the same section and different radii of the pipeline.
Since the asymmetry of oil distribution is related to the non-uniform radial distribution of velocity, Equation (1) can be expressed in different forms in different regions of the pipe cross-section. For the fully developed turbulent flow in the pipe, the flow region can be divided into a viscous sublayer, buffer layer and turbulence core, as shown in Figure 1.
Energies 2020, 13, x FOR PEER REVIEW 3 of 26 a method to measure the asymmetry and quantify the relevant influencing factors. It is very useful for guaranteeing the precise operation of multiproduct pipelines.

Flow-Control Equation
The pipe flow of oils is actually three-dimensional (3D). However, oil concentration varies mainly in the radial and axial directions, and therefore the transport of the oil concentration can be considered as a two-dimensional (2D) convection-diffusion problem. The basic form of the equation is where, u is the flow velocity, m/s; c is the oil concentration; D is the turbulent diffusion coefficient, m 2 /s; r is the radial distance, m; x is the axial distance, m; t is the time, s. This equation shows the general law of convection-diffusion of oil concentrations in radial and axial directions. The left-hand side of the equation contains the partial derivative of concentration with respect to time and convection terms, while the right-hand side of the equation contains the diffusion term. In the twodimensional model, the values of oil concentration, flow velocity and diffusion coefficient are different at the same section and different radii of the pipeline.
Since the asymmetry of oil distribution is related to the non-uniform radial distribution of velocity, Equation (1) can be expressed in different forms in different regions of the pipe cross-section. For the fully developed turbulent flow in the pipe, the flow region can be divided into a viscous sublayer, buffer layer and turbulence core, as shown in Figure 1. The influence of convection and diffusion in different regions varies greatly, so different empirical formulas are needed to calculate the velocity field and turbulent diffusion coefficient [20]. The simplified formula for each region is as follows.
(1) Viscous sublayer According to the no-slip boundary condition, the viscosity of the fluid near the pipe wall is extremely significant, the oil flow velocity is zero, and there is always a thin viscous sublayer on the pipe wall. This area satisfies the following conditions where, y is the distance from pipe wall, m;  is the kinematic viscosity of oil, m 2 /s; ρ is the density of oil, kg/m 3 ; 0  is the wall shear stress, Pa.
The velocity distribution in viscous sublayer is as follows The influence of convection and diffusion in different regions varies greatly, so different empirical formulas are needed to calculate the velocity field and turbulent diffusion coefficient [20]. The simplified formula for each region is as follows.
(1) Viscous sublayer According to the no-slip boundary condition, the viscosity of the fluid near the pipe wall is extremely significant, the oil flow velocity is zero, and there is always a thin viscous sublayer on the pipe wall. This area satisfies the following conditions where, y is the distance from pipe wall, m; ν is the kinematic viscosity of oil, m 2 /s; ρ is the density of oil, kg/m 3 ; τ 0 is the wall shear stress, Pa. The velocity distribution in viscous sublayer is as follows The diffusion mechanism is only molecular diffusion. Therefore, the diffusion coefficient D in this area can be simplified to molecular diffusion coefficient D L , i.e., D = D L . Then, the two-dimensional convective-diffusion equation is simplified as follows Equation (4) is used to calculate the oil concentration in the viscous sublayer.
(2) Buffer layer The buffer layer refers to the transition area from the viscous sublayer to the turbulent core area. The flow area of the buffer layer is divided as follows The velocity distribution in the buffer layer is as follows [21] The influence of both molecular diffusion and turbulent diffusion in this area cannot be ignored. Therefore, the diffusion coefficient is taken as the sum of molecular diffusion coefficient and turbulent diffusion coefficient Therefore, the convection-diffusion equation in this region can be simplified to the following formula Equation (8) is the equation to calculate the oil concentration in the buffer layer.
(3) Turbulent core region When the fluid flows turbulent in a circular tube, most of the fluids are in the turbulent core region [22], and the fluids in the above two regions only account for a small part. The turbulence core region contains the following regions [21] y > 30ν/ τ 0 /ρ (9) The velocity distribution in this area is as follows In the turbulent core, as the turbulent diffusion is far greater than the molecular diffusion, the turbulence diffusion coefficient is about an order of magnitude greater than the molecular diffusion coefficient. Therefore, it is considered that the diffusion coefficient in the equation is equal to the turbulent diffusion coefficient, and the calculation formula is as follows Energies 2020, 13, 6398

of 24
This diffusion coefficient can be substituted into Equation (1) so that the oil concentration in turbulence core region can be calculated.
For all the above equations, the initial conditions are The boundary conditions are as follows where c represents the oil concentration, %; x represents the axial position of the calculated pipe section, m.

Model with Polar Adsorption Layer
The equations in Section 2.1 did not consider the adsorption effect. Diesel and gasoline can form a polar membrane on the pipe wall. The existence of the membrane is easy to cause a longer trace of the mixed oil tail. Thus, it is important to consider the polar adsorption layer in the calculation of the mixed oil concentration. The effect of adsorption is reflected by the interaction between the pipe wall and oil molecules. In the actual adsorption process, there are usually 3-4 layers of molecules in the polar adsorption layer, among which, the first layer of molecules closest to the pipe wall has the strongest adsorption capacity. For the convenience of calculation, the polar adsorption layer is considered as monolayer adsorption. The dynamic process of monolayer adsorption can be expressed by Equation (14), while the desorption process can be expressed by Equation (15) where c is the oil concentration in the pipeline, %; c r is the concentration when equilibrium is reached in the polar adsorption layer; K a is adsorption constant, s −1 ; c r * is the maximum concentration that can be reached in the adsorption layer without desorption, % where K d represents the desorption constant, s −1 . By comprehensively considering the dynamic equilibrium process of adsorption and desorption, the following equation can be obtained According to the assumed solution conditions, the concentration of B oil product is c; At the initial time, the adsorption concentration in the adsorption layer is 0, which can be obtained as follows It is worth noting that only when desorption does not occur can the concentration in the limiting adsorption layer reach the limit adsorption concentration c r * , where K d = 0. When the adsorption and Energies 2020, 13, 6398 6 of 24 desorption processes occur simultaneously and reach dynamic equilibrium after a long enough time, the concentration in the adsorption layer is where the values of parameters a and b are to be determined. If the time required for adsorption and desorption to reach dynamic equilibrium is much less than the time required for oil replacement, the adsorption concentration in the polar adsorption layer is The derivative is Substitute Equation (20) into the convection-diffusion equation The above equation is a two-dimensional model considering a single layer of polar adsorption on the pipe wall. Equation (21) can be used to calculate the oil concentration in the polar adsorption layer.

Model Considering Temperature Change
The density and viscosity of oil are related to temperature, so that temperature has an important role in batch transportation of products [23]. There is a linear relationship between the density and temperature. For liquids, in the absence of phase transitions, the higher the temperature, the less dense the liquid is. The specific conversion formula of density is as follows [5] ρ t = ρ 20 − ε(t − 20) (22) where ρ t , ρ 20 are the oil density at t • C and 20 • C, kg/m 3 ; ε is the temperature coefficient ε = 1.825 − 0.001315ρ 20 , Kg/(m 3 · • C). The viscosity of oil at any temperatures can be obtained by converting the viscosity-temperature index of oil to the known viscosity at a certain temperature. In the absence of phase transitions, the higher the temperature, the smaller the viscosity is for a Newtonian fluid. The specific formula for viscosity conversion is as follows [5] v where v t represents the viscosity of the oil at t • C, m 2 /s; v 0 represents the viscosity of the oil at the temperature of t 0 • C, m 2 /s; u is the viscosity-temperature index of the oil.
The influence of temperature on Reynolds number is indirectly caused by the influence of temperature on viscosity. The influence of temperature on Reynolds number is as follows The influence of temperature on friction coefficient The influence of temperature on the shear stress of the pipe wall The effect of temperature on the diffusion coefficient is the same as that on the Reynolds number, so it will not be listed separately here. In the model considering temperature, the influence of temperature can be obtained only by replacing the parameters affected by temperature in the formula with the function expression related to temperature.

Computational Domain
As the external node method is adopted in this paper, the grid division under the column coordinate system is shown in Figure 2. In the figure, the horizontal axis is the axial direction of the pipeline, containing a total of Nx nodes, and the vertical axis is the radial direction of the pipeline, containing a total of Nr nodes. When the boundary layer mesh is encrypted, the mesh between the Nr -1 and Nr layers is encrypted, and the outermost encrypted mesh contains Nr b nodes.
0 v e  The influence of temperature on friction coefficient The influence of temperature on the shear stress of the pipe wall The effect of temperature on the diffusion coefficient is the same as that on the Reynolds number, so it will not be listed separately here. In the model considering temperature, the influence of temperature can be obtained only by replacing the parameters affected by temperature in the formula with the function expression related to temperature.

Computational Domain
As the external node method is adopted in this paper, the grid division under the column coordinate system is shown in Figure 2. In the figure, the horizontal axis is the axial direction of the pipeline, containing a total of Nx nodes, and the vertical axis is the radial direction of the pipeline, containing a total of Nr nodes. When the boundary layer mesh is encrypted, the mesh between the Nr-1 and Nr layers is encrypted, and the outermost encrypted mesh contains Nrb nodes. In the discrete format of the calculation area, the calculation formula of the average oil concentration at the pipe section is as follows where c represents the average oil concentration at a certain axial node of the pipeline crosssection; i c represents the oil concentration value at the i-th radial node on the section; i V represents the control volume and circumferential direction at the corresponding node i c . In the discrete format of the calculation area, the calculation formula of the average oil concentration at the pipe section is as follows where c represents the average oil concentration at a certain axial node of the pipeline cross-section; c i represents the oil concentration value at the i-th radial node on the section; V i represents the control volume and circumferential direction at the corresponding node c i .

Convection Equation Discretization
The form of the pure convection equation is as follows. Based on the basic form of the pure convection equation, the characteristic line method is used to solve the pure convection equation.
Energies 2020, 13, 6398 8 of 24 The above equations are discretized along the direction of characteristic lines. Figure 3 shows the principle of characteristics method. The slanted dotted line in the figure represents the direction of the characteristic line. The n-time layer is related to the concentration value on the (n + 1)-time layer, t n represents the n-time layer, t n+1 represents the (n + 1)-time layer, and x represents the axial coordinates of the pipeline. j, k is an integer, a is a real number, and ∆t represents the time interval between the n-time layer and the (n + 1)-time layer.

Convection Equation Discretization
The form of the pure convection equation is as follows. Based on the basic form of the pure convection equation, the characteristic line method is used to solve the pure convection equation.
The above equations are discretized along the direction of characteristic lines. Figure 3 shows the principle of characteristics method. The slanted dotted line in the figure represents the direction of the characteristic line. The n-time layer is related to the concentration value on the n + 1-time layer, tn represents the n-time layer, tn+1 represents the n + 1-time layer, and x represents the axial coordinates of the pipeline. j, k is an integer, a is a real number, and t  represents the time interval between the n-time layer and the n+1-time layer. When only convection terms are considered, along the direction of the characteristic line Then, to calculate the concentration value of node j at n + 1-time layer, it is only necessary to calculate the concentration value at point j-a at n-time layer. Because a is a real number greater than zero, j-a may not fall on the divided grid nodes. Here, the concentration value at the j-k point of the grid node is used to approximate the concentration value at the j-a node. Suppose the relationship between the concentration distribution and the distance at n-time layer is as follows Through the above formula, the concentration value of node j at n + 1-time layer can be obtained, and so on. The concentration value at each node, except the beginning of the pipeline (boundary condition) at n + 1-time layer, can be calculated.
The solution of the two-dimensional convection term is the same as the one-dimensional principle, except that the velocity in the two-dimensional pure convection equation is no longer the Then, to calculate the concentration value of node j at (n + 1)-time layer, it is only necessary to calculate the concentration value at point j-a at n-time layer. Because a is a real number greater than zero, j-a may not fall on the divided grid nodes. Here, the concentration value at the j-k point of the grid node is used to approximate the concentration value at the j-a node. Suppose the relationship between the concentration distribution and the distance at n-time layer is as follows Through the above formula, the concentration value of node j at (n + 1)-time layer can be obtained, and so on. The concentration value at each node, except the beginning of the pipeline (boundary condition) at (n + 1)-time layer, can be calculated.
The solution of the two-dimensional convection term is the same as the one-dimensional principle, except that the velocity in the two-dimensional pure convection equation is no longer the average velocity of the pipe cross section, but the specific velocity at a certain point of the pipe cross-section. In the calculation, the velocity can be considered in different regions.
In order to reduce the numerical diffusion effect of Equation (32), Casulli adopts one upwind point and two downwind points to carry out quadratic difference, and its calculation formula is as follows Energies 2020, 13, 6398 9 of 24 When k = 0, the above formula is transformed into Lax-Wendroff format Equation (34) does not introduce numerical diffusion effect, but there is a dispersion phenomenon. Equations (33) and (34) are weighted and averaged to obtain Equation (35), which is effective in most cases.

Diffusion Equation Discretization
The pure diffusion equation form of the one-dimensional mixed oil model is as follows Compared with the explicit format, the implicit format solution has the advantage of absolute stability. The following describes the method of implicitly solving the diffusion term. The format of implicit discretization is as follows The recursive formula after organizing is as follows From the above formula, the relationship between the concentration of a certain point at the n-time layer and the concentration of three points at the (n + 1)-time layer can be obtained. The concentration value marked n above represents the concentration value of the known time layer, while the concentration value at (n + 1)-time layer has yet to be solved.
When n = 1 and i = 2, the above equation becomes where c 2 1 is the boundary condition at the beginning of the pipeline, which is a known value. Move this term to the column vector at the right end of the equation, and use the same method to process the points at the end of the pipeline to obtain the following matrix. When the number of axial nodes in the pipeline is Nx, the tridiagonal matrix on the left is ( where, c n+1 i is approximated by c n+1 i = c n+1 i−1 according to the boundary conditions at the end of the pipeline. Then, the matrix becomes the following form; pay attention to the change in the coefficient in the lower right corner of the tridiagonal matrix and the change in the last term of the column vector on the right side of the equation With the help of the above tridiagonal matrix, the concentration value of each node at the next time layer can be solved according to the known concentration value of each node at the previous time layer.
The solution method of the two-dimensional diffusion term is similar to that of the one-dimensional one, which uses the form of diagonal matrix to relate the concentration value at the previous time layer with that of the next time layer to obtain a series of closed equations, and solve the whole equation at the same time step. Because the two-dimensional diffusion term includes both axial and radial directions, the diagonal matrix is more difficult to solve than the one-dimensional case. Since the solution is a global solution at the whole-time layer, all the internal nodes corresponding to the axial and radial directions at the whole-time layer need to be transformed into a column vector, so that the global solution can be performed at a time layer.
Due to the complexity of the penta-diagonal matrix, the penta-diagonal matrix is divided into two matrices with the same size, and the penta-diagonal matrix is the sum of the two matrices, namely Q = Q 1 + Q 2 . Now, Q 1 and Q 2 are described, respectively.
Firstly, the elements to be used in the matrix are defined In the expressions of the above three elements m, n and k, D is the diffusion coefficient. The diffusion coefficient of the model used in this paper is related to the position of the node in the pipeline, so the value of D at different nodes is different, which is uniformly represented by D here.
The matrix Q 1 is composed of smaller matrices A and B. The expression for Q 1 is as follows, where the value of unlisted elements is zero, the same as below where A and B are matrix composed of m, n, and k.
In the matrix Q 2 , the values of elements on two lines parallel to the main diagonal are not equal to zero, the values of the elements on the two lines are the same, and the two lines are symmetric regarding the main diagonal of Q 1 . The distance between the non-zero element m in the first row and the first element in the first row is Nx − 2 to determine the position of the two lines.
The matrix P on the right side of the equation represents the known concentration value corresponding to each node at the time layer. The matrix on the right side of the equation (denoted as P here) has the following form P = P 1 + P 2 (49) P 1 and P 2 can be expressed by column vectors, respectively, as T (50) where, c i,j represents the oil concentration value of the rear oil corresponding to the radial node i and the axial node j.
The matrix M to be solved represents the concentration value of each node at the new time layer, and its expression is as follows Then, the whole system is expressed in the form of matrix product as Energies 2020, 13, 6398 12 of 24

Boundary Condition Dispersion
At the pipe axis and the pipe wall The discretized format is t ≥ 0 : where c represents oil concentration, %; x represents the axial position of the calculated pipe segment, m.

Proposal and Calculation Method of Deviation Volume
At present, there is no complete method to measure the overall asymmetry of the mixed oil concentration curve. In this paper, the concept of deviation volume is proposed to characterize the degree of asymmetry in the concentration curve of mixed oil. Figure 4a shows the curve of mixed oil concentration; the horizontal axis is the length of mixed oil. The vertical axis is the concentration of the rear oil. The interface of the mixed oil section is asymmetric with respect to the 50% concentration.
To measure the degree of asymmetry of the mixed oil concentration curve, the curve is divided into two parts: the front part and the back part. The front part, i.e., the part where the mixed oil head is located, is the left part of the curve, where the oil concentration before it is greater than or equal to 50%. The back part, i.e., the part where the mixed oil tail is located, is the right half of the curve, where the oil concentration of rear oil is greater than 50%.
The volume of the rear oil products contained in the first half of the mixed oil section is V f represents the volume of the rear oil contained in the first half of the mixed oil; the lower index f represents the front-part of the mixed oil; the upper index 0.5 represents that the concentration of the front oil is greater than or equal to 0.5; R represents the radius of the pipeline; L f represents the length of the first half of the mixed oil; c i represents the concentration of the rear oil. Similarly, the volume of the front oil contained in the second half of the mixing oil section is V t represents the volume of the front oil contained in the second half of the mixed oil; the lower index t represents the tail-part of the mixed oil; the upper index 0.5 represents that the concentration of the rear oil is greater than 0.5; R represents the radius of pipeline; L t represents the length of the second half of the mixed oil.
Flip the first half of the oil mixture concentration curve, flip and translate the second half, as shown in Figure 4b. The two halves are similar in shape but do not overlap. The area of the gap formed between the two curves is used to represent the degree of asymmetry of the two oil-mixed curves, as shown in the shaded part of Blue in Figure 4c. The area is the same as the length unit, and the unit is the volume after multiplying by the pipe section area. The area of the blue shaded part in Figure 4c is defined as the deviation volume. Its physical meaning is the difference between the volume of the front oil in the second half of mixing oil and the volume of the rear oil in the first half of mixing oil. The calculation method is where V d represents the deviation volume, V f represents the volume of the rear oil in the first half of the mixed oil, and V t represents the volume of the front oil in the first half of the mixed oil. The deviation volume can reflect the asymmetry of the entire mixed oil concentration curve, and has a clear physical meaning.
The absolute value of the deviation volume indicates the degree of asymmetry of the concentration curve. When the deviation volume value is positive, it means that the second half of the mixed oil contains more front oil than the rear oil contained by the first half. When the deviation volume is negative, it means that the second half of the mixed oil contains less front oil than the rear oil contained by the first section.

Basic Data
The relative parameters of oil density and viscosity used in the calculation process are shown in Tables 1 and 2.   When analyzing the difference between the mixing oil head and the mixing oil tail, the volume of the rear oil contained in the mixing oil head is The volume of the front oil contained in the mixing oil tail is Then, the difference between mixed oil head tail can be described by The upper index 0.1 indicates that the comparison is the volume difference between the front oil contained in the mixed oil tail and the rear oil contained in the mixed oil head. From the graph of the concentration curve, it represents the deviation volume corresponding to the green part in Figure 4d. The dimensionless deviation volume is used in this paper. The characteristic length of the pipe is the pipe diameter D, and the dimensionless deviation volume is defined as The calculation method is The absolute value of the deviation volume indicates the degree of asymmetry of the concentration curve. When the deviation volume value is positive, it means that the second half of the mixed oil contains more front oil than the rear oil contained by the first half. When the deviation volume is negative, it means that the second half of the mixed oil contains less front oil than the rear oil contained by the first section.

Basic Data
The relative parameters of oil density and viscosity used in the calculation process are shown in Tables 1 and 2.

Model Validation
To validate the above models, we test the grid independence of numerical results, and compare the numerical results with physical experiments and famous empirical formulas as follows.

Grid Independence Test
The grid independence test of the model is shown in Figure 5. Theoretically, the calculation results are more accurate using a denser mesh, but too dense a mesh may lead to a too-long computational time and insufficient computer memory. As shown in Figure 5, when dx and dr continue to decrease, the mixed oil concentration curve of the products pipeline changes little, and the length of the mixed oil section of the products pipeline changes very slightly. To improve the calculation efficiency, this paper finally selects the model grid parameters dx = 0.025 m, dr = 0.007 m for simulation. results are more accurate using a denser mesh, but too dense a mesh may lead to a too-long computational time and insufficient computer memory. As shown in Figure 5, when dx and dr continue to decrease, the mixed oil concentration curve of the products pipeline changes little, and the length of the mixed oil section of the products pipeline changes very slightly. To improve the calculation efficiency, this paper finally selects the model grid parameters dx = 0.025 m, dr = 0.007 m for simulation.

Experimental Verification
The simulation results of the two-dimensional model were compared with the loop data obtained by the research team. A water-sodium chloride aqueous solution was used as the experimental medium. The length of the experimental pipe is 200 m. The length of the numerical simulation pipe is 200 m, and the simulation medium and pipe parameters are the same as the experiment. The length of the mixed liquid measured in the experiment is 9.2 m. The numerically calculated length of the mixed liquid is 9.5 m. Their deviation is as low as 3.3%. The three solid dotted lines in Figure 6

Experimental Verification
The simulation results of the two-dimensional model were compared with the loop data obtained by the research team. A water-sodium chloride aqueous solution was used as the experimental medium. The length of the experimental pipe is 200 m. The length of the numerical simulation pipe is 200 m, and the simulation medium and pipe parameters are the same as the experiment. The length of the mixed liquid measured in the experiment is 9.2 m. The numerically calculated length of the mixed liquid is 9.5 m. Their deviation is as low as 3.3%. The three solid dotted lines in Figure 6 Figure 7 shows the comparison result of the mixed oil quantity obtained by the numerical simulation and the empirical formula. The total length of the simulated pipeline is 30 km, and the simulated transportation temperature is 20 °C, regardless of temperature changes along the pipeline. The order of oil delivery is 0# diesel first and 90# gasoline later. The pipeline flow rate is 1 × 10 3 m 3 /h. The numerical simulation results are the closest to Jablonskij's empirical formula [23].  Figure 7 shows the comparison result of the mixed oil quantity obtained by the numerical simulation and the empirical formula. The total length of the simulated pipeline is 30 km, and the simulated transportation temperature is 20 • C, regardless of temperature changes along the pipeline. The order of oil delivery is 0# diesel first and 90# gasoline later. The pipeline flow rate is 1 × 10 3 m 3 /h. The numerical simulation results are the closest to Jablonskij's empirical formula [23]. Figure 7 shows the comparison result of the mixed oil quantity obtained by the numerical simulation and the empirical formula. The total length of the simulated pipeline is 30 km, and the simulated transportation temperature is 20 °C, regardless of temperature changes along the pipeline. The order of oil delivery is 0# diesel first and 90# gasoline later. The pipeline flow rate is 1 × 10 3 m 3 /h. The numerical simulation results are the closest to Jablonskij's empirical formula [23].

Comparison of One-Dimensional Model and Two-Dimensional Model
The solution results of the one-dimensional model and the two-dimensional model are shown in Figure 8. The red dashed line in the figure is the calculation result of the one-dimensional model. After testing, the deviation volume is 0, which is a center-symmetric figure. The blue and black solid lines, respectively, indicate the changes in the concentration of oil at the axis of the pipeline and the pipe wall when the two-dimensional model is used. Through observation, it can be found that the mixed oil concentration curve calculated by the two-dimensional model is asymmetric. This asymmetry is caused by uneven flow velocity across the pipe section. The calculation results of the one-dimensional model and the two-dimensional model are different in the length of the mixing section, the position of the mixed oil section, and the degree of symmetry of the mixed oil concentration curve. This is due to the large model error of the one-dimensional model itself, which can only describe the change in the mixed oil concentration in the axial direction of the pipeline. In

Comparison of One-Dimensional Model and Two-Dimensional Model
The solution results of the one-dimensional model and the two-dimensional model are shown in Figure 8. The red dashed line in the figure is the calculation result of the one-dimensional model. After testing, the deviation volume is 0, which is a center-symmetric figure. The blue and black solid lines, respectively, indicate the changes in the concentration of oil at the axis of the pipeline and the pipe wall when the two-dimensional model is used. Through observation, it can be found that the mixed oil concentration curve calculated by the two-dimensional model is asymmetric. This asymmetry is caused by uneven flow velocity across the pipe section. The calculation results of the one-dimensional model and the two-dimensional model are different in the length of the mixing section, the position of the mixed oil section, and the degree of symmetry of the mixed oil concentration curve. This is due to the large model error of the one-dimensional model itself, which can only describe the change in the mixed oil concentration in the axial direction of the pipeline. In the calculation process, the average concentration used in the one-dimensional model is only related to the axial position, and the concentration, average density and average viscosity are considered unchanged on the same section of the pipeline. The two-dimensional model continuously updates the average density and average viscosity at each node by calculating the concentration of each node.   Figure 9 shows the numerical simulation results of the two-dimensional model. At the axis of the pipeline, the replacement speed of oil is faster than that at the pipe wall, because the flow rate of the oil at the center of the pipeline is greater than that at the wall.  Figure 9 shows the numerical simulation results of the two-dimensional model. At the axis of the pipeline, the replacement speed of oil is faster than that at the pipe wall, because the flow rate of the oil at the center of the pipeline is greater than that at the wall.  Figure 9 shows the numerical simulation results of the two-dimensional model. At the axis of the pipeline, the replacement speed of oil is faster than that at the pipe wall, because the flow rate of the oil at the center of the pipeline is greater than that at the wall.
(a) Pipeline axial-radial oil concentration (b) Pipeline axial oil concentration

Comparison of the Effects of Convection and Diffusion on the Mass Transfer Process
In this paper, the time-division method is used to solve the convection-diffusion equation. The convection-diffusion equation is divided into a convection equation and diffusion equation to be solved separately. One of the advantages of this is that it can solve the convection equation and the diffusion equation separately, and compare the influence of the two through the numerical simulation results, and observe the concentration of oil when only the convection equation or only the diffusion equation is calculated. However, in the actual oil-mixing process, convection and diffusion always exist simultaneously. This paper compares the solution results of pure convection

Comparison of the Effects of Convection and Diffusion on the Mass Transfer Process
In this paper, the time-division method is used to solve the convection-diffusion equation. The convection-diffusion equation is divided into a convection equation and diffusion equation to be solved separately. One of the advantages of this is that it can solve the convection equation and the diffusion equation separately, and compare the influence of the two through the numerical simulation results, and observe the concentration of oil when only the convection equation or only the diffusion equation is calculated. However, in the actual oil-mixing process, convection and diffusion always exist simultaneously. This paper compares the solution results of pure convection equation and pure diffusion equation, aiming to compare the strength of the effect of convection on diffusion, and deepen the intuitive understanding of the effects of convection and diffusion.
The left side of Figure 10 shows the simulation results of simultaneous calculation of mixed oil concentration and diffusion, and the right column shows the distribution of oil mixture concentration when only convection terms are calculated. The abscissa in the figure represents the axial direction of the pipeline, and the ordinate represents the radial direction of the pipeline. The length of the pipe used for calculation is 0.8 m and the pipe diameter is 250 mm. The simulated conveying temperature is 20 • C. Run 0# diesel first, then 90# gasoline. The average flow velocity is 1.8 m/s. Figure 10a In the initial stage of oil mixing, the oil mixing speed is faster. After the occurrence of the mixed oil section, the mixed oil section acts as an isolation for the two product oils, and the growth rate of the mixed oil amount slows down. At this flow rate, when the running time is 0.5 s, the mixed oil section has basically left the calculation pipe section. When only convection is calculated, the product oil concentration distribution is basically consistent with the flow velocity distribution at the pipe interface, which is not much different from the product oil concentration distribution when the convection diffusion equation is solved.
the mixed oil amount slows down. At this flow rate, when the running time is 0.5 s, the mixed oil section has basically left the calculation pipe section. When only convection is calculated, the product oil concentration distribution is basically consistent with the flow velocity distribution at the pipe interface, which is not much different from the product oil concentration distribution when the convection diffusion equation is solved.  Figure 11 shows the comparison of the oil concentration distribution between the twodimensional model and only the diffusion term. The left side is the calculation result when the convection term and the diffusion term are calculated at the same time, and the right side is the result when only the diffusion term is calculated. Figure 11a-e corresponds to running times of 0.1 s, 0.2 s, 0.3 s, 0.4 s and 0.5 s, respectively. The speed of diffusion is obviously slower than that of convection. When calculating the diffusion coefficient, the distance between the calculation position and the pipe wall is also considered, and different formulas are used for calculation in different areas. At the  Figure 11 shows the comparison of the oil concentration distribution between the two-dimensional model and only the diffusion term. The left side is the calculation result when the convection term and the diffusion term are calculated at the same time, and the right side is the result when only the diffusion term is calculated. Figure 11a-e corresponds to running times of 0.1 s, 0.2 s, 0.3 s, 0.4 s and 0.5 s, respectively. The speed of diffusion is obviously slower than that of convection. When calculating the diffusion coefficient, the distance between the calculation position and the pipe wall is also considered, and different formulas are used for calculation in different areas. At the bottom layer of laminar flow, the diffusion coefficient is equal to the molecular diffusion coefficient, and a constant is used in the calculation process. In the buffer layer and the turbulent core area, the size of the turbulent diffusion coefficient is positively related to the distance from the pipe wall, so the concentration field under the action of diffusion also shows that the amount of oil mixed at the pipe axis is more than that at the pipe wall.
According to the formula, the diffusion coefficient in the model is only related to the radial position, and the diffusion coefficient at the same radial position is considered the same. Figure 12 shows the distribution of the calculated turbulent diffusion coefficient and molecular diffusion coefficient at different radial positions of the pipe section. In the calculation process, the molecular diffusion coefficient is a fixed value. It can be seen from the figure that the turbulent diffusion coefficient takes the largest value at the pipe axis and the smallest value at the pipe wall. This is also the reason why when only the diffusion term is calculated in Figure 11, the oil spreads faster at the axis of the pipe and slower at the pipe wall.
bottom layer of laminar flow, the diffusion coefficient is equal to the molecular diffusion coefficient, and a constant is used in the calculation process. In the buffer layer and the turbulent core area, the size of the turbulent diffusion coefficient is positively related to the distance from the pipe wall, so the concentration field under the action of diffusion also shows that the amount of oil mixed at the pipe axis is more than that at the pipe wall. According to the formula, the diffusion coefficient in the model is only related to the radial position, and the diffusion coefficient at the same radial position is considered the same. Figure 12 shows the distribution of the calculated turbulent diffusion coefficient and molecular diffusion coefficient at different radial positions of the pipe section. In the calculation process, the molecular diffusion coefficient is a fixed value. It can be seen from the figure that the turbulent diffusion coefficient takes the largest value at the pipe axis and the smallest value at the pipe wall. This is also the reason why when only the diffusion term is calculated in Figure 11, the oil spreads faster at the axis of the pipe and slower at the pipe wall.

The Impact of the Viscous Sublayer
To measure the impact of the viscous sublayer on the asymmetry of mixed oil, and to describe the oil velocity distribution more accurately at the viscous sublayer, the boundary layer grid was dense, as shown in Figure 13

The Impact of the Viscous Sublayer
To measure the impact of the viscous sublayer on the asymmetry of mixed oil, and to describe the oil velocity distribution more accurately at the viscous sublayer, the boundary layer grid was dense, as shown in Figure 13. Among them, (a) is the calculation result when the boundary layer mesh is not dense, (b) is the effect of encrypting the original outermost grid, and (c) is the enlarged display of the red box in (b).

The Impact of the Viscous Sublayer
To measure the impact of the viscous sublayer on the asymmetry of mixed oil, and to describe the oil velocity distribution more accurately at the viscous sublayer, the boundary layer grid was dense, as shown in Figure 13. Among them, (a) is the calculation result when the boundary layer mesh is not dense, (b) is the effect of encrypting the original outermost grid, and (c) is the enlarged display of the red box in (b).  Figure 14 compares the distribution of oil concentration in the pipe section before and after the boundary layer mesh refinement. The left column is the simulation effect when the boundary layer mesh is not refined, and the right column is the refined effect. (a-c) correspond to a running time of 0.2 s, 0.4 s and 0.5 s respectively. Comparing the left and right columns of Figure 14a, in the initial stage of oil mixing, the oil concentration distribution after the boundary layer mesh refinement is roughly similar. Comparing the left and right columns of Figure 14b, as the calculation time increases, the effect of boundary layer mesh refinement gradually manifests. The refinement of the boundary layer mesh not only affects the oil concentration distribution close to the pipe wall, but also affects the oil concentration distribution at the center of the pipeline. Comparing the left and right columns  Figure 14 compares the distribution of oil concentration in the pipe section before and after the boundary layer mesh refinement. The left column is the simulation effect when the boundary layer mesh is not refined, and the right column is the refined effect. (a-c) correspond to a running time of 0.2 s, 0.4 s and 0.5 s respectively. Comparing the left and right columns of Figure 14a, in the initial stage of oil mixing, the oil concentration distribution after the boundary layer mesh refinement is roughly similar. Comparing the left and right columns of Figure 14b, as the calculation time increases, the effect of boundary layer mesh refinement gradually manifests. The refinement of the boundary layer mesh not only affects the oil concentration distribution close to the pipe wall, but also affects the oil concentration distribution at the center of the pipeline. Comparing the left and right columns of Figure 14c, after considering the refinement of the boundary layer mesh, the mixing of oil close to the pipe wall is significantly slower. In the right column of Figure 14c, close to the pipe wall, you can see the mixed oil tail trace left by the front oil.  Figure 14c, after considering the refinement of the boundary layer mesh, the mixing of oil close to the pipe wall is significantly slower. In the right column of Figure 14c, close to the pipe wall, you can see the mixed oil tail trace left by the front oil.

Numerical Analysis of Influencing Factors of Asymmetry Distribution
It can be seen from the analysis of the distribution characteristics of the concentration field that the concentration distribution on the pipeline interface is not uniform. To further investigate the effect of influencing factors on the asymmetry of the concentration curve, we analyze them one by one in this section. The factors are adsorption coefficient, viscosity, diffusion coefficient, temperature, pipe

Numerical Analysis of Influencing Factors of Asymmetry Distribution
It can be seen from the analysis of the distribution characteristics of the concentration field that the concentration distribution on the pipeline interface is not uniform. To further investigate the effect of influencing factors on the asymmetry of the concentration curve, we analyze them one by one in this section. The factors are adsorption coefficient, viscosity, diffusion coefficient, temperature, pipe diameter, velocity. The details of the results and analyses can be found in the Supplementary Materials. Based on the detailed analyses, we further carried out a sensitivity analysis of these factors, as shown in the following four figures.
In Figure 15, when the changing rate of influencing factors is +10%, the changing rate of deviation volume V d 0.5 is: +75.8%, −18.7%, −9.8%, +5.6%, +2.7%. The degree of influence of each factor on the deviation volume V d 0.5 , in descending order, is adsorption coefficient, pipe diameter, average flow velocity, temperature, and diffusion coefficient. Among them, the adsorption coefficient and the thermometer diffusion coefficient have a positive effect, while the pipe diameter and flow velocity have a negative effect. The oil viscosity has non-linear relationship with the deviation volume. When the viscosity is increased, the variation rate of the deviation volume first decreases slightly, and then increases. The sensitivity analysis of the deviation volume Vd 0.1 is shown in Figure 16. When the changing rate of the influencing factor is +10%, the changing rate of the offset volume Vd 0.1 is: −40.0%, −21.4%, +14.1%, +8.5%, +5.7%. Thus, the degree of influence of each factor on the deviation volume Vd 0.1 in descending order is pipe diameter, average flow velocity, adsorption coefficient, temperature, and diffusion coefficient. Among them, the adsorption coefficient and the thermometer diffusion coefficient have a positive effect, while the pipe diameter and flow velocity have a negative effect. The oil viscosity has no linear relationship with the offset volume Vd 0.1 and its influence are small.
The sensitivity analysis of the ratio of the two deviation volumes (Vd 0.1 /Vd 0.5 ) is shown in Figure  17. When the changing rate of the influencing factor is +10%, the rate of change in the ratio of the two deviation volumes is: −13.8%, −9.6%, +3.3%, +3.0%, +2.6%. The change in the influencing factor does not exceed 4% when the changing rate is +10%. The degree of influence of each factor on the ratio of the two in descending order is pipe diameter, average flow velocity, diffusion coefficient, adsorption coefficient, and temperature. Among them, the adsorption coefficient and the thermometer diffusion coefficient have a positive effect, while the pipe diameter and flow velocity have a negative effect. The oil viscosity has no linear relationship with the ratio of the two and its influence is small. The influence of various factors on the ratio of the two is like the influence on the deviation volume Vd 0.1 , but the influence of the contrast value is smaller. The sensitivity analysis of the deviation volume V d 0.1 is shown in Figure 16. When the changing rate of the influencing factor is +10%, the changing rate of the offset volume V d 0.1 is: −40.0%, −21.4%, +14.1%, +8.5%, +5.7%. Thus, the degree of influence of each factor on the deviation volume V d 0.1 in descending order is pipe diameter, average flow velocity, adsorption coefficient, temperature, and diffusion coefficient. Among them, the adsorption coefficient and the thermometer diffusion coefficient have a positive effect, while the pipe diameter and flow velocity have a negative effect. The oil viscosity has no linear relationship with the offset volume V d 0.1 and its influence are small.
The sensitivity analysis of the ratio of the two deviation volumes (V d 0.1 /V d 0.5 ) is shown in Figure 17.
When the changing rate of the influencing factor is +10%, the rate of change in the ratio of the two deviation volumes is: −13.8%, −9.6%, +3.3%, +3.0%, +2.6%. The change in the influencing factor does not exceed 4% when the changing rate is +10%. The degree of influence of each factor on the ratio of the two in descending order is pipe diameter, average flow velocity, diffusion coefficient, adsorption coefficient, and temperature. Among them, the adsorption coefficient and the thermometer diffusion coefficient have a positive effect, while the pipe diameter and flow velocity have a negative effect. The oil viscosity has no linear relationship with the ratio of the two and its influence is small. The influence of various factors on the ratio of the two is like the influence on the deviation volume V d 0.1 , but the influence of the contrast value is smaller.   Figure 18 shows the sensitivity analysis of the sensitivity analysis of the ratio of mixed oil head to mixed oil tail length. When the changing rate of influencing factors is +10%, the changing rate of the ratio (mixed oil head/mixed oil tail length) is: +40.4%, −6.2%, +3.5%, −3.2%, +1.0%, +0.9%. The degree of influence of each factor on the ratio of the two, in descending order, is adsorption coefficient, pipe diameter, temperature, average flow rate, oil viscosity, diffusion coefficient. Among them, it is positively correlated with adsorption coefficient, temperature, oil viscosity and diffusion coefficient, and negatively correlated with pipe diameter and average flow velocity.   Figure 18 shows the sensitivity analysis of the sensitivity analysis of the ratio of mixed oil head to mixed oil tail length. When the changing rate of influencing factors is +10%, the changing rate of the ratio (mixed oil head/mixed oil tail length) is: +40.4%, −6.2%, +3.5%, −3.2%, +1.0%, +0.9%. The degree of influence of each factor on the ratio of the two, in descending order, is adsorption coefficient, pipe diameter, temperature, average flow rate, oil viscosity, diffusion coefficient. Among them, it is positively correlated with adsorption coefficient, temperature, oil viscosity and diffusion coefficient, and negatively correlated with pipe diameter and average flow velocity.  Figure 18 shows the sensitivity analysis of the sensitivity analysis of the ratio of mixed oil head to mixed oil tail length. When the changing rate of influencing factors is +10%, the changing rate of the ratio (mixed oil head/mixed oil tail length) is: +40.4%, −6.2%, +3.5%, −3.2%, +1.0%, +0.9%. The degree of influence of each factor on the ratio of the two, in descending order, is adsorption coefficient, pipe diameter, temperature, average flow rate, oil viscosity, diffusion coefficient. Among them, it is positively correlated with adsorption coefficient, temperature, oil viscosity and diffusion coefficient, and negatively correlated with pipe diameter and average flow velocity.

Conclusions
A mathematical model representing the asymmetry of mixed oil concentration, modified by polar adsorption near pipe wall, is established. The model is validated by experiments and empirical formula. Factors influencing the asymmetry are analyzed in detail. Based on this, a sensitivity analysis is made to quantify the influencing content of each factor. The concept of deviation volume is proposed for the first time to represent the asymmetry of mixed oil. It can be used to quantitatively define the asymmetry of mixed oil concentration curve. Thus, the study on mixed oil concentration is put forward from the qualitative phenomenon to the quantitative analysis. This study is the first report on the local characteristics of mixed oil concentration rather than previous studies to treat the development of mixed oil as a whole, so that more accurate mixed oil cutting in batch transportation of multi-product oils in pipeline engineering has a powerful tool.