This paper investigates the problem of gas-liquid flow SO2 removal in nonquiescent flue gas. Using venturi scrubber as the prototype, the population balance model (PBM) combined with the CFD is implemented to characterize the droplets behaviors. Discrete methods as class model (CM) and various quadrature-based moment models (QBMMs) are applied to numerically solve the population balance equations (PBEs). Taking NaOH solution as the reaction kinetics, the sulfur removal efficiency simulation with CM and different QBMMs methods is validated through the operation measurements. The comparison results show that the CM can achieve better accuracy with more bins, which showed the minimal error 3.6%, consisting 30 bins. However, the computational time of the CM is approximately 19.3 times as long as QBMMS. Among the QBMMs, the ECQMOM approach enjoys the best balance between the simulation efficiency and accuracy, while EQMOM shows the least computational load and CQMOM wins the minimal calculation precision. This result will provide sufficient reference for engineers working in the field of the droplets distribution in the venturi scrubber design.

1. Introduction

As the main pollutant in the petroleum industry, exhaust gas from fuel contains a significant amount of SO2 owing to the sulfur contained in fossil fuels [1]. The emissions of SO2 severely cause atmosphere pollution and ecosystem damage, which significantly threatens the human life. Also, the increase of the emission is one of the principal causes of the acid rain [2]. With the stringent emission standards of SO2, the improvement of currently used methods is essential. Currently, wet flue gas desulfurization coupled with selective catalytic reduction is popular for the simultaneous SO2 removal [3, 4].

The venturi that is investigated in this paper is a device used to absorb sulfur dioxide, which is installed at the end of production line as an equipment to clean the flue gas. Venturi scrubber is a two-phase reactor with low-energy input and good-mixing effect. Now most of the relevant studies are based on experiments, which require a long time and a lot of cost to verify the correctness of the designs. The flow field inside venturi is complex, and the fluid velocity is fast, which make the experiment difficult to observe. These above make numerical simulation indispensable in venturi’s research. So this paper aims at finding the most suitable simulation method of the venturi. Due to the high ratio of gas to liquid, the result of CFD with multiphase model cannot meet the rate measured in factory. Thus, PBM is needed to compute a correct result. And most reactors use low-density medium as the discrete phase to study the dynamics characteristics of particles. Besides, the discrete phase of the venturi selected in this paper is liquid. In this case, the behavior of particles such as growth, breakup, and aggregation is more affected by its own stress, which is different from the situation of gas as the discrete phase. So the paper studies the condition of the internal flow field and adaptability of different discrete methods to the venturi firstly. The condition of the flow can show the process of the mass and energy transfer, which could indicate the direction of redesigns. And the choices of the discrete methods should consider both accuracy and time cost, which could provide the basement with subsequent investigations and designs.

As for energy transfer between two phases, Bhutada et al. [5] investigated the relationship between throat structure and the change of pressure and velocity. Dopkin et al. [6] considered that longer multiphase-mixed tube could make consumption of momentum increased due to friction force. In contrast, shorter tube which cannot meet requirement of complete exchange of energy is also averse to mixture. Sun et al. [7] presented an analysis model to optimize venturi based on constant pressure mixing theory. Galanis et al. [8] proposed a thermodynamic mode to calculate critical pressure and structure parameters of venturi. To study the process of reaction happening in venturi, Guerra et al. [9] thought that there is a strong dependence of the dispersed phase distribution on the parameter values used in the droplet size distribution model. Majid et al. [7] studied gas velocity in throat, volume faction, and removal efficiency to analyze performance of venturi scrubber. As for the multiphase flow model, Lehr et al. [10] fully characterized and described the behavior of the discrete phase. Ishii et al. [11] studied the interaction of discrete phases with surrounding media. Ameras et al. [12] discussed turbulence and mixing caused by discrete phases. Due to a series of studies on multiphase flow, numerical simulations can be used for satisfactory simulation work in many cases. The volumetric mass transfer coefficient is the key parameter in mass transfer model in venturi scrubber. However, there are some deficiencies in the direct prediction of , since its theoretical correlation of the microscale hydrodynamics and microscale interfacial phenomena is very complex.

