Numerical Analysis of Influence of Roughness of Material Surface on High-Speed Liquid Droplet Impingement

Fluid machinery is frequently significantly damaged by high-speed liquid droplet impingement on material surfaces. The liquid droplet impingement erosion occurs at the area where droplets impinge on fluid machinery at high speed, such as the last stage of blades in steam turbines and the elbows in steam pipes. In particular, the pinholes in piping generated by a pipe wall thinning in aging nuclear power plants can result in serious accidents, such as leakage of radioactive materials [1]; hence, it is critical to understand the influences of high-speed droplet impingement on pipe wall thinning rates. Therefore, many studies have been conducted on the liquid droplet impingement erosion. Heymann [2] showed quantitatively the influence of the droplet impingement velocity and droplet diameter on the liquid droplet impingement erosion rate. Sanchez-Caldera [3] derived the predict equation of erosion volume per droplet impingement. We also reproduced the power index of droplet impingement velocity for erosion rates earned from existing experimental studies by our in-house numerical analysis [4].

The erosion process advances from the incubation stage, to the acceleration stage, maximum rate stage, deceleration stage, and to the terminal or final steady-state stage over time [5]. These erosion processes are thought of the influence of the liquid film and roughness on the material surfaces.

In an actual fluid machinery, as many droplets exist and impinge on material surfaces constantly, it is assumed that the humidity and liquid film may exist there [6,7]. Regarding the influence of a liquid film on droplet impingement, it is thought that the impingement pressure is damped since the erosion of materials is suppressed. In existing studies, Shi et al. [8] reported that the droplet impingement is damped by the existence of a liquid film experimentally. Fujisawa et al. [9] conducted an experimental study on the effect of a liquid film on material under droplet impingement and proposed an equation of the erosion rate considering the effect of the liquid film. We have also numerically demonstrated that the pressure caused by a droplet impingement is damped by the liquid film [10].

On the other hand, the actual wall surfaces are not completely smooth, as the wall surfaces are eroded by repeated droplet impingement, and these material surfaces also have initial roughness from processing during manufacture. Therefore, it is considered that droplets impinge on the material surface on which roughness exists. It has also been reported that an elbow surface with a pinhole has numerous pits and that the surface is significantly rough as shown in Fig. 1. Regarding the influence of roughness on a material surface, it was observed that the erosion increases as the material surface roughed, and it was suggested that the impingement pressures close to a pit could be greater than those on a plane surface [11] and the erosion rate also increased [12]. However, few analyses explain the mechanism of which erosion progresses with keeping a pit shape even in liquid droplets impinging uniformly on the entire material surface.

In this study, by using an in-house fluid/material two-way-coupled numerical method that considers reflection and transmission on the fluid/material interface, the numerical analysis of the phenomenon of liquid droplet impingement on a rough surface with a single wedge and a single bump is conducted. The influence of roughness of the material surface in the liquid droplet impingement and the mechanism of erosion progress with keeping a pit shape is discussed.

In this study, the locally homogeneous model of compressible gas–liquid two-phase medium with phase changes [13–15] is used for the interface capturing method in a droplet interface. The governing equations in the fluid are the continuity and momentum equations, the total energy conservation equation of the gas–liquid mixture phase, and the continuity equation of the gas phase. In this numerical analysis, the continuum surface force model [16] is used in order to consider the surface tension. In addition, the phase change model of which the phase equilibrium theory of a flat gas–liquid interface [17] is extended to a homogeneous gas–liquid two-phase medium [18] is used. The saturated vapor pressure of water is calculated by the empirical formula given by Sugawara [19].

The governing equation system is closed by satisfying the equation of state of homogeneous medium and the total energy equation of two-phase medium.

In this study, to analyze the stress wave propagation in a material by droplet impingement, a material is treated as an elastic body, and the displacement and density change of the material is not considered. Here, the considered material deformation and density change are important to investigate the erosion process that develops by the impingement of many droplets. However, in this study, such long-term time change and transient phenomenon are not paid attention but the effect of a very short time like a high-speed single droplet impingement is analyzed. As it is assumed that the displacement and density change of a material are very small in the microscale single droplet impingement, the entire material is calculated with a fixed mesh system.

