MODELING OF MULTIPHASE FLOW IN AN AIR-COOLING SYSTEM USING THE CFD-FSCA APPROACH

The deposition of crystallized ammonium salt particles presents serious damage to the hydrogenation reactor effluent air cooler (REAC) in petrochemical engineering. The correlation technology process was analyzed for predicting the particle deposition situation in air coolers. Moreover, the particle crystallization and deposition were revealed, and a mathematical model for flow field simulation was established by computational fluid dynamics (CFD). In addition, the fast search clustering algorithm (FSCA) realized in Matlab software was applied to detect tube anomaly. The relatively greater risk of deposition or blockage in the tube area was detected, and this finding is consistent with the distribution of the surface temperature field in air coolers. Meanwhile, good qualitative agreement was found between the distribution rule of particle deposition rate as simulated by CFD and the situation of opened header box on the spot. Therefore, the combined method of CFD simulation and FSCA detection can be utilized for predicting the condition of particle deposition in the air-cooling systems.


INTRODUCTION
The petroleum and chemical industries, which prioritize crude oil refining, are the pillar industries of the national economy in China (Lin et al., 2015;Dion et al., 2012).As an important secondary technical process in refinery, the hydrogenation reactor effluent air cooler (REAC) is one of the most important equipments in hydrogenation units.In recent years, the increasing proportion of inferior crude oil processing issues, such as high sulfur, high acid, and presence of chlorine, have caused failure incidents of tube deposition or constant clogging in REAC systems, leading to many stoppage accidents (Radwan et al., 2015;Jin et al., 2010).Moreover, the mixture of NH 4 Cl, NH 4 HS, FeS, and others from the reaction in the system often causes deposition, blockage, or corrosion in air coolers, and this failure type exhibits evident characteristics of locality and risk (Toba et al., 2011;Popova et al., 2015;Fei et al., 2014).
In recent years, several experimental investigations and numerical simulations have been conducted to study this failure behavior.The presence of daily varied corrosive agents (insoluble salt, chlorides, and sulfidic compounds) was studied in crude oil mixture to corrode and choke the tubes, and the origin of the corrosion damage of air-cooler tubes was identified in oil refineries (Jenabali et al., 2005).The crystal and deposition behavior of ammonium chloride salt (NH 4 Cl) and multiphase flow simulation was Brazilian Journal of Chemical Engineering investigated in refinery hydrocracking REAC by using Aspen software and computational fluid dynamics (CFD) technology (Zhu et al., 2015).The root cause of tube explosion was studied in refinery hydrocracking REACs, and the structure impact on the ammonium salt deposition location and the corrosion behavior of wet ammonium salt was analyzed in REAC tubes (Ou et al., 2013).The experimental and CFD modeling was studied on heat transfer, friction factor, and thermal performance of an air-cooled heat exchanger, which was equipped with three types of tube inserts, namely, butterfly, classic, and jagged twisted tape (Shabanian et al., 2010).
However, the risk of deposition or clogging caused by the existence of ammonium salt crystal particles in the hydrogenation REAC system has yet to be comprehensively evaluated given that only crystallization temperature or the under-deposit corrosion rate has been measured (Shammerl et al., 2010).In recent years, numerous analytical and numerical methods have been developed and adopted for forecasting in various engineering fields (Massinaei et al., 2010;Pourtousi et al., 2015;Shamshirband et al., 2015;Karami et al., 2012;Lemoine et al., 2008;Shamshirband et al., 2013;Lai et al., 2009;Abbas, 2008;Petković et al., 2014).Soft computing methods provide several tools to show nonlinear input-output mapping and solve a variety of nonlinear and complex problems.Among these methods, the clustering algorithm is applied for classifying elements into categories or clusters on the basis of their similarity and can be also used for anomaly detection in industrial applications (Rodriguez et al., 2014;Elhamifar et al., 2013;Birant et al., 2007;Assidjo et al., 2008;Kazemi-Beydokhti et al., 2015).Besides, the combination of clustering algorithm and CFD method is rarely researched.Here, a new approach is attempted for anomaly detection in air cooler tubes, during which the data of average fluid velocity and liquid fraction in every tube was exported from CFD results; then, the soft computing technique was implemented by MATLAB software for anomaly detection with the fast search clustering algorithm (FSCA) in two dimensions.
In this paper, the mathematical model of particle deposition in the REAC system is established.This work includes analyzing the technology correlation process, revealing the formation process of ammonium salt crystallization and deposition, and simulating the inner flow field.Meanwhile, the relatively greater risk of deposition or blockage in tube areas was detected by FSCA, for predicting the particle deposition situation in air coolers.These results provide significant reference in terms of in-service inspection and risk assessment of inferior oil processing in the hydrogenation REAC system.