The mass transfer coefficient is commonly predicted by the penetration model [10] and eddy cell model [11]. depends on the behavior of dispersed phase, such as breakup and aggregation, which will change the area between two phases and the transfer of mass and energy. As shown in Figure 1, the liquid in the venturi reactor atomizes small particles under pressure after injection into the concentration. The particles will grow, break up, and aggregate due to the influences by internal and external factors. The external stress of continuous phase leads to deformation to particles. Nevertheless, the surface stress of particles and viscous stresses result in a stable shape of particles [13]. In multiphase, the particles would grow by adsorbing deposited liquid and spreading itself [14]. At the same time, the collision of droplets and vortex in flow field would make droplets breaking under turbulent conditions. And the viscous stress of continuous phase can cause the velocity gradient at the interface around the particles, which would also lead to deformation and even breakage. As for aggregation, the eddy will cause the collision of particles less than 8 mm in size. The difference in speed on particles and wake eddy of big droplets would also cause aggregation happening (Figure 1).

In the simulation of the multiphase flow, the population balance model (PBM) have a good performance to describe the process of reaction and the Sauter Mean Diameter (SMD) of droplets, which can be used to analyze the behavior of the population of droplets from the single droplet in their local environment [1517]. Danckwerts [18] first used population balance theory to describe the residence time distribution (RTD) of the process of flowing and mixing. Hulburt et al. [19] built PBM matrix equations, used to show the behaviors of particle about nucleating, growing, and aggregating in chemical engineering. Coulaloglou [20] set up phenomenological model about breakup and aggregation of droplets. Lee et al. [21] presented functions of breakup and aggregation, which made probability density distribution replaced by number density distribution. On the basis of this, Lehr et al. [22] investigated the distribution of particle size in bubble-bed reactor and presented aggregation function basing on the experiment. Olmos [23] used bubble-bed reactor to study SMD and gas volume fraction. Venneker et al. [24] researched local gas holdup and volume mass transfer coefficient in stirred tank reactors by PBM.

To simulate the discrete phase distribution, the population balance equation (PBE) is used to describe the particles behaviors as breakage and aggregation. Several numerical methods are commonly applied, such as the class method (CM) [25, 26] and method of moments (MOM) [27, 28]. The CM can descript the description of the particles with a highly computational resource demanding. It needs a number of bins as the amount of discretized sections, which is larger leading to a more accurate model with a longer calculate time. The MOM is able to predict the SMD as accurately as the CM, which cannot drawback the description about diameters of particles. The MOM calculates the PBE by tracking the time dependence of the lower-order moments of the SMD rather than the distribution of particle size. Because of the balance between complexity of calculation and accuracy, the simulation methods that are based on the MOM are very popular, of which the methods that used Gaussian quadrature to close the source terms are grouped into the so-called quadrature-based moment methods (QBMMs).

To extend the area adopted by the MOM, McGraw [29] proposed the quadrature method of moments (QMOM) to solve the problems in simulation by PBE. In the QMOM, the moments of the SMD are approximated by the n-point Gaussian quadrature, which converts the integral form of the number density function to the sum of characteristic abscissas and weights that can be computed by product-difference (PD) algorithm [30]. And the SMD can be approximated by the abscissas with a delta function [30]. However, there are some restrictions in use of the QMOM; for example, it is only applicable to univariate distributions due to the limitations of Gaussian quadrature. As the result, the conditional quadrature method of moments (CQMOM) was came up to solving the multivariate system. It uses conditional probability function theory to calculate the multivariate SMD. In the CQMOM, the velocity of particles is treated as a separate “internal coordinate,” and the quadrature approximation is used to overcome the closure problem, which consists of different abscissas or nodes that can be considered as separate particle classes. This method was mainly used to predict particle trajectory crossing [31]. Another method is the extended quadrature method of moments (EQMOM), which was developed by Yuan et al. [32]. The QBMMs that were referenced above only predict a discontinuous value of SMD, which would lead to problems in the situation where the PBE cannot explain the loss of particles of zero-size well. To evaluate the problem, the value of the particle size for a zero-size particle is required [33]. One possible solution is to increase substantially the number of abscissas so that the phase space is adequately discretized. Besides, the EQMOM was presented to simulate the distribution of the particles from a moment set using continuous kernel density functions, which has heavier computation compared to QMOM and CQMOM [34]. The ECQMOM was presented by Chalons [35] et al., of which calculations of the source terms and the moment transport equations are similar to the CQMOM and the EQMOM.

