An LES study of secondary motion and wall shear stresses in a pipe bend

Supercritical carbon dioxide (⁠$sCO2$⁠) is an ideal working fluid for power generation in concentrated solar, fossil fuel, and nuclear power plants as it is a nontoxic, nonflammable, and low-cost fluid. The high-pressure (200–350 bar), high-temperature (550–$750 °C$⁠), and high-density fluid enables extremely compact and high efficiency turbomachinery designs. However, severe erosion has been observed under several occasions in the sCO₂ power cycle test loops, particularly in the turbine nozzles and blades. Since the pipelines usually bear erosion and corrosion simultaneously in multiphase conditions, it was conjectured that erosion-corrosion (E-C) type of behavior, i.e., small oxide particles and impurities in sCO₂, caused the erosion and failure of the blades. Recently, Finn and Doğan¹ performed large-eddy simulations of turbulent flows over a turbine nozzle in real operating conditions and modeled solid particle collision on the turbine blades using one-way coupling. Based on the simulation results and the Finnie model,² the erosion magnitudes on the turbine blade surface were estimated. It was found that the trailing edge erosion is roughly two orders of magnitude larger than that on the leading edge, consistent with the observations reported by Fleming et al.,³ suggesting that mechanical erosion, caused possibly by the small solid oxide particles spalled from the internal walls of the pipe, could be a source of severe erosion in the power generation systems.

However, erosion of the material was also observed even with a fairly pure sCO₂ flowing through the power cycles.⁴ It is, thus, critical to identify the potential sources of erosion in this pure flow scenario with minimal particulates. The wall thinning has long been an important safety issue for complex pipeline systems including elements such as elbows and orifices,⁵ and the so-called flow accelerated corrosion (FAC) is the major cause for this phenomenon.⁶ FAC usually happens at the oxide layer of the carbon steel wall where the mass transfer of iron ions into the turbulent bulk flow occurs through a diffusive process,^7,8 which could further promote the E-C type of erosion on a carbon steel wall. A wall thinning issue caused by the flow turbulence has been observed particularly in pipe bends and orifices.^9–11 Hence, it is reasonable to speculate that FAC first contributes to the generation of small iron particulates at complex geometries (bend or orifice) in the pure sCO₂ power cycles by thinning the carbon steel wall; then, the particles entrained in the flow are responsible for the erosion of the turbine parts. This speculation aligns well with the observation reported by Fleming and Kruizenga⁴ that high iron particles were detected in various areas on the flow path, and the numerical study of the mechanical erosion on turbine blades reported by Finn and Doğan.¹ The current work focuses on the first half of the conjuncture to explore the potential cause of the erosion in relatively pure sCO₂ systems.

It is noteworthy to mention the physical underlying of the mass transfer associated with FAC has not been fully understood due to complex geometry and the secondary flows generated by the curved pipe.⁵ Limited studies were published to investigate such phenomenon. Recently, Takano et al.¹² experimentally studied how the swirling flow could impact the mass and momentum transfer of a pipe with elbow and orifice. It was found that the nonaxisymmetric flow structure is the major mechanism responsible for the increased mass transfer downstream of the pipe. Utanohara and Murase¹³ conducted experimental work on the influence of velocity and temperature on FAC in the pipe bend by measuring the FAC rate at various velocity and temperature settings. The results indicate that the combined effects are case dependent and particularly strong at the elbow. Some researchers argued the FAC may be attributed to large fluctuations of wall shear stress.¹¹ Lister and Lang¹⁴ proposed that the shear stress is directly leading to the mechanical shearing of the outer oxide layer.