THEORETICAL Formation processes of ammonium salt particles
As shown in the process flow diagram of a hydrogenation REAC system (Figure 1), crude oil is dehydrated, and then flows into buffer tank with the pressure raised by a high-pressure pump.The raw material and hydrogen (including circulating hydrogen and new hydrogen) mix and then transfer heat with the reaction effluent in the heat exchanger.After the pre-heating process in the furnace, this mixture enters into the hydrogenation reactor (DC101A/B) and the hydrocracking reactor (DC102) in proper order for desulfurization, denitrification, deoxidization, and cracking reaction.The reaction effluent drains out from the bottom of the hydrocracking reactor (DC102), and flows through the heat exchanger (E01), the air coolers (EC101) and the heat exchanger (E02) in turn, and then it is separated into circulating hydrogen (reused in the loop), oil and sour liquid through the high pressure separator (FA103).Meanwhile, liquid injection points A and B are set before and after the heat exchanger (E01) in upstream of air coolers, and the temperatures in inlet and outlet are 140ºC and 52ºC, respectively, and the operating pressure is 15.8 MPa.
Crude oil containing N and Cl will react with H 2 , and then produce NH 3 and HCl.The temperature of air coolers will decrease under the action of a fan, and the NH 4 Cl crystallized particles can be generated with the entry of reaction effluent according to the formula (1). (1) In general, the temperature of a reaction effluent in tubes undergoes a continuous decrease.When the temperature decreases to below the NH 4 Cl crystallization temperature, NH 4 Cl will go through a sustained reversible crystallization reaction.The value of K p (product of NH 3 and HCl partial pressure) in the gas phase obeys the change trend of the equilibrium curve (Ou et al., 2011).
As shown in Figure 2, if the temperature decreases by ΔT, the K p value exceeds the corresponding equilibrium constant and NH 4 Cl particles start to crystallize and separate out.Then, the K p value in the gas phase rapidly returns to the position of the crystallization equilibrium curve after NH 4 Cl crystallization is finished, and this system restores the equilibrium process.If the temperature decreases " + by again, n processes of crystallization are repeated as such.When lim(ΔT)→0, the K p value changes from the supersaturation position in the gas phase to the equilibrium position; therefore, the reduction of HCl or NH 3 in this process is correlated with the reaction rate of NH 4 Cl salt.The total amount of NH 4 Cl crystallization can be expressed as the accumulation of each step.Based on the above principle, the amount of NH 4 Cl crystallization Q can be obtained by formula (2). (2) (3) where K the equilibrium coefficient of NH 4 Cl and R is the ideal gas constant (= 8.31×10 -3 kJ/(mol.K)).All values were obtained by calculation (Jin et al., 2015).Therefore, the particle crystallization rate over a certain time period can be acquired according to the above equations.