In this contribution, four type of QBMMs, as QMOM, CQMOM, EQMOM and ECQMOM as well as the CM, have been presented to solve the multiphase reaction with multivariate SMD and mass transfer, which are described in Section 2. In Section 3, the apparatus setup and the numerical strategy are introduced about the proposed venturi scrubber. In Section 4, the experiment and the CFD-PBM simulation results are compared with CM and QBMMs, which the balance between the calculation efficiency and model accuracy are analyzed. Finally, conclusions are drawn in Section 5.

2. Numerical Methods

The sulfur dioxide content before and after the flue gas passes through the venturi scrubber with the NaOH solution injection, where the mass transfer efficiency of wet desulfurization is simulated throughout the venturi scrubber chamber [36]. The CM and QBMMs are adopted in the simulation of this paper to simulate the change of droplets size inside the venturi scrubber, so as to calculate the average specific surface area of droplets and mass transfer efficiency and compare with the measured data to calculate the error and measure the accuracy of the model in this work.

The population balance model is adopted for modeling the droplets behaviors in the venturi scrubber, which describes the processes of the natural growth and decrease of droplets and the birth and death of particles occur due to aggregation and breakage. And the breakup kernel which was presented by Hengel et al. [36] and aggregation kernel which was presented by Vigil and Ziff [37] was applied in the simulation. The governing equations for PBM model [38, 39] are shown in Table 1.

After the procedure given above, the SMD is worked out by the simulation, which is used to calculate the contact surface between gas and liquid phases. Then the removal rate of sulfur dioxide is reckoned based on the mass transfer equation, which are shown in Table 2 [4447].

3. Experiment and Simulation

This section may be divided into subheadings. It should provide a concise and precise description of the experimental results and their interpretation as well as the experimental conclusions that can be drawn.

3.1. Apparatus Setup

In this work, a sulfur removal scrubber with a capacity of 4000 tons/year is taken as the prototype (refer to Figure 2). In venturi, sodium hydroxide solution is used to absorb sulfur dioxide from the flue gas, thus achieving the purpose of cleaning the gas. And at the same time, the liquid transmits some energy to the gas by crashing. The venturi scrubber is composed of gas-liquid two-phase inlet, receiving chamber, throat, and diffusion section. The inlet diameter of gas that the device involved in this work is d1, which is 750 mm. The diameter of the liquid sparger is d2, which is designed by requirement of spray. The convergence has a diameter of dcon, which is 1.3 times as long as d1, with a length of l2 (1.9d1). The throat has a diameter of dthroat (0.85d1), and a length of l4 (0.48d1). The diameter of the outlet of the diffusion section is dout, which is the same length as d1. The contraction angle from the convergence to the throat is θ1, and the expansion angle after the throat is θ2.

In field equipment, SGA94-SO2 flue gas analyzer made by KANE corporation (refer to Figure 3) is used to measure sulfur dioxide concentration at the gas phase inlet and outlet of the scrubber respectively, refer to Figure 3. After the analyzer is calibrated in normal air, the probe of the analyzer is dipped into the venturi through the three holes on the wall at the monitoring locations (inlet, throat, and outlet) marked in Figure 2(a). The sensor measures the concentration of sulfur dioxide through the degree of oxidation reaction on the electrode, and converts it into a current signal and demonstrates the result on the portable device. Field data measured by sensors are shown in Table 3 (refer to Table 3).

3.2. Numerical Strategy

In this section, the model in the previous section is simulated and verified. Three-dimensional transient and steady simulations are carried out with the commercial CFD software Ansys FLUENT 18.0. The iteration time step is set to be 1e − 5 s. In this study, Tomiyama’s model is chosen to describe the drag force [48] between the multiphase. The Universal-drag model is used to the lift force [49]. The wall lubrication force [50] chose Antal et al model. And the turbulent dispersion force [51] chose Burn et al model. The realizable k-epsilon model is added to simulate the turbulence kinetics [52]. Besides, the stagey of convergence adopts a logarithmic descent direction algorithm, which was presented by Wu et al. [53]. The proposed algorithm is based on the Karush–Kuhn–Tucker necessary optimality condition and the damped Newton method. At last, the iteration residual errors are set to be 1e − 5.

