A Novel Method for the Prediction of Erosion Evolution Process Based on Dynamic Mesh and Its Applications

Abstract

Particle erosion is a commonly occurring phenomenon, and it plays a significantly important role in service life. However, few simulations have replicated erosion, especially the detailed evolution process. To address this complex issue, a new method for establishing the solution of the erosion evolution process was developed. The approach is introduced with the erosion model and the dynamic mesh. The erosion model was applied to estimate the material removal of erosion, and the dynamic mesh technology was used to demonstrate the surface profile of erosion. Then, this method was applied to solve a typical case—the erosion surface deformation and the expiry period of an economizer bank in coal-fired power plants. The mathematical models were set up, including gas motion, particle motion, particle-wall collision, and erosion. Such models were solved by computational fluid dynamics (CFD) software (ANSYS FLUENT), which describes the evolution process of erosion based on the dynamic mesh. The results indicate that: (1) the prediction of the erosion profile calculated by the dynamic mesh is in good agreement with that on-site; (2) the global/local erosion loss and the maximum erosion depth is linearly related to the working time at the earlier stage, but the growth of the maximum erosion depth slows down gradually in the later stage; (3) the reason for slowing down is that the collision point trajectory moves along the increasing direction of the absolute value of θ as time increases; and (4) the expiry period is shortened as the ash diameter increases.

1. Introduction

Erosion is mechanical damage resulting from the impact of particles carried by fluid [1]. When a particle with a certain speed strikes a solid surface, the impact region of the surface will be deformed [2,3,4]. Subsequently, the surface deformation translates into material removal, which shortens the service life, such as the SCR (selective catalytic reduction) catalyst and the economizer in Figure 1, among others. Erosion is a complex problem caused by numerous factors, such as particle velocity, impact angle, particle size, and particle shape.

Considerable studies have been performed on erosion. Finnie [6,7] and Bitter et al. [8,9] presented the earliest work focusing on theoretical erosion mechanism. Finnie’s erosion model was based on the experimental phenomenon of erosion and the models of the micro-cutting mechanism. However, at the impact angle of 90°, it predicted no erosion rate at all, and the analyzed data clearly indicated that it is incorrect. Bitter’s erosion model was proposed with respect to deformation and cutting wear, despite lacking the support of any physical models. Huang et al. [10] addressed these shortages, deriving a phenomenological erosion model by analyzing the normal and tangential forces acting on the abrasive particle. Huang’s erosion model has been increasingly acknowledged by researchers. In addition to the theoretical erosion model, numerous erosion equations have been developed based on experimental tests. Most of these erosion equations have been presented regarding the function of the velocity exponent and the impact angle. Grant and Tabakoff [11] developed an empirical erosion equation. This model pointed out that the coefficient of restitution should be incorporated into the equation, owing to multiple times of impingement by the particle. Okal et al. [12,13] conducted erosion tests for a wide range of materials and erodent particles, considering both material hardness and the load relaxation ratio. Many parameters were involved in this model, although the parameters are hard to acquire. Besides these, Haugen [14], Ahlert [15], Mclaury [16], and Zhang [17] at E/CRC of the University of Tulsa also proposed different empirical erosions, and Shamshirband et al. [18] adapted the ANFIS (adaptive neuro-fuzzy inference system, a type of neural network) model to precisely predict the total and maximum erosion rate.

The aforementioned work focused on the erosion ratio. However, the evolution process of erosion has been little studied in the research, mainly because no methods or models have described long periods of erosion or shape evolution, except for the experiments. However, such experiments take too much time and material power. If feasible methods or models were available to describe the evolution process, the research would make progress.

In this study, the material removal of erosion was investigated by CFD (computational fluid dynamics software) coupled with the UDF (user defined function) of the particle motion and erosion model. It is worth highlighting that the detailed evolution process of erosion will be realized by the dynamic mesh technique. This DPM (discrete phase model) approach, coupling the erosion model and dynamic mesh, combines the material removal and the mesh deformation quantitatively, which can demonstrate the node displacement of the mesh in the erosion region and the new flow field. Accordingly, the flow field and the particle trajectory are transformed at all times. Based on those concepts, the evolution of the erosion loss and erosion depth was investigated by the CFD-DPM approach firstly. The erosion profile was achieved secondly.