In the current work, it is hypothesized that the substantial shear stress on the pipe walls and its spatiotemporal variation due to turbulence and secondary flow patterns are the major causes of erosion in the pure sCO₂ power systems. To test the hypothesis and to gain better understanding of the stress distributions on the inner walls of pipe bends, predictive numerical simulation of a single-phase, turbulent flows without heat transfer in pipe bends at representative conditions and flow Reynolds numbers are performed. Considerable work exploring the dynamics of laminar and turbulent flow through pipe bends has been conducted.^15–17 Flow in curved pipes with bends involves flow turning that results in centrifugal forces on the fluid elements and corresponding pressure field to balance these forces. Formation of counter-rotating vortices, commonly termed as Dean vortices (Fig. 1), are observed as the fluid elements with higher velocity are forced to the outside of the bend and lower velocity elements are forced toward the inside.^15,16 These curvature-induced swirling motions and their strength is described in terms of the Dean number, De, where $De=D/2RcReD, ReD=DUb/ν$ is the Reynolds number, D is the diameter of the pipe, R_c is the radius of curvature of the bend measured to the pipe centerline, U_b is the bulk velocity, and ν is the kinematic viscosity of the fluid. Under unsteady flow conditions, the Dean vortices also exhibit transient behavior. The two counter-rotating swirling motions start to be dominant alternately over space and time. This unsteady behavior of Dean vortices is called swirl-switching.¹⁸ The transient behavior of swirl-switching, time scales related to oscillations of these vortices, and resultant variations in the shear stresses on the walls in turbulent flow have not been thoroughly investigated, especially at large Reynolds numbers.

Rütten et al.¹⁶ investigated turbulent flows in a $90°$ pipe bend using large-eddy simulations at Re_D of 5300–27 000. Their work indicated that the overall forces on the pipe walls showed a distinct peak in the power spectra related to the vortex shedding at the inner side of the bend. At large Re_D, the power spectra also exhibited a frequency of oscillation much lower than the vortex shedding from flow separation. This was attributed to two Dean vortices whose strength varies in time dominating the flow field. This also suggests that local variation in the wall shear forces, especially after the bend, can be important in magnitude as well as time scales. Such variations in local shear stresses may be relevant to the erosion problem at very large Re_D.

Furthermore, Hellström et al.¹⁷ used time-resolved stereoscopic particle image velocimetry (PIV) method to study turbulent flows at multiple Re_D from 25 000 to $ReD=115 000$ through a $90°$ bend. Snapshot proper orthogonal decomposition (POD) method was applied to identify the Dean motion and swirl-switching, as well as extracting the corresponding POD modes. The time coefficient of the first POD mode was found to be oscillating at the frequency, which is close to that of the overall force reported by Rütten et al.¹⁶ Although Re_D of the two cases are not exactly the same, this leads to an open question as to whether the swirl-switching can cause the oscillation of the shear forces on the pipe downstream of the bend and be responsible for the material erosion observed in the power cycles.

More recently, Carlsson et al.¹⁹ performed large-eddy simulations of turbulent flows at $ReD=3.4×104$ through pipe bends with four different curvature ratios $γ=0.16,0.25,0.35, and 0.5$⁠, which is defined as $γ=Rc/D$⁠. Two distinct types of swirl-switching were observed at high (Strouhal number $St=0.5∼0.6$⁠) and low (⁠$St≤0.01$⁠) frequencies, respectively. The Strouhal number is a nondimensional frequency and defined as $St=fD/Ub$⁠. The authors argued that the high frequency switching roots from the bend pipe; while the low frequency signal results from the upstream turbulent pipe flows. However, Hufnagel et al.²⁰ disagrees with this conclusion, where direct numerical simulations of flow at Re_D = 11 700 in a bend pipe with $γ=0.3$ were conducted. Three dimensional POD analyses were applied and the traveling waves were identified to be the cause of the swirl-switching rather than the straight pipe section upstream of the pipe bend. One important distinction of the simulation settings between Carlsson et al.¹⁹ and Hufnagel et al.²⁰ is that the recycling inflow boundary was used in the former study, while a synthetic eddy method was employed by the latter. Hufnagel et al.²⁰ claims that the artificial frequency introduced by the recycling inflow boundary attributes to the connection between swirl-switching and upstream pipe flow established by Carlsson et al.¹⁹ Besides, Ikarashi et al.²¹ measured the velocity fields of very high Re_D flows (⁠$1×105∼3×105$⁠) through a pipe with an short elbow (⁠$γ=1.0$⁠) using PIV. They inferred that the periodic velocity fluctuations are generated by the oscillation of flow structures downstream of the pipe bend. However, no shear stress/force data were reported in any of these recent publications.