On the basis of the simulation, the PBM model is adopted to use the CM. According to Kazakis et al. [54], the mean Sauter diameter of droplets performances at the liquid sparger is initialed as equation (37), and the distribution is compiled in UDF. In the process of CM simulation, different numbers of bins and different ratio exponents are set to change the diameter range and accuracy of droplets calculated by simulation, so as to find the most suitable simulation parameters for the reactors and reaction conditions involved in this paper. The QBMMs are implemented for the case comparison study. The maximum and minimum particle sizes of the inlet droplets are set as 2 mm and 0.05 mm, respectively. And the calculation results of the above model are compared for simulation accuracy of different QBMM models.where d32 is the mean Sauter diameter and ds is the diameter of sparger. And the distribution of particle is calculated as follows:

When the atomized medium is water, a = −0.097, b = 0.392, and h = −0.177.

In the above simulation model, sodium hydroxide is added to absorb sulfur dioxide, and the reaction model proposed by Uchida [55] is used:

As the results of simulation, the concentration of sulfur dioxide is compared with the measure value of experiment to obtain the errors of different models and the simulation model of venturi reactor which is the most suitable for the working environment mentioned in the paper is found by the comprehensive consideration of calculation time and calculation accuracy.

4. Results and Discussion

In this section, the results of above experimental and simulation are introduced, and the results are compared and analyzed to draw relevant conclusions.

4.1. Mesh Size Analyze

Firstly, this paper discusses the influence on the results by mesh size, which compares the simulation results of velocity under different sizes of mesh. In the simulation, more accurate results need the smaller mesh, which would increase the computation and prolong the simulation time. Bhutani et al. firstly applied mesh adaptivity to population balance equation for modeling the multiphase flow and presented the mesh at two-phase mixing area should be smaller [56]. This paper adopts a fixed mesh strategy. Because the process of transfer mainly occurs at the injection area, the mesh of the venturi central area is set to half of the maximum mesh size. When smaller size does not improve accuracy, the size should be chosen. The velocity distribution at the outlet of venturi is studied when the size of mesh is 8, 10, 15, and 20 mm, respectively. It can be seen from the results of simulation that when the size of mesh is 15 mm and 20 mm, the velocity at center of the outlet is larger and the circumference is smaller, which has more errors. When the mesh size is less than 10 mm, there is no significant fluctuation in the simulation results. So, 10 mm is adopted as the diameter of mesh in the subsequent simulation in this paper (refer to Figure 4).

4.2. Velocity Distribution

After determining the mesh size, the process of reaction mass transfer is added on the basis of the cold mode flow field to calculate the concentration of sulfur dioxide, and the PBM is also added to investigate the distribution of droplet size. The velocity and pressure distribution of the flow field in the reactor is shown in Figure 5, which shows the velocity at the top, middle, and bottom parts of the reactor and the pressure distribution on the axis. The change of liquid pressure and velocity reflects the change of liquid energy. It can be seen that on the central axis of the reactor, the liquid entrances convergence at a high speed to atomize into droplets, and the speed gradually decreased during the mixing process with the gas. The speed on the top line is a little faster than the speed on the bottom line due to the gas that flows into convergence from top inlet. In convergence, the pressure of liquid drops slightly due to diffusion and friction with the gas. When the flow passes the contraction section, the pressure decreases due to the narrowing pipe with the liquid accelerating. In this section, the pressure drops nearly 300 kPa, and the velocity changes to 19.8 m/s from 9.5 m/s. According to the Bernoulli equation, nearly 50% of pressure decrease does work to accelerate. Wu et al. [57] thought that reduction in pressure due to friction with the pipe is 40% of that resulting in acceleration. So the remaining 30% of the pressure loss works to transfer energy to the gas. And there is sudden change of pressure decrease and velocity increase at the entrance of the throat because of the geometrical mutation between the tail of convergent and throat entrance, which was also mentioned by Wu et al. [57]. Near the exit of throat, the velocity begins to decrease because of backflow, which is made by negative pressure. In the diffuser, the velocity of the mixed phase decreases gradually and the pressure increase due to the wider pipe. The process of change about energy of gas and liquid can be seen by the change of pressure and velocity from Figure 5. The loss of liquid energy is partly due to the conversion of its energy. As the pipe narrows, the pressure drops with an increased velocity. At the same time, the work done by the liquid in the process of collision with the gas also causes the drop of liquid pressure and the increase of gas energy (refer to Figure 5).

4.3. Phase Volume Faction