Moreover, this paper applies the method to discuss the expiry period of the economizer in modern coal-fired power plants. The economizer is used to improve the overall thermal efficiency in the gas-water exchangers, exacting residual heat energy from the flue gas and transferring the energy to the feed water [19]. However, the flue gas includes ash particles, corrosive gas, temperature oscillation, and so on. Those factors respectively contribute to ash erosion, corrosion fatigue, thermal fatigue, and so forth, which may ultimately lead to bursting of economizer tubes [20]. Among these, ash erosion plays a significantly important role in the expiry period of the economizer. Therefore, this work will benefit the design of the economizer and guide the boiler operation.

2. Method Description

Figure 2 demonstrates the calculative process of erosion evolution. The process includes the flow calculation, particle tracking, material removal calculation coupled with the UDF of the erosion model, erosion profile calculation with the UDF of the dynamic mesh, and the ALE (arbitrary Lagrange–Euler) calculation. ANSYS FLUENT calculates the trajectory of the discrete particle by integrating the force balance on the particle, which is done in the Lagrangian frame. Then, the collision point of the particle-wall collision is tracked, and the particle-wall collision is translated into the material removal in the UDF of the erosion model. Then, the material removal is converted to wall-mesh deformation in the UDF of the dynamic mesh. Dynamic mesh is universally applied for the fluid–structure interaction [21], flapping flight [22] and so forth. However, this dynamic mesh technology is rarely introduced to the erosion profile that is formed by wall-mesh deformation. Figure 3 shows the mesh deformation before and after the particle-wall collision by the dynamic mesh. The value of the material removal is stored in the center of the boundary mesh. However, the dynamic mesh moves the node of the boundary mesh, and the value in the node is interpolated by the value in the center. In this interpolation, this paper keeps the elementary volume constant, as follows:

$$h_{i,j}=\frac{1/4\Delta t{\displaystyle \sum _{k=0}^1{\displaystyle \sum _{l=0}^1{Q}_{i-1/2+k,j-1/2+l}^n}}}{1/4{\rho }_m{\displaystyle \sum _{k=0}^1{\displaystyle \sum _{l=0}^1{A}_{i-1/2+k,j-1/2+l}^{n-1}}}}$$

(1)

$Q_{i-1/2,j-1/2}$ is the rate of the material removal, the function of the particle mass flow rate $m_p$ and the erosion rate $E_r$ in the [$i-1/2$, $j-1/2$] element mesh; $A_{i-1/2,j-1/2}$ is the area of the [$i-1/2$, $j-1/2$] element mesh; n is the time step; $\Delta t$ is the time step size.

The ALE calculation is used to correct the flux after the mesh deformation. In the ALE description, the nodes of the computational mesh may be moved with the continuum in a normal Lagrangian fashion, or be held fixed in a Eulerian manner [23,24]. At last, this paper adopts the unsteady calculation, and the erosion effect of the time step amounts to 1day. It is possible to regard the erosion rate $E_r$ as the constant during 1day. Additionally, when computation time t is at the maximum time t_max (2 years in this paper) given, the calculation is terminated.

When the erosion profile is transformed, the flow field, especially the boundary layer, will be changed also. The specific mathematical relation between the material removal and the wall-mesh deformation is in Formula (1). The wall-mesh deformation leads conversely to the modification of the particle trajectory. This is the coupled calculation of multi-physics, including the flow, particle motion, and erosion.

3. Application

3.1. Mathematical Model

3.1.1. Flow Configuration

The economizer of a 600MW coal-fired power plant was studied in this work. It is composed of 135 rows of S-shaped circular tubes and made from SA 210 GRA1(N); the other detailed parameters are shown in

Table 1. Structural parameters of the economizer tube [25].

Material	Rows	Diameter/mm	Transversal Pitch a/mm	Longitudinal Pitch b1/mm	Longitudinal Pitch b2/mm
SA 210 GrA1(N)	135	Φ51×6	144	102	69

. Considering the periodic distribution of the tube bank in Figure 4, this paper takes the shaded region as the research unit.

Figure 5 gives the three-dimensional diagram of the research unit and the boundary conditions. Boundary conditions were set as follows: velocity-inlet, pressure-outlet, and wall. As mentioned earlier, there exists periodicity in the Y direction. Meanwhile, along the tube (Z direction), the flow is similar. Thus, the periodic boundary is loaded onto the Y direction and the Z direction.