In the present work, we simulate the turbulent flows at three different Re_D of 5300, 27 000, and 45 000 through a $90°$ bend to first validate our results against those available data from experiments as well as other simulations and then investigate the nature of shear force variations at different locations after the pipe bend. In addition, the snapshot POD approach is applied to detect the dominant POD modes and any frequency associated with the swirl-switching. The power spectral density functions are computed for both shear force and POD modes to reveal the connections between them. The rest of the paper is arranged as follows. Section II describes the mathematical formulations and numerical scheme used to solve the governing equations. Section III briefly goes through the simulation setups including geometry and inlet boundary condition as well as grid resolutions. The results are presented and discussed in Sec. IV. Section V states the summary and conclusions of the study.

The governing equations, the dynamic Smagorinsky model, and the numerical solver used in this work are briefly described below.

In the present work, the standard mathematical formulations for the single-phase, incomprehensible large-eddy simulations (LES) are employed as

where

and $ũi$ is the filtered velocity, $μ¯$ is the viscosity, and $p¯$ is the pressure. The subgrid scale unclosed transport terms in the momentum equation are grouped into the residual stress $qijr$⁠. The dynamic Smagorinsky model by Mahesh et al.,²² Moin and Apte²³ is used. The dynamic model has no adjustable constants, and thus allows evaluation of the predictive capability of the model. The subgrid scale stress, $qijr$ is modeled as

where the isotropic part of the subgrid stress q² is absorbed into the pressure term, and the subgrid viscosity ν_t is given as

Here, $Δ¯$ is the filter width and the coefficient $Cμ$ is obtained using least squares approach following the dynamic procedure.²⁴ 

An energy-conserving, finite-volume scheme for unstructured, arbitrary shaped grid elements is used to solve the fluid-flow equations using a fractional step algorithm.^22,23,25 The velocity and pressure are stored at the centroids of the volumes. The cell-centered velocities are advanced in a predictor step such that the kinetic energy is conserved. The predicted velocities are interpolated to the faces and then projected. Projection yields the pressure potential at the cell-centers, and its gradient is used to correct the cell and face-normal velocities. A novel discretization scheme for the pressure gradient was developed by Mahesh et al.²² to provide robustness without explicit numerical dissipation on grids with rapidly varying elements. This algorithm was found to be imperative to perform LES at high Reynolds number in complex flows. The overall algorithm is second-order accurate in space and time for uniform orthogonal grids. A numerical solver based on this approach was developed and shown to give very good results for both simple²⁶ and complex geometries²³ and is used in the present study.^27,28 A thorough verification and validation of the algorithm was conducted^29,30 to assess the accuracy of the numerical scheme for test cases involving two-dimensional decaying Taylor vortices, flow through a turbulent channel and duct flows,²⁹ and particle-laden turbulent flows.^23,26

The computational geometry and the boundary conditions are provided in detail in this section. Three different Re_D flows were simulated at 5300, 27 000, and 45 000. Among them, the lowest Re_D flow (5300) mostly serves as a validation case, whose results were compared with the published numerical and experimental data set. The results from the rest of the cases with higher Re_D were analyzed thoroughly to investigate the characteristics of such flows. A simulation in a periodic pipe was performed at each Re_D to generate proper inlet flow data for the pipe bend. The grid resolutions are also listed and discussed.

The computational setup used consists of a circular cross section pipe with a $90°$ bend as shown in Fig. 2. The pipe consists of a straight section of length 8D, followed by the bend section and another straight section of 8D length. The bend with the curvature ratio of 1 is investigated. Owing to the bend geometry an unstructured grid with predominantly hexagonal mesh is used for calculations. The so-called O-grid mesh is generated using a commercial software ANSYS ICEM-CFD.