Figure 6 shows the concentration distribution of liquid and SO2 in venturi, which can directly reflect the simulation of mixing. It can be seen that the liquid is atomized and vaporized under high pressure after entering venturi and the mixing begins to occur at the interface of gas-liquid. In the contraction section, due to the dramatic change of pressure and velocity, the gas distribution tends to be uniform rapidly. When the liquid accelerates, the break frequency of droplets would increase. It is easier to mix the gas with dispersed droplets. So, there is only a little difference in concentration because of the reaction. In the throat, the fastest speed and the effect of negative pressure make the turbulence intensity stronger and the flow field more complex. And under the influence of the above reasons, the droplets are broken more violently and dispersed more evenly. So the throat is where most of the reaction happens in venturi. Finally, the flow entrances the diffuser. In the process of pressure and velocity change, the velocity around the pipeline decreases fast and the velocity in the center decreases slowly due to the influence of friction on the reactor wall. The difference in velocity would create the difference in pressure, which was called collapse pressure by Zhang [58] et al. Also, the collapse pressure makes the cloud of droplets moving inward, which can be seen from Figure 6(a). On the whole, when the gas enters venturi, the concentration of sulfur dioxide is 2 v%. The reaction is basically balanced when the multiphase is mixed. As a result, the concentration of sulfur dioxide at the outlet is 1.097 v%. According to the simulation results, it can be seen that the reaction mainly occurred in the throat section. After entering the diffuser section, the concentration of sulfur dioxide basically does not change, so the throat should be the focus of venturi reactor research (refer to Figure 6).

4.4. Experimental Verification of Different Methods

In the experiment, because the venturi reactor mentioned in this paper is a working equipment, it is impossible to directly obtain the state of the internal flow field by using PIV and other devices. Therefore, only the measured value of sulfur dioxide can be used to verify the correctness of the simulation. So sulfur dioxide concentrations at different locations in the reactor are measured using the sensors mentioned in Figure 3. During normal operation of venturi, the concentration of SO2 is measured on the inlet, throat, and outlet sections of the venturi scrubber in 10 groups. Each group is measured for 10 minutes, and the average value of the meter reading is calculated within 10 minutes as the measurement result. According to the results, the amounts of mean concentration at inlet, throat, and outlet are 633,550 and 392 mg/Nm3, respectively. It indicates that the removal reaction of sulfur dioxide mainly occurs in the convergence and throat, while the interaction between the two phases in the contraction section is mainly dominated by collision and mixing. In the convergence, the concentration of sulfur dioxide decreases by 132 mg/Nm3, which may be because the two phases have higher concentrations of sulfur dioxide and sodium hydroxide at the beginning of the reactor entry, so the reaction happening is larger. In the contraction section, the concentration of sulfur dioxide dropped by 83 mg/Nm3, but the mixture is more uniform. In the throat, where the reaction is most concentrated due to the degree of mixing and the intensity of turbulence, the concentration of SO2 drops by 158 mg/Nm3. Because of the other component in the gas (e.g., H2S and CO) and the accuracy of sensor, there is a 5% error in the measurement. The removal rate by measurement is 48.76%. And the sulfur dioxide is removed 20.66% during the throat (refer Figure 7).

In this paper, venturi is simulated by CM and QBMMs, respectively. And the results of simulation are compared with those of experiment to obtain the most suitable method (refer to Figure 8).

As for the CM, the quantity of bins would influent on the accuracy and time cost consumption. Bannari [59] et al. studied the relation between bins and simulation about gas hold-up and liquid velocity. And Bannari thought 25 bins would a good choice for simulating gas-liquid flow, because using more than 25 bins gave also good results but the computational efforts became prohibitive. But Laakkonen [46] et al. found that more than 80 bins should be used to minimize errors. In this case, the CM method in 5, 10, 20, and 30 bins are used to simulate the flow in venturi. The results of simulation are shown in Figure 8(a). It can be seen that the more the bins, the more accurate the result of simulation is. But a more accurate result could cost more time, which can be seen from Table 4. When the number of bins is chosen as 30, the error is 3.2% (removal rate is 45.59%). However, the simulation by 30 bins spends nearly 10 times as much as that by 5 bins (error of 14.99%) (refer to Table 4).