For the simulation of stress wave propagation in a material, the governing equations composed of the equation of motion and the time-differential constitutive equations of homogeneous isotropic elastic medium are simultaneously solved [20].

The magnitude of the stress in the material is determined by the equivalent stress, which is based on the concept of shear strain energy theory. In this study, the von Mises's equivalent stress is used as the material evaluation.

The cell-centered finite volume formulation is used to discretize the governing equations of both the fluid and material. The convective term is estimated AUSM type scheme with interpolation by using the third-order MUSCL-TVD method with a minmod limiter. The fourth-order Runge–Kutta method is used for the time integration.

The algorithm of two-way coupled numerical method in fluid/material interface [10] is applied as follows: the pressure and normal stress in the vertical direction on the material surface are obtained by considering the reflection and transmission of pressure and stress waves with acoustic impedance. Nonslip condition is adopted on fluid/material surface, where the vertical velocity on the boundary surface has values obtained by considering reflection and transmission. As explained in the previous Sec. 2.2, in this numerical analysis, as the displacement of the material is insignificant with a microscale single droplet impingement, the entire material, that includes its surface, is calculated with a fixed mesh system, although the displacement velocity is taken into account on the material surface. On the other hand, the tangential velocity on the boundary interface of the fluid has the same velocity of material which is calculated by the material analysis.

In order to investigate the influence of the roughness of material surfaces on droplet impingement and the influence of the difference of impingement points against the roughness, a three-dimensional (3D) simulation of vertical droplet impingement on wedge and bump surfaces is performed. The initial conditions, the calculation area, the shape, and the size of the wedge and bump are shown in Figs. 2 and 3, respectively. The impingement velocity is V = 200 m/s and the droplet diameter is d = 100 μm. The initial temperature and pressure of the water droplet in vapor are T = 293.15 K and p = 0.1 MPa, respectively. In addition, this droplet includes the mass fraction of the gas phase Y = 1.0 × 10⁻⁹. The material is assumed to be steel and the initial stress of the material is a compressive stress of 0.1 MPa. The material properties for steel are density ρ_s = 7800 kg/m³, Young's modulus E = 200 GPa, and Poisson's ratio ν = 0.3. In the boundary of the calculation area, the temperature is fixed at 293.15 K.

Under these calculation conditions, because of the assumption of a symmetrical phenomenon for one aspect, only half of all the space in the depth direction is calculated. The calculation grid is a 3D orthogonal grid and the calculation area is 2.0d × 1.5d × 1.0d ((x-direction) × (y-direction) × (z-direction)) on the fluid side, and 2.0d × 1.0d × 1.0d on the material side. The buffer area is placed around the primary calculation area. The grid numbers are 200 × 200 × 100 on the fluid side and 200 × 100 × 100 on the material side. The grid spacing in the horizontal direction on the material surface has equidistant grid spacing of 0.01d on both the fluid and material sides. For the grid spacing in the vertical direction on the material surface, the calculation area on the fluid side has a minimum grid spacing of 0.0025d and a maximum grid spacing 0.02d. The grid spacing is gradually decreased as it approaches the material boundary, and the material side has equidistant grid spacing of 0.01d. In addition, the grid spacing of the liquid film has an equidistant grid spacing of 0.0025d per mesh in the vertical direction on the material surface.

The positional relationship between the droplet and the wedge or bump is shown in Figs. 2 and 3. The width 2w of both wedge and bump is 20 μm (0.2d), and the depth of the wedge or height of the bump is 2 μm (0.02d). In Figs. 2 and 3, l is the offset from the wedge or bump center to the impingement center of the droplet. Hereafter, the dimensionless distance l/w, normalized l by w, is used for the offset distance.

Figures 4–6 show the comparison of time variations of the pressure in the fluid and the equivalent stress in the material at each offset l/w = 0 (Fig. 4), 1 (Fig. 5), and 2 (Fig. 6) of (a) wedge and (b) bump. As a reference, the result for a flat surface is also shown in Fig. 4(c). The black line is assumed as the gas–liquid interface at a void fraction of α = 0.5. In the calculation results, time t = 0 indicates the moment when a droplet contacts somewhere on the material surface.