When conducting large-eddy simulations, it is important to impose consistent, continuity-satisfying fluctuations in the velocity field at the inlet³¹ to create proper turbulence structure and velocity fluctuations in downstream. The use of random fluctuations scaled with turbulence intensity values has been proposed. However, this approach can lead to convergence issues in a conservative numerical solver. As a result, a separate periodic pipe flow (length 5D as used in Rütten et al.,¹⁶ see Fig. 2) is simulated at the desired mass-flow rate and Re_D using a body-force technique.³¹ The turbulent flow field in this cycling pipe was verified by comparing the first-and second-order turbulent statistics with those from straight pipe flow at the same Re_D. Then, the instantaneous velocity field at the outlet face of the periodic pipe is stored at every time step and is used as the inlet flow condition for the pipe bend. This provides a fully developed, instantaneous, turbulent velocity profile at the inlet of the pipe bend. A convective outlet condition is enforced at the outlet. To avoid the possible artificial frequency associated with this process,²⁰ the time span of the inflow data were chosen to be long enough (>60 flow-through times for the flow in the pipe bend). In this study, only an isothermal flow field is simulated and thus effects of temperature and density variations and properties of sCO₂ are not included. According to Utanohara and Murase,¹³ both velocity and temperature have strong influence on the erosion rate, and the combined effects are rather complicated than a simple linear combination from the two aspects. Later study involves turbulent flow with heat transfer where these effects do become important.

The flow fields simulated are for three different Re_D, 5300 and 27 000 similar to the work by Rütten et al.¹⁶ and 45 000 which is one of the experimental cases conducted by Hellström et al.¹⁷ The grid resolutions for the three cases are listed in Table I. Stretched grids are generated based on the hyperbolic tangent function. The grids in the wall normal (radial) direction have a minimum grid resolution near the wall in wall units of $Δr+∼0.8$ for Re_D = 5300, less than 2 for Re_D = 27 000 and 45 000. Nearly uniform grids are used in the axial direction with a resolution of $Δz+∼20$ for both Re_D = 5300 and 27 000, while above 30 for Re_D = 45 000. The maximum resolution in the circumferential direction is on the order of $(rΔθ)+∼6$ for Re_D = 5300, around 20 for Re_D = 27 000 and 31 for Re_D = 45 000. Grid resolutions in the azimuthal direction used in this work for $ReD=5000 and 27 000$ are finer than those used for LES by Rütten et al.¹⁶ (14.1 for Re_D = 5300 and 41.9 for Re_D = 27 000); Re_D = 45 000 is comparable to the LES case (Re_D = 37 600) by Chin et al.³² In the present wall-resolved large-eddy simulation, no-slip conditions are applied at walls.

In this section, the simulated results are presented and discussed. The mean and fluctuating velocity profiles at multiple locations of both upstream and downstream of the pipe bend are shown and compared with available reference data. Then, snapshot POD method is applied to analyze the flows at $ReD=27 000 and 45 000$⁠. Three most energetic modes at three xy-planes (⁠$z/D=1, 3, and 5$⁠) are visualized. The power spectral densities of the associated time coefficients are computed. Besides, the shear forces on the pipe wall were analyzed in the frequency domain to identify any distinguished frequencies, which are found to be identical to those for the swirl-switching phenomenon.

To validate the large-eddy simulation results in the present study, the inflow data obtained in the straight pipe section are analyzed first. The RMS velocity components at $ReD=5300 and 27 000$ are compared with data from multiple references.^16,33,34 Both the mean and fluctuation velocity profiles on multiple cross sections of the pipe bend are shown to demonstrate the effect from pipe bend to the upstream flow. Also, the velocity profiles at one diameter (1D) downstream of the pipe bend for Re_D = 5300 are compared with data from both experimental³⁵ and numerical¹⁶ studies. Finally, the pressure distribution on pipe wall and pipe center is reported and set side by side with the reference data.