Governing equations and CFD modeling
The multiphase flow medium is the gas, oil and water phase in the hydrogenation air cooler, so the Mixture model is chosen to solve the equations.The main phase is the gas phase, and the second phase is the oil and water phase.
The continuity equation of the mixed phase is (Safaei et al., 2016;Akbari et al., 2016): The momentum equation of the mixed phase is: (5) where n is the number of phase; F is the body force; µ m is the dynamic viscosity of the mixture; v dr,l is the drift velocity of k phase; ρ l is the density of l phase; φ l is the volume fraction of l phase; σ s is the additional stress modification.
The RNG k -ε turbulent model is adopted to solve the internal flow field parameters (Yakhot et al., 1986).
Turbulent kinetic energy k transport equation: Turbulent dissipation rate ε transport equation: where the turbulent viscosity µ t is a function of turbulent kinetic energy k and turbulence dissipation rate ε; In addition, µ t = ρC µ k 2 / ε /, C 1 = max{0.43,η / (η+5) }, and η = Sk / ε.The model constant values C 1ε , C 2 , σ k , σ ε are set as 1.44, 1.9, 1.0, and 1.2, respectively (Yakhot et al., 1986).S is the source item.The equation of motion for the solid phase is given by: where F p (u-u p ) is the drag force per unit particle mass; g(ρ p -ρ)/ ρ p is gravity source; u p is velocity of the particle phase; F other is the other force (not be considered in this research).F p is defined as formula (9).( 9) where d p is the particle diameter; C D is the drag coefficient; Re is the Reynolds number, defined as (Morsi et al., 1972): (10) The deposition rate of solid particles is defined as: (1 1) where m p is the particle mass; A face is the area of the cell face at the wall, N particles is the total amount of particles.
The half-structural diagram of a hydrogenation air cooler is shown in Figure 3.The whole structure can be obtained with a symmetrical arrangement and contains four parts, namely, inlet bend, header box, distribution plate, and tubes.A total of 24 tubes with length 30 times bigger than its diameter (25mm) were placed in a single row to ensure a stable flow state.The numerical orders of tubes from left to right along the direction of X-axis are T1 to T24, and the rows from up to down along the negative direction of the Y-axis are R1 to R3. Multiphase flow medium, such as oil, vapor, and liquid, enters into the header box through an inlet bend, and a rectangular distribution plate is set directly facing to the outlet of the bend for uniform flow in the header box and tubes.
The simulation domain is built on the basis of the above structure features.As shown in Figure 4, the structured mesh arrangement with hexahedral elements is utilized in the inlet bend and tubes, and the tetrahedral unstructured mesh is used in the header box and distribution plate.In this work, the grid independency study shows that a larger number of grids will not influence the accuracy of this case when the grid quantity is divided into 0.92 million, 1.16 million, 1.23 million, and 1.45 million, and the relative flow parameters, such as flow velocity or pressure drop, do not change evidently.The grid number used in this case is approximately 1.16 million through integration taking into account calculation precision and time.In this case, the Reynolds number is 8.4×10 5 as the gas volume fraction is approximately 93%, and the RNG k -ε turbulent model is adopted to solve the internal flow field parameters.The computational domain is set as velocity inlet and pressure outlet.The shape of solid particles is assumed to be spheres, which are injected uniformly at the inlet at the same speed as the continuous phase.The roughness constant is set to the default value of 0.5.The physical parameters of inlet multiphase flow are shown in Table 1.
The software of ANSYS Fluent is used for numerical simulations.This process was divided into two steps: first, the Euler method was adopted for the continuous phase, and then particle trajectories were tracked by the Lagrange approach.In this work, the SIMPLE algorithm (Cui et al., 2014) was employed to solve the coupling of the pressure and velocity field.Standard discretization schemes were used for the pressure terms, and the second-order upwind discretization schemes were used for the convection terms and divergence terms.The convergence criteria for all steady simulations were set to 10 −5 for each equation.The standard wall function was set in combination with the previously discussed RNG kε model.The number of continuous phase iterations per discrete phase model (DPM) iteration was set to 10 ( Rahmanian et al., 2014;Goodarzi et al., 2014 andChen et al., 2016).A total of 24733 particles would be tracked in this simulation, and the density of particles was set as 1550 kg/m 3 .

Fast Search Clustering Algorithm (FSCA) Modeling
The FSCA is a one-time clustering algorithm based on the calculation of the sample distribution density and the sample distance.The algorithm assumes N no-label sample points and the local density of every sample point j is defined as: where d ji is the distance of samples and d c is the cutoff distance with characteristic robustness with regard to large data; then, the distance for each sample j is built as: The δ j value is determined by computing the minimum distance between point j and any other point with higher density.Thus, the ρ j and δ j of every sample point can be obtained.The point with large ρ j and δ j can be regarded as the clustering center, whereas the point with small ρ j and large δ j can be considered as an abnormal point, which does not belong to any of the clustering centers.The diagram of the algorithm process is shown in Figure 5.
The clustering process of the FSCA in two dimensions (Rodriguez et al., 2014) is shown in Figure 6.According to the above principle, results of the anomaly detection show that points 1 and 10 can be taken as the clustering centers, and points 26, 27 and 28 can be regarded as abnormal points.
The complexity analysis of the FSCA is described as follows.The algorithm assumes N sample points, and n-dimensional features in each data such as the calculation complexity of the distance of two points is n 2 .Every data need to be calculated for the distance with the other data; thus, the complexity is N×n 2 .The  Brazilian Journal of Chemical Engineering calculation of sample local density ρ j is completed by N×(N−1) times cutoff operation; then the minimum distance δ j was calculated.The local density and distance need to be subjected to linear scanning and comparison under the worst situation, and it is recorded as a 2N times' operation.Therefore, the calculation complexity of the overall algorithm can be defined as N×n 2 + N×(N − 1) + 2N = N×(n 2 + N + 1).
The only parameter determined by a person is the cutoff distance d c in the FSCA.In general, the maximum distance between samples is calculated first, and then a percentage reduction is conducted according to this distance.Meanwhile, the percentage of abnormal points in the overall samples is determined by combining the actual situation and manual experience.Finally, the number of abnormal points can be determined.The cutoff distance d c is usually taken as the maximum distance D max from 10 % to 30%.In the next section, the comparison of abnormal points is conducted by considering 10%, 20%, and 30%, respectively, and combining with CFD results to determine the final abnormal points.