The velocity is selected at 8.12 m/s, which represents 100% of the boiler THA (turbine heat acceptance) load. The diameters of the particles are selected at 50 μm, 100 μm, and 250 μm, with the same mass flow. The mass flow of ash is set to 1.52 × 10⁻² kg/s, equal to 32 g/m³ of the gas.

Considering the effect of the boundary layer on the flow computation, the grid of the tube wall is refined in Figure 6 (enlarged from the red frame of Figure 5). Then it will capture the internal flow field information on the boundary layer. Moreover, the grid-independency study gives the appropriate grid number, whose mean error is controlled at less than 2%. The grid number ensures the accuracy and rapidity of the calculation.

3.1.2. Governing Equations

Gas Motion Model

Gas from the furnace flows around the tube bank of the economizer. The flow field changes, but it meets Navier–Stokes equations. Considering the flow as a turbulent flow, the two-equation turbulent model (standard k-ε model) is introduced. The Navier–Stokes equations coupled with the standard k-ε turbulence model is employed to solve the gas motion. Navier–Stokes equations of gas are as follows [26,27]:

$$\frac{\partial ({\rho }_f{u}_i)}{\partial x_i}=0$$

(2)

$$\frac{\partial ({\rho }_f{u}_i{u}_j)}{\partial x_j}={\rho }_{f}g-\frac{\partial p}{\partial x_i}+\frac{\partial {\tau }_{ij}}{\partial x_j}$$

(3)

$$\frac{\partial \left({\rho }_{f}e{u}_i\right)}{\partial x_i}=-\frac{\partial \left(u_{i}p\right)}{\partial x_i}+\frac{\partial }{\partial x_i}\left(k\frac{\partial T}{\partial x_i}\right)+\frac{\partial \left(u{\tau }_{ij}\right)}{\partial x_i}$$

(4)

${\tau }_{ij}$ is the stress tensor, excluding the surface pressure. The stress tensor is calculated as:

$${\tau }_{ij}=\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\frac{2}{3}\mu \frac{\partial u_k}{\partial x_k}{\delta }_{ij}$$

(5)

$e$ is the total energy, consisting of the internal and kinetic energy. The total energy is calculated as:

$$e=C_{p}T+\frac{1}{2}{V}_f^2$$

(6)

The standard k-ε turbulence model is as follows [28]:

$$\begin{array}{cc}\hfill \frac{\partial \left(\rho k_f{u}_i\right)}{\partial x_i}=& \frac{\partial }{\partial x_j}\left(\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right)+G_k+G_b-\rho \epsilon \hfill \\ & -Y_M+S_k\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \frac{\partial \left(\rho \epsilon u_i\right)}{\partial x_i}=& \frac{\partial }{\partial x_j}\left(\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right)\hfill \\ & +C_{1\epsilon }\frac{\epsilon }{k_f}\left(G_k+C_{3\epsilon }{G}_b\right)-C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}+S_{\epsilon }\hfill \end{array}$$

(8)

where:

$${\mu }_t=\rho C_{\mu }\frac{k_f^2}{\epsilon }$$

(9)

Lagrangian Formulation for Particle Motion Model

The particle motion is traced by Newton’s Second Law in Formula (10), driven by the drag force $F_D$, thermophoretic force $F_T$, gravity $F_G$, and Saffman lift force $F_{SL}$ [29].

$$\frac{d{u}_p}{dt}=F_D+F_T+F_G+F_{SL}$$

(10)

where,

$$F_D=\frac{18\mu }{{\rho }_p{d}_p^2}\frac{C_{D}Re}{24}\left(u_f-u_p\right)$$

(11)

$$F_T=-\frac{6\pi d_p{\mu }^2{C}_S\left(k_f/k_p+C_{t}Kn\right)}{\rho \left(1+3C_{m}Kn\right)\left(1+2k_f/k_p+2C_{t}Kn\right)}\frac{1}{m_{p}T}\nabla T$$

(12)

$$F_G=\frac{g\left({\rho }_p-{\rho }_f\right)}{{\rho }_p}$$

(13)

$$F_{SL}=\frac{2K{v}^{1/2}\rho d_{ij}}{{\rho }_p{d}_p{\left(d_{lk}{d}_{kl}\right)}^{1/4}}\left(u_f-u_p\right)$$

(14)

Particle-Wall Collision Model