Figure 3 provides the RMS velocity profiles in the straight pipe, which was employed to generate inflow data, together with data from multiple reference, where Re_D = 5300 is illustrated in Fig. 3(a) and Re_D = 27 000 in (b). The LES results by Rütten et al.¹⁶ and DNS results by Unger³³ are given as well for comparison in Fig. 3(a). The data from two LES studies^16,34 are illustrated in Fig. 3(b). All three components of the RMS velocity are matching well with the reference data, suggesting the good credibility of the current simulation settings. In addition, the mean velocity and RMS velocity profiles of Re_D = 5300 at various upstream locations, e.g., −0D, −1D, and −2D as illustrated in Fig. 2, are plotted in Fig. 4. Here, “inner side” means the profile is plotted along the inner side line of the pipe, and likewise the “outer side” and “center” [see Fig. 4(g)]. It is observed that the mean velocity profiles before 1D of the bend all collapse together, and coincide with that of a fully developed straight pipe. This is consistent with what has been observed in Vester et al.³⁶ that the flow is not influenced by the bend before one diameter of the pipe. The same characteristic can be obtained from RMS velocity profiles as well. The velocity profiles in Fig. 4 are in good agreement with those in Rütten et al.,¹⁶ where the same flow was simulated, though they are not included here.

The velocity profiles along the centerline at 1D downstream of the bend are shown in Fig. 5. The LES data are compared with PIV data from Brücker³⁵ and the LES results by Rütten et al.¹⁶ It is observed that the mean velocity profile matches well with both experimental and numerical data. The maximum relative difference was estimated to be 7%. The profile of the axial component of RMS velocity in Fig. 5(b) also shows a good comparison with experimental and published numerical data. The RMS velocity near the pipe walls agrees well with the LES study by Rütten et al.¹⁶ 

Besides, the pressure distribution along the pipe wall and pipe center at Re_D = 5300 is investigated and plotted in Fig. 6(a) together with the LES data from Rütten et al.³⁷ The abscissa represents the nondimensional distance along the wall starting at 3 diameter before the bend. The ordinate describes the nondimensional pressure drop computed by $Δpn=Δp/(ρ2U¯b2)$⁠. The results again match well with the reference data.

The flow field is visualized to provide an overall understanding of the flow features. Figures 7(a) and 7(b) show the contours of the mean velocity magnitude field nondimensionalized by the bulk velocity (⁠$|Um|=ux¯2+uy¯2+uz¯2Ub$⁠) on the symmetry plane at $ReD=27 000 and 45 000$⁠, respectively. No significant difference of the mean velocity magnitude was observed between the two Re_D flows. The high velocity region from the straight pipe first gets closer to the inner wall starting about 1D before the bend. Then after about $30°$ of curvature angle, it deviates from the inner wall, moves toward the outer wall and keeps this way downstream. This flow feature is consistent with what was reported by Sudo et al.³⁸ 

Besides, the nondimensional in-plane mean velocity magnitude field (⁠$|Uxy|=ux¯2+uy¯2Ub$⁠) overlaid with streamlines is depicted at the cross sections of pipe at $z/D=1 and 3$ [see the dashed line in Fig. 7(a)]. Figures 7(c) and 7(d) are for $z/D=1$⁠; Figs. 7(e) and 7(f) $z/D=3$⁠, where the left-hand side is for Re_D = 27 000 and right-hand side Re_D = 45 000. Symmetric distribution of the in-plane mean velocity was observed as expected. Similar to the observation from mean velocity magnitude field, limited difference was spotted for the in-plane mean velocity distributions between the two Re_D flows. A stagnation point is present at the centers of both inner (right-hand side in the figure) and outer (left-hand side) sides of the pipe. Two counter-rotating vortices are generated due to the centrifugal force, with a secondary vortex existing solely at $z/D=1$⁠. The high velocity regions are mostly concentrated near the outer wall, where large magnitudes of shear forces are guaranteed. Since the velocity fields are time-averaged, the unsteady swirl-switching is not detected in these plots and a more detailed temporal analysis is needed.

In order to identify the swirl-switching in the turbulent flows through the pipe bends, POD is applied to extract the flow features. POD is capable of capturing the optimal modes from any field being examined bases on the energy content. In other words, the optimal basis vectors are computed that can best represent the given data. The snapshot POD method can greatly reduce the computational effort compared to the original POD approach. Hence, this method has been widely used in the fluid dynamics to detect the coherent structures and flow patterns.^17,39 In this work, the mathematical formulations for the snapshot POD method described in Taira et al.⁴⁰ are adopted and briefly discussed in the following. For any given vector field $q(x,t)$ with the temporal mean $q¯(x)$ removed, the remaining fluctuation part can be decomposed as