The distribution of continuous phase
The numerical simulation of the inner flow field in air coolers was conducted, and the velocity magnitude distribution of the horizontal center section in each row of tubes is shown in Figure 7.The plane of y = 54 mm represent the horizontal center section of the first row (R1), and the planes of y = 0 mm and y = −54 mm represent the second and the third row, respectively.As shown in Figure 7(a), the velocity is basically symmetrically distributed, and the multiphase flow velocity from the center to both ends decreases gradually owing to the increase in flow space in the header box.As a rectangular distribution plate is set directly facing to the outlet of the bend, the structure mutation intensifies the turbulent state of this area, inducing severe mass and momentum exchange of the continuous phase.A symmetric distribution of a high velocity area (zone A) forms near the distribution plate, where the maximum velocity reaches 7.5 m/s.When the multiphase flow reaches both ends of the header box, the velocity decreases to below 2.0 m/s.The velocity in the leeward side of the distribution plate is 1.0-2.5 m/s because of the blocking effect, whereas that of the other tube area is 3.0-4.0m/s.In Figure 7(a), 7(b), and 7(c), the velocity distribution trend is basically consistent in three rows.The velocity vector distribution of the horizontal center section in R2 is shown in Figure 8.The biased flow phenomenon appears after multiphase fluid impact to the distribution plate, and a symmetric distribution of high velocity area and two circinate recirculation zones form on both sides of the distribution plate, making a portion of fluid diffuse along the terminal direction of the header box, and the other fluid continues to flow into tubes.
Figure 9 shows the velocity distribution in the longitudinal center section of the header box (plane of z = 0 mm), and the delaminating phenomenon of fluid velocity forms due to the action of the distribution plate, which acts as flow guide, thus making the whole fluid expand in an orderly manner along the terminal direction of the header box.
The average liquid-phase fraction distribution curve in every tube (Figure 10) presents a shape of "W".Meanwhile, the relatively low liquid-phase fraction distribution in the corresponding tube area