In the case of an offset of l/w = 0 for the wedge, as the droplet impinges on the wedge surface, two contact edges occur and move to the inside and outside (Fig. 4(a), (i)). The inside contact edge concentrates at the bottom of the wedge and a significant pressure results there (Figs. 4(a), (ii) and (iii)). This result qualitatively agrees with the study of Field et al. [21]. This high pressure also causes the equivalent stress to increase considerably.

On the other hand, in the case of an offset l/w = 0 for the bump, the droplet impinges on the top of the bump, and the contact edge moves on the downslope of the material surface. While the contact edge is moving on the downslope (Figs. 4(b), (i) and (ii)), the pressure increase is smaller than that on a flat surface (Fig. 4(c), (ii)). This is because the speed of the continuous impingement between the droplet interface and the material surface is slower compared to the droplet impingement on a flat surface because of the downslope, and the speed of the initiation of pressure waves also decreases. Therefore, as the concentrations of pressure waves become suppressed, the pressure does not increase as much as that of the impingement on a flat surface. After the contact edge reaches the flat surface, a pressure concentration is observed, however, this pressure is significantly smaller than that of an impingement on a wedge surface.

Even when comparing immediately after the droplet impingement, it is shown that the impingement pressure of the bump is less than those of other two surfaces. In addition, in the initial stage of the droplet impingement, it can be also seen that the behaviors of droplet impingement differ from depending on the surface condition.

In the case of an offset of l/w = 1 on a wedge, while the contact edge of the wedge direction is moving on the downslope, the pressure increase is smaller than that of the contact edge moving on a flat surface (Fig. 5(a), (i)). For the same reason as discussed above, the pressure does not increase as much as that of the impingement on a flat surface because of the downslope. However, when this contact edge reaches the bottom of the wedge, a pressure concentration occurs because of changing from a downslope to an upslope (Fig. 5(a), (i)). This is the opposite of the case of the downslope, as the speed of the continuous impingement between the droplet interface and the material surface is increased by the upslope, the pressure concentration is not difficult to achieve. In addition, the occurrence of side jets is also suppressed by the upslope; therefore, the pressure concentration occurs during the contact edge moving on the upslope.

In the case of an offset of l/w = 1 on a bump, the droplet impinges on the slope of the bump and the contact edges move on each upslope and downslope (Fig. 5(b), (ii)). When the contact edge on the right side reaches the corner, the pressure concentration occurs with the significant change in the slope (Fig. 5(b), (iii)). In both cases, for the wedge and the bump, the pressure concentration initiates from the depressed point where the gradient increases in the moving direction of the contact edge.

Regarding the results of an offset of l/w = 2 on a wedge, in the initial stages of impingement, the behavior of the impingement is similar to the case of the flat surface (Figs. 6(b), (i) and (ii)). However, in this case, the pressure increase at the bottom of the wedge is insignificant as the pressure waves decay before reaching the bottom (Figs. 6(a), (iv) and (v))). Therefore, as the offset distance increases from a wedge center, the behavior of the impingement approaches that of a flat surface because of the effect of the wedge.

In the case of an offset of l/w = 2 of a bump, in the initial stages of the impingement, the behavior is similar to that of a flat surface (Figs. 6(b), (i) and (ii)). However, as the impingement advances, the contact edge on the side of the bump reaches the upslope of the bump and then a pressure concentration occurs as discussed in the previous paragraph (Fig. 6(b), (iii)). This pressure is large as compared to the maximum pressure of impingement on a flat surface.

With existence of an offset, as the contact edge due to the droplet impingement reaches the recessed portion of wedge or bump, the high pressure and equivalent stress occur at the time corresponding to the offset distance.

In the case of droplet impingement on the wedge surface, for the offset l/w = 0, as the pressure concentration occurs at the wedge center, the p_max and (σ_eq)_max are observed at the wedge center and the increasing ratios of p_max and (σ_eq)_max are approximately 1.7 times and 1.6 times greater than those of the flat surface, respectively. For the offset l/w = 1, as the pressure concentration is generated by the upslope of the wedge, the p_max and (σ_eq)_max occur from the bottom to the upslope of the wedge. At this time, the increasing ratios of p_max and (σ_eq)_max are 1.1 and 1.4, respectively. These values are smaller than those of offset l/w = 0 but greater than those of the flat surface. For the offset l/w = 2, as the offset has little effect of increasing p_max and (σ_eq)_max from that of a flat surface, the increasing ratio is approximately 1 and the values of p_max and (σ_eq)_max are similar to those of a flat surface. Therefore, as the offset increases, the increasing ratio γ decreases and the impingement behavior of the wedge surface approaches that of a flat surface. In the case of droplet impingement on the bump surface, for the offset l/w = 0, as the pressure concentration is inhibited by the downslope, p_max and (σ_eq)_max are marginally smaller than those of a flat surface. For the offsets l/w = 0, 1, and 2, as the pressure concentration occurs at the corner of the depression point, the increasing ratios are greater than those of a flat surface. From these results of increasing p_max and (σ_eq)_max for both wedges and bumps, the roughness of material surfaces has an effect of increasing p_max and (σ_eq)_max on the droplet impingement.