where $ϕj$ represents the POD modes, and a_j the corresponding time coefficients. The POD mode is defined as a function of location only; while the temporal coefficient depends on time. In this way, the unsteady components of the vector field $q(x,t)$ can be decomposed into temporal and spatial domains, respectively.

In this work, the velocity fields at three different locations: $z/D=1, 3, and 5$ after the bend are investigated using the snapshot POD method. The snapshots of the velocity fields at these locations are stored in a specific frequency that is much higher than the target frequency of the POD modes to avoid aliasing. Table II provides the detailed time interval and span for generating the data sets used to compute the POD modes. Figure 8 shows the relative and integrated energy levels for the first 20 modes at three cross sections for both high Re_D flows. It shows that the first three modes are the most energetic ones, which takes over 50% of the total energy and are illustrated in Figs. 9 and 10. The modes are distinguished by black solid lines and are overlaid by the normalized streamwise velocity field that is pseudo-colored.

The first mode is a single cell covering the whole pipe cross section regardless of mode number and Re_D. This mode, consistent with what is reported by Hellström et al.¹⁷ for flows at $ReD=25 000and 115 000$⁠, represents the swirl-switching phenomenon. Since the modes extracted by POD are orthogonal, the single cell mode observed can represent the two counter-rotating vortices which share the same basis, as pointed out by Hellström et al.¹⁷ 

The Dean-like vortex is recovered in the second and third modes, where two vortices can be identified, indicating that the Dean vortices have less energy and are weaker than the swirl-switching. In fact, the Dean-like motions could be seen as the transient period of the swirl-switching. It is noticed that the second mode at $z/D=1$ is not a typical Dean motion compared to that at $z/D=3 and 5$ for both $ReD=27 000 and 45 000$⁠. This may be attributed to the fact that $z/D=1$ is too close to the bend. The Dean motion is still under development at this specific location. Another interesting observation is, the second and third modes are nearly perpendicular to each other at both $z/D=3 and 5$⁠, demonstrating the bimodal characteristic of such type of flows. This feature was not present at $z/D=1$⁠, possibly due to the incompleteness of the Dean motion development. Hellström et al.¹⁷ did not report such abnormal behavior because all of the data reported were at or after 5D of the bend. In addition, no substantial difference in the POD modes' structures was observed between Re_D = 27 000 and Re_D = 45 000.

Moreover, the power spectral density functions for the time coefficients a_j are plotted in Figs. 11(a)–11(f), representing the temporal characteristics of the POD modes. The first row is for PSD at xy-plane of $z/D=1$⁠, second $z/D=3$⁠, and third $z/D=5$⁠. The figures in the left column are for Re_D = 27 000, right column Re_D = 45 000. Aiming to reduce the noise and focus on the meaningful frequencies, the PSD of the time coefficients are filtered using Bartlett window with the degree of freedom of 60. The 95% confidence interval is included in each figure as reference. The abscissa of the figure is the Strouhal number. A very well pronounced peak at St = 0.25 is observed for Re_D = 27 000, and it shifts to St = 0.28 at Re_D = 45 000. Those distinguished peaks suggest that the motion of swirl-switching is oscillating at such frequencies. The sole dependency of the peaking frequency on Re_D indicates that this oscillation happens over the whole pipe section after the bend, thus is not spatially dependent.

The PSD for mode 2 has a low frequency peak at $St∼0.01$ at all locations regardless of Re_D. Carlsson et al.¹⁹ claimed the root of this low frequency peak is the very-large-scale motions (VLSMs) created in the upstream pipe flow; while Hufnagel et al.²⁰ argued it results from the artificial signal created by the recycling inflow data. As mode 2 exhibits the low frequency peak disregarding the location and Re_D, the latter explanation seems plausible. However, with a closer look, the low frequency signal in mode 3 gradually becomes more intensified over the distance away from the bend. In addition, such tendency is more clearly observed at higher Re_D flow (Re_D = 45 000), where mode 3 has almost the same magnitude of low frequency peak as that for mode 2 at $z/D=5$⁠. Such observation suggests the peaking low frequency signal may be related to both location and Re_D at least for mode 3, which contradicts with the arguments from Hufnagel et al.²⁰ but aligns with the conclusion by Carlsson et al.¹⁹ Last but not least, due to the very long time span for the inflow data used in the present work (>60 flow-through time), the low frequency is very unlikely to come from the recycling inlet condition.