In particle motion, the particle–particle collision and the particle-wall collision change the motion state, including velocity and direction of motion. In consideration of the economizer bank in the tail shaft flue and low concentration of ash, there are relatively few possibilities of particle–particle collision. This research mainly focuses on the particle-wall collision, neglecting the particle–particle collisions. Therefore, this study adopted the DPM approach instead of the DEM (discrete element method) approach. In the particle-wall collision, the normal and tangential components of reflected velocity are achieved by the coefficient of velocity restitution. The coefficient of velocity restitution is caused by kinetic energy loss, a function of the incidence angle $\alpha$ in the particle-wall collision model, as follows [30]:

$$e_n=0.993-1.76\alpha +1.56{\alpha }^2-0.49{\alpha }^3$$

(15)

$$e_t=0.998-1.66\alpha +2.11{\alpha }^2-0.76{\alpha }^3$$

(16)

Erosion Model

Particle-wall collision not only changes the motion state, but also removes the material. Huang’s erosion model was adopted to estimate the abrasion loss from the following Formula (17). Huang’s erosion model includes mechanisms of deformation damage removal and cutting removal [10]. This model also demonstrates the effect of particle size on erosion. The erosion rate $E_r$ is the ratio of abrasion loss and particle mass.

$$E_r=C{\rho }_p^{0.15}{\left(V_p\sin \alpha \right)}^{2.3}+D{m}_p^{0.1875}{d}_p^{-0.0625}{V}_p^{2.375}{\left(\cos \alpha \right)}^2{\left(\sin \alpha \right)}^{0.375}$$

(17)

Numerical Procedures

The spatial discretization applies the second-order upwind scheme to discretize the whole convective terms. In addition, the numerical calculation adopts the Eulerian–Lagrangian frame, and the SIMPLEC algorithm (a semi-implicit algorithm) is applied to pressure/velocity coupling in the Eulerian frame. This paper is divided into two numerical procedures, as follows:

• Step 1: Verification of erosion model on SA210 GrA1(N);
• Step 2: Erosion calculation of economizer unit.

3.2. Verification of Erosion Model

C and D are related to the properties of the particle and economizer material in Formula (17). C and D are determined by fitting the curve. By comparing with the erosion experimental data by the ASTM standard-tested bed in the India National Institute of Technology [31,32], the accuracy of Huang’s erosion model was verified and found to be appropriate for SA 210GrA1(N) in Figure 7 and Figure 8. Figure 7 gives the relation between the incidence angle and the erosion rate $E_r$. Compared to Finnie’s erosion model, Huang’s erosion model can not only estimate the effect of the small impact angle, but also the large angle. Figure 8 shows the relation between the particle diameter $d_p$ and the erosion rate $E_r$. From Figure 7 and Figure 8, the predicted results obtained are consistent and in good agreement with experimental results reported elsewhere. For the erosion curve of the impact angle or the ash size, the mean relative error is less than 5%, which satisfies the research requirements.

4. Results and Discussion

4.1. Erosion Loss

4.1.1. Global Erosion Loss

Figure 9 shows the global erosion loss under different ash sizes. From Figure 9, the global erosion loss is linearly related to the working time. Meanwhile, as the ash size increases, the growth amplification of global erosion loss decreases ($\Delta 1>\Delta 2$). This coincides with the fact that it is hard to change the motion of the larger particle, which has the larger inertia. In other words, the motion of the larger particle becomes smaller, and the collision point of a single particle from the injection changes little. Thus, the mass flow of ash impacting the economizer increases slowly and becomes stable gradually, which contributes to the slow growth amplification. In addition, the particle rebound or collision also impacts the global erosion loss, illustrated in Section 4.1.2.

4.1.2. Local Erosion Loss

Figure 10 shows the local erosion loss of different rows under various ash sizes. From this figure, the local erosion loss is consistent with the global erosion loss. The local erosion loss is approximately linear to the working time. The erosion loss of the first row occupies more than half of the global erosion loss. Ash impacts the first rows, and the kinetic energy decreases. Then most of the ash particles follow the gas flow, which hinders the particle collision in the other rows.

Similarly, as the ash size increases, the growth amplification of local erosion loss decreases in the first row. Meanwhile, with the ash size increased, the mass flow of ash impacting the economizer rises in the first row, which raises the probability of particle rebound and collision in other rows (compared with Figure 10a–h). For an ash particle diameter of 200 µm, these rebound phenomena occur in all rows. Second place is 150 µm, and 100 µm is the least. The probability of the particle rebound ultimately influences the growth amplification of global erosion. As the ash size increases to a certain numerical value, the probability of the particle rebound reaches the limit, which leads the growth amplification to decrease.