As for the QBMMs, the most significant advantage is the fast and accurate potential, which becomes more important as the number of bins used in the CM is increased [60]. The QMOM served as the simplest of QBMMs only need few abscissas, which can change the weights to minimize the committed error [61]. And the error of QMOM is 7.75%. Compared with the QMOM, the EQMOM could be slower than QBMMs [62], because an additional parameter is calculated to reconstruct a continuous number density function, and most of the computational resources are spent on the calculation of the source terms. So the EQMOM spends 1.5 times as much as the QMOM with the error of 4.35%, with is the most arcuate method among QBMMS. The distinctive advantage of CQMOM is the speed of computation, which is suitable for simulating the multiphase system with mass transfer [63, 64]. The time, the CQMOM costs, is eight-tenths of what QMOM does with the error of 6.34%. Finally, the ECQMOM has a balance between accuracy and consumption, of which the result is 43.48%, with an error of 5.28% (refer to Table 5).

5. Conclusion

The SMD of droplets plays a fundamental role in the CFD study of multiphase flow via venturi scrubber under industrial operating conditions. PBM is coupled with CFD simulations to account for dispersed phase behavior. The results of velocity, pressure, and volume faction can illustrate the accuracy of the simulation. And it also can be seen that the CM and QBMMs approaches can predict the gas/liquid flow dynamics, mass transfer, and the SO2 removal rate in the venturi scrubber.

The CM can compute SMD of droplets in venturi and describe the distribution about diameter of droplets. And the accuracy of CM heavily depends on the quantity of bins. There is 3.2% gap between the results of experiment and simulation by the CM with 30 bins, while the simulation by CM with 5 bins has an error of 14.99%.

However, the QBMMs do not pay attention to the process of aggregation and breakage of each droplet, which are effectively able to save the time consumed in the simulation of multiphase, which only costs about 5% of time spent by CM. Among the various QBMMs, the CQMOM computes the result of simulation in a shortest time and the EQMOM adapts venturi most satisfactorily. The ECQMOM can calculate the result of simulation fast under a higher accuracy (5.28%). So the ECQMOM can become the best choice for the optimization design.

The application of QBMMs to simulate the condition of flow in the venturi could be an important evolution of investigation in this work. And the QBMMs methods could calculate an accurate result under a faster speed, especially ECQMOM. But this method cannot calculate the distribution of the particle size. So if you want to study quantitatively the size of particles, such as the size of the daughter particle after breakup or the proportion of small particles in the discrete phase, it could not be used. And the error of the QBMMs can be accepted in engineering projects. However, if a higher accuracy is required in the experimental study, the CM method with a large number of bins would be the only choice. In this paper, the PBM-CFD is proved to predict the flow field in venturi at different methods and provide a basis for optimization. A comparison on cost and accuracy between CM and QBMMs in this paper and the reduction of high gas-liquid ratio flow field would be useful, in order to highlight advantages and drawbacks.

Data Availability

(1) The data used to support the findings of Table 3 have been deposited in the figshare repository (https://figshare.com/articles/work_environment_xlsx/12171741). (2) The data used to support the findings of Figure 4 have been deposited in the figshare repository (https://figshare.com/articles/results_in_differnet_mish_size/12172314). (3) The data used to support the findings of Figure 5 have been deposited in the figshare repository (https://figshare.com/articles/Untitled_Item/12173205). (4) The data used to support the findings of Figure 6 have been deposited in the figshare repository (https://figshare.com/articles/concentration_distribution/12173367). (5) The data used to support the findings of Figure 7 have been deposited in the figshare repository (https://figshare.com/articles/measure/12173361). (6) The data used to support the findings of Figure 8 have been deposited in the figshare repository (https://figshare.com/articles/simulation_results/12173220). (7) The data used to support the findings of Tables 4 and 5 have been deposited in the figshare repository (https://figshare.com/articles/simulation_accuracy_and_cost_xlsx/12173328).

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

All authors have read and agreed to the published version of the manuscript. Shuo Zhang and Wenyue Cui were responsible for numerical model. Shuo Zhang and Wenyue Cui provided the experimental device. Shuo Zhang and Wenyue Cui performed simulation. Shuo Zhang and Wenyue Cui performed validation. Tao Wu and Xiaohang Zhao performed data analysis. Shuo Zhang performed investigation. Wenyue Cui wrote the original draft. Shuo Zhang reviewed and edited the article. Shuo Zhang supervised the study. Wenyue Cui performed project administration. Shuo Zhang was responsible for funding acquisition.


The validation measurements were provided by CNPC Northeast Refining & Chemical Engineering Co. Ltd in a certain CNPC refinery plant CFF boiler. This research was funded by National Natural Science Foundation of China (Grant no. 61803071), China Postdoctoral Science Foundation (Grant nos. 2018M631786 and 2019T120205).