One significant difference of the PSD between the two Re_D cases is that modes 1 and 2 for Re_D = 45 000 share the same peak at all three xy-planes; while the overlapping of the peak frequency is only heavily pronounced at $z/D=5$ for Re_D = 27 000. Such different behavior has not been reported by other research, to the best of the authors' knowledge. One possible reason is that, the length of pipe for flow to achieve the fully developed Dean motions is depending on Re_D. This is yet to be confirmed by further studies with even higher Re_D flows.

The shear forces caused by the fluid motions are obtained everywhere on the pipe surface. The overall forces are analyzed in the frequency domain to identify any oscillation. The shear forces have been nondimensionalized by the stagnation pressure of the bulk velocity $12ρUb2$ and the squared pipe diameter D² following Rütten et al.¹⁶ It is observed in Figs. 12(a) and 12(b) that the PSD has the peak frequency at St = 0.25 for Re_D = 27 000 and St = 0.28 for Re_D = 45 000. These two frequencies are exactly matching with those of the POD modes. This is a strong evidence that the swirl-switching is highly correlated with the oscillation of the total shear force. This suggests that the turbulent flow induces oscillations of wall shear forces which can be a cause of the erosion, especially if the force magnitudes are large. For real sCO₂ power plant, the operating conditions can lead to extreme high Re_D (⁠$∼20$ million), and hence much higher shear forces. These oscillating shear forces can remove material from the internal surface of pipe such as oxide scale formed as a result of a high-temperature reaction between CO₂ and steel surface.⁴¹ In many instances, oxide layer growing on the steel surface is loosely adhered to the surface and prone to spallation under external forces such as oscillating shear forces demonstrated above or even due to internal growth stresses.⁴² Once the spalled oxide particles become entrained in the CO₂ flow, they can accelerate erosion through hard particle impact on the internal surfaces of downstream components.¹ In addition, changes in temperatures and thermal properties of the fluid can influence the forces and their spectra as well.¹³ 

Large-eddy simulations of turbulent flows through a $90°$ pipe bend with unit radius of curvatures (⁠$Rc/D=1$⁠) were performed using an unstructured grid solver. Flows at three different Re_D (5300, 27 000, and 45 000) were simulated to understand the behavior of the flow field within the pipe bends, the dynamics of the Dean vortices, and variations in the shear stress. First, the solver was validated by comparing the mean and RMS velocity profiles before and after the bend with available experimental and numerical data. Then, the snapshot POD method was applied to detect the optimal modes and identify the swirl-switching motion. The associated time coefficients were analyzed in the frequency domain to obtain its oscillating frequency. In addition, the shear forces on the pipe were computed and the power spectral density function was calculated as well. It was found that the swirl-switching has the same frequency as that of the oscillation of the shear stress. This represents strong evidence that the swirl-switching could possibly be the cause of the erosion by contributing the oscillating shear forces on the pipe.

This work was performed in support of the U.S. Department of Energy's Fossil Energy Crosscutting Technology Research and Advanced Turbine Programs. The research was executed through NETL Research and Innovation Center's Advanced Alloy Development Field Work Proposal. This research was supported in part by an appointment (X.H.) to the NETL Research Participation Program sponsored by the U.S. Department of Energy and administered by the Oak Ridge Institute for Science and Education. The use of the NETL's supercomputer Joule is also acknowledged. This research was partially supported by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research, Environmental System Science (ESS) Program. This contribution originates from the River Corridor Scientific Focus Area (SFA) project at Pacific Northwest National Laboratory (PNNL).

This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.