4.2. Evolution Process of Erosion

4.2.1. Maximum Erosion Depth

Figure 11 shows the maximum erosion depths of the different rows under various ash sizes. The maximum erosion depth is defined as the most critical region, as it represents the thinnest location of the tube wall. From this figure, the maximum erosion depth depends linearly on the working time at the initial time. However, compared to the linear trend (gray dashed line in Figure 11a,b), the growth of the maximum erosion depth slows down gradually. This growing trend forms a favorable flow passage, which prolongs the service life of the economizer effectively.

Since maximum erosion depth in the first row is maximal, this row is the most dangerous and of the greatest concern. Maximum erosion depth of the first row is double or even several times that of other rows. Moreover, not in all rows, the 200 µm ash particle impacts the most powerfully, and the maximum erosion depth is maximal. In the third, fourth, fifth, and eighth row, the maximum is under an ash particle diameter of 150 µm. In the sixth row, the maximum appears under an ash particle diameter of 100 µm. This is mainly caused by particle rebound in the first row, which drives other rows to remove more material as the rebounded particles impact them with different kinetic energies and incidence angles.

4.2.2. Erosion Profile

Considering that erosion in the first row is the most serious, this paper focuses on the evolution process of the erosion profile for the first row. Figure 12 shows the evolution process of the erosion profile with different ash sizes in the first row. From this figure, as time proceeds, a “V-shape” is gradually formed, which is similar to the erosion profile in Figure 1. This confirms that the dynamic mesh technology is advisable to describe the erosion of economizer tubes.

In addition, θ is introduced as the angle between the position of the erosion profile and the reverse direction of incoming flow, and θ_max represents θ in the location of maximum erosion depth in Figure 12a. From Figure 12, larger ash impacts the first row with a smaller angle of θ_max. As mentioned earlier, larger ash has greater inertia and less motion, ensuring that the location of the maximum erosion rate appears with a smaller θ.

4.3. Evolution Process of Particle Motion

Figure 13 shows the particle motion from the specific position in the inlet. From this figure, the particle collision point obviously moves to the center of the economizer tube over time. Meanwhile, due to more material loss with a larger particle, the larger particle moves across a longer distance. In order to describe the particle motion along the radial and tangential direction, Figure 14 is adopted to illustrate the collision point trajectory from the specific position in the inlet, which is the top view of Figure 13a. From this figure, with time increasing, the collision point trajectory moves along the increasing direction of the absolute value of θ. Just because of this, the probabilities of impingement and the growth of the maximum erosion depth slow down gradually in later stages. Although the trend of erosion loss is unlike the maximum erosion depth after 2 years, the erosion loss based on this result is predicted to slowly increase in the future.

4.4. Expiry Periods

In Figure 11a, the safe region is defined as the value below 1.99 mm (equal to the critical value of 4.01 mm in wall thickness [25], which is the minimum wall thickness by the specific theoretical calculation). When the maximum erosion depth is up to 1.99 mm, the first row of the economizer tube easily bursts. From this figure, the expiry periods are respectively 710, 530, and 440 days with 100 µm, 150 µm, and 200 µm. These data indicate the expiry period is shortened as the ash diameter increases. Moreover, as the ash size increases, the declining amplification of expiry periods is retarded.

5. Conclusions

The present research aims to investigate a novel method for the solution of the evolution process and its applications to the expiry period of an economizer bank. The CFD-DPM approach coupled with the UDF of the erosion model and the dynamic mesh technology is utilized to solve this complex issue. Based on the results of the applications, the following conclusions are made:

(1)

The CFD-DPM approach coupled with the UDF of Huang’s erosion model and the dynamic mesh technology can describe the evolution process of erosion on an economizer bank, especially the erosion profile; by comparing the simulation results with the erosion profile on-site, the correctness is verified for this proposed CFD-DPM approach.

(2)

The global/local erosion loss and the maximum erosion depth are linearly related to the working time, but the growth of the maximum erosion depth slows down gradually in the later stage; as the ash size increases, the growth amplification of global erosion loss, local erosion in the first row, and the maximum erosion depth decrease.

(3)

With increasing time, the collision point trajectory moves along the increasing direction of the absolute value of θ in the first row, which explains why the growth of the maximum erosion depth slows down in later stages.

(4)

The expiry period is shortened as the ash diameter increases; moreover, as the ash size increases, the declining amplification of expiry periods is retarded.