Typically, for a material like steel, it is known that the erosion rate increases instantaneously at the transition from the incubation stage to the acceleration stage. In these calculations, the fact that p_max and (σ_eq)_max of the droplet impingement on the rough surface become greater than those of a flat surface can explain the transition from the incubation stage to the acceleration stage in actual erosion. The mechanism is as follows: first the concave and convex phases occur without mass loss, and once the concave and convex phases occur, even though the droplets impinge evenly on the material surface, the equivalent stress generated by the droplet impingement becomes greater, and as a result the mass loss could occur.

In all cases, except for the bump with an offset of l/w = 0, the greater pressure compared to that of a flat surface occurs in the depressed area where the gradient increases in the moving direction of the contact edge, and the equivalent stress also increases there. In particular, as shown in the case of the wedge with an offset of l/w = 0, the significantly greater pressure occurs in the depressed area where a number of contact edges concentrate simultaneously. On the other hand, the pressure concentration does not occur at positions like the top of the bump. Therefore, there is a possibility that the erosion speed increases in the depressed area, and on the contrary, is smaller at the protruding point. These results suggest that the impingement point of the droplet does not erode, but a certain spot, depending on the shape of the material surface, erodes. Figure 8 shows the summation of the individual distribution of (σ_eq)_max of five droplet impingements on wedge and bump surfaces (l/w = −2, −1, 0, 1, and 2). Although for actual droplet impingements it is considered that a number of droplets impinge simultaneously or over time, and that the residual stress should also be present in the material, these influences are not considered in Fig. 8. From Fig. 8, the possibility exists that, in the case of a wedge, the center part of the wedge erodes deeply, and in the case of a bump, the bump base will erode deeper, therefore, the concave and convex shapes become greater. It is considered that this result indicates that erosion is formed in the shape of pits even though droplets impinge evenly with the material surface.

In this study, in order to investigate the influence of roughness of the material surface on droplet impingement, by using an in-house fluid/material two-way coupled numerical method that considers reflection and transmission on the fluid/material interface, the numerical analysis of the phenomenon of liquid droplet impingement on a rough surface with a single wedge and a single bump was conducted and the distributions of elastic stresses were compared.

From the numerical results of droplet impingement with different offsets from the center of the wedge and bump, the maximum equivalent stress becomes greater both on wedge and bump surfaces than on a flat surface, no matter where the droplet impinges, except for the center of the bump. The results explain the transition of the erosion from the incubation stage without mass loss to the acceleration stage with mass loss, where once the roughness appears, the erosion suddenly accelerates.

The possibility exists that the center part of the wedge is deeply eroded but the bump head remains free from erosion. The result explains the mechanism why erosion pitting occurs with droplet impingement, even though the droplets impinge homogeneously across the material surface.

This work was supported by AMED under grant Number JP17he1302008 and J. MORITA MFG. CORP. The numerical simulations were performed on the Supercomputer system “AFI-NITY” at the Advanced Fluid Information Research Center, Institute of Fluid Science, Tohoku University.

d =

droplet diameter (μm)

D =

depth of wedge or bump (μm)

E =

Young's modulus (GPa)

l =

offset from wedge or bump center to droplet impingement point (μm)

l/w =

dimensionless offset

p =

pressure (MPa)

T =

temperature (K)

Y =

mass fraction of gas

w =

width of wedge or bump (μm)

α =

void fraction

γ_p =

increasing ratio of maximum pressure

γ_σ =

increasing ratio of maximum equivalent stress

σ_eq =

equivalent stress (MPa)