The distribution of solid particles
Statistical analysis of the amount of solid particles in tubes is shown in Figure 12, and a total of 24733 particles was tracked.The amount of solid particles seems to be relatively symmetric in R1 and R2, whereas the distribution of particles in R3 is asymmetric because of the non-uniform distribution of the particles under the carrying action of the liquid phase and the bypass effect of the distribution plate.The multiphase fluid in R2 and R3 exceeds that in R1 because of the effect of gravity; thus, a relatively small amount of solid particles in R1is formed.If the captured amount of solid particles in R2 and R3 in tubes is little, the concentration of NH 4 HS crystalline particle would be low, resulting in low risk of deposition or blockage in the area and normal flow of liquid in tubes.With a large captured amount of solid particles in R2 and (the leeward side of the distribution plate) is also in a relatively low level of velocity magnitude compared with Figure 7 because of the blocking effect caused by the distribution plate.In conclusion, a relatively high risk of ammonium salt crystallization deposition or blockage is present in the tube areas of low fluid velocity and liquid fraction.
In general, the distribution situation of fluid velocity and liquid fraction has a great influence on ammonium salt crystallization deposition or blockage.The data of average fluid velocity and liquid fraction from the CFD simulation in every tube was exported from Fluent, and then the MATLAB software was implemented for anomaly detection by the FSCA.The average fluid velocity and liquid fraction were taken as two input variables, and the minimum distance δ was taken as the output variable; thus, the FSCA was two-dimensional with 72 total data.The comparison with different cutoff percentages in tubes is shown in Figure 11; the detection result shows that the area of T12-T15 is relatively abnormal compared with the other tubes.Meanwhile, the anomaly detection result with 10% cutoff corresponds to the surface temperature distribution of air coolers on the spot (Figure 13).R3, as shown in the areas of T6-T8 and T11-T15 in Figure12, the deposition rate is larger in this tube area compared with the others (Figure 14) and increases the risk of deposition or blockage.As an explanation, the low liquid fraction or fluid velocity cannot dissolve crystalline particles or induce normal flow in tubes; then, more ammonium salt will separate as the wall temperature of tubes decreases gradually under the fan action; finally deposition or blockage in tubes forms.As seen through comparison and analysis, the particle distribution situation in tubes not only basically agrees with the surface temperature field (Figure 13) measured by the infrared camera, but also with the distribution rule of solid particle deposition rate (Figure 14).The   The distribution of solid-phase particles deposition rate by the CFD simulation is shown in Figure 14.The particle deposition rate is higher in the bottom of the header box, and a serious deposition area is observed at the end of the header box, and finally, the particles deposit on the wall.The tube areas with low liquid fraction or fluid velocity are at higher risk of deposition or blockage; moreover, if particle deposition occurs in the junction of header box and tubes, the phenomenon of tube blockage will become increasingly serious until complete blockage occurs.The morphology of corrosion product deposition in the opened header box of air coolers is also shown in Figure 14.The main types of corrosion products after assaying are ammonium salt and FeS, and the deposition position occurs mainly at the end of the header box on the spot, basically in agreement with the distribution rule of solid particle deposition rate as simulated by CFD in the header box.

CONCLUSIONS
Based on the correlation technology process of the hydrogenation REAC system, ammonium salt crystallization and deposition in air coolers was revealed, and the inner flow field was simulated by using the mixture model and the discrete phase model.When a multiphase medium flows through the rectangular distribution plate, the phenomenon of biased flow appears and makes the distribution of flow velocity and liquid phase fraction unbalanced.Then, the MATLAB software was used for detecting the anomaly of tubes by the FSCA in two dimensions.In agreement with the distribution of the surface temperature field in air coolers, a relatively higher risk of deposition or blockage in the tube area is detected.The deposition rate of particles is higher at the bottom of the header box relative to other areas, and serious deposition areas are observed at the ends of the header box.If particle deposition occurs in the junction of the header box and tubes, the phenomenon of tube blockage will become increasingly serious until complete blockage occurs.With a large captured amount of solid particles in R2 and R3, as shown in the areas of T6-T8 and T11-T15, the deposition rate is larger in this tube area compared with the others and increases the risk of deposition or blockage.Good qualitative agreement was observed between the distribution rule of solid particles deposition rate as simulated by CFD and the deposition situation of an opened header box on the spot.These results can verify the prediction result of particle deposition by using the CFD-FSCA approach in air-cooling systems.

ACKNOWLEDGEMENTS
The project is supported by the National Natural Science Foundation of China( 51876194

Figure 1 .
Figure 1.Process flow diagram of a hydrogenation REAC system.

Figure 2 .Figure 3 .
Figure 2. NH 4 Cl crystallization equilibrium curve and sketch of the crystallization process.

Figure 4 .
Figure 4. Sketch of computational grids in hydrogenation air coolers.

Figure 5 .Figure 6 .
Figure 5.The diagram of the algorithm process

Figure 7 .
Figure 7. CFD result of the velocity magnitude distribution in the horizontal center section of each row.

Figure 8 .
Figure 8. CFD result of velocity vector distribution in the horizontal center section of R2.

Figure 9 .
Figure 9. CFD result of velocity magnitude distribution in the longitudinal center section of the header box

Figure 12 .
Figure 12.CFD result of percentage of particle amount distributions in tubes.

Figure 13 .
Figure 13.The surface temperature distribution of air coolers on the spot.

Figure 14 .
Figure 14.The comparison of particle deposition situation in the header box (left: CFD result right: on the spot)

Table 1 .
Physical parameters of inlet multiphase flow.