Critical Appraisal of Integrated Computational Fluid Dynamics/Surface Roughness Models for Additive Manufactured Swirl Burners

Additive manufacturing (AM) has been identified as a disruptive technology that enables the creation of complex geometries and structures that were previously impractical to manufacture using conventional methods. This novel technology is increasingly being recognized in the gas turbine industry as it offers numerous benefits including time-efficiency, cost-effectiveness, and unprecedented potential for improving the use of renewable low- and zero-carbon fuels like hydrogen, ammonia, and biofuels [1]. For instance, Ansaldo has developed a new sequential burner, called the Center Body Burner, for implementation into the GT36 H-class gas turbine using AM technology. The new burner surpassed state-of-the-art hardware regarding emission reduction, fuel flexibility, and load flexibility [2]. Other additive solutions for industrial gas turbines (GTs) are rapid prototyping, on-site repair service, developing advanced-cooling structures, and mass production [3].

The adoption of AM technology in the gas turbine industry has been limited due to a few challenges, despite its numerous benefit [4]. One of the major challenges is the rigorous design requirement for gas turbine parts, which demands a comprehensive understanding of AM process and material properties. AM processes produce typically higher surface roughness compared to conventional processes due to the layer-upon-layer manufacturing technique and the complex nature of particle deposition and fusion [5]. A review study demonstrated that rough surfaces can significantly affect the flow and heat transfer by modulating boundary layer flows [6]. This can compromise the aerodynamic efficiency, structural integrity, and overall performance of gas turbine parts. Another study has reviewed the interactions of turbulent flows with rough surfaces, highlighting roughness-induced effects of increased pressure drop, induced boundary layer transition, and enhanced heat transfer [7]. A series of experiments were conducted on both AM and traditionally machined swirler inserts in a representative gas turbine combustor [8]. The findings demonstrated the modification of mean velocity, turbulence statistics, and NOx emissions with surface roughness height. To gain a thorough understanding of the rough wall-turbulent flow interactions involved, it is essential to conduct comprehensive investigations using both experimental and numerical methods. While experiments are crucial, computational fluid dynamics (CFD) can provide more detailed information, particularly of interest in cases where experiments are not feasible.

Computational fluid dynamics is widely used for designing and optimizing gas turbine components, employing a range of turbulence models for predicting the key features of heat and flow transfer within reasonable accuracy. Among all other models, Reynolds-averaged Navier–Stokes (RANS) models are still being adopted in industrial applications due to its relatively low computing power requirements and ease of use. Many researchers have demonstrated good agreement with experimental data to predict swirling flow structures, typically encountered in GTs, within the RANS approach [9–11]. However, surface roughness adds further complexity and uncertainty to numerical flow simulations due to the variations of roughness geometry and scale in near-wall regions [7]. In the literature, three main approaches consider surface roughness effects in CFD models, each with limitations and requirements.

The first approach is to modify the boundary condition on the walls to ensure the downward shift in the logarithmic velocity due to roughness elements [7,12]. This approach is based on the modification of the standard law-of-the-wall for smooth surfaces. Many researchers have used this approach with RANS turbulence models in a variety of applications [13,14]. However, this method has its challenges as it assumes a correlation between measured surface roughness parameters (e.g., the measured peak-to-valley roughness heights, R_z) and equivalent sand-grain roughness (k_s). The lack of a universal correlation makes it difficult to accurately apply this method, despite many proposed correlations [15,16].

The second approach is the “discrete element” approach that has shown promise in overcoming these limitations [6]. This approach introduces an extra term into the governing equations to account for the flow restriction caused by surface roughness, as well as the drag and heat transfer on roughness elements [17]. One of the advantages of this approach is that it is not correlated to the Reynolds analogy, making it applicable to both uniform and nonuniform surface roughness [18]. However, it is not well-suited for use in three-dimensional unsteady flow fields, which has hindered its use in GT-related flows.

The third approach involves fully resolving surface roughness within the computational grid, which theoretically offers the ultimate way to investigate the effects of surface roughness. However, the computational requirements of the simulation domain often limit the applicability of this method due to a high ratio between the associated geometry and roughness length scales [19,20].

In order to meet the design requirements of the GT combustors made through AM, it is very important to predict turbulent flow-rough wall interactions. There is a clear need to develop reliable and robust models that require less computational demand but can still accurately predict AM-induced surface roughness affects. In this study, two different roughness approaches were compared within the RANS framework applied to an unconfined, atmospheric premixed burner with different AM swirl inserts. The paper describes the process of applying wall roughness approaches to CFD simulations. The study conducted a mesh independence analysis to ensure that the results were not dependent on the grid. Finally, the paper discusses the ability of the selected approaches to predict the effects of roughness elements on the mean characteristics of swirling flows.

The CFD simulations were performed for a swirling premixed burner equipped with AM swirl inserts of different surface roughness heights. The study used a commercial software ansysfluentv.2023.r1. Two wall roughness modeling approaches were compared and validated in this study: the modified law-of-the-wall and roughness resolving approaches. The RANS approach with the realizable k–ɛ closure model was used to predict the time-averaged motions of turbulent swirling flows in the computational fluid domain. This model is widely used in the research studies of turbulent swirling flows with a good prediction of measured velocity profiles [21,22]. Scalable wall function and enhanced wall treatment were selected based on the wall roughness approaches adopted and its requirements for wall-bounded turbulent flows. In order to maintain the boundary layer mesh entirely within the log-law region and avoid the singularity issues arising from finer mesh for the modified law-of-the-wall approach, the scalable wall function was used in the CFD simulations. Moreover, enhanced wall treatment was applied to the wall-resolved RANS simulations, to ensure the resolution of the viscous layer on the rough surfaces and the application of the wall functions to the rest.

Cardiff University Gas Turbine Research Center’s high pressure optical combustor used in the numerical simulations houses a swirl burner that consists of a modular solid body with radial-tangential inserts giving a geometrical swirl number of 0.8, as shown in Fig. 1.

The terms A_noz and A_tan refer to the exit area of the burner nozzle and the area of the tangential inlet, respectively. The variables r_tan and r_noz represent the effective radius of the tangential inlet and the radius of the burner exit nozzle, respectively. Additionally, Q_tan indicates the tangential flowrate, while Q_tot signifies the total flowrate.

The setup involves the use of turbulent swirling flows, which emerge from the swirl inserts and then stabilize on an annular bluff body with an outer diameter of 18 mm. The flow then expands into the rig through a nozzle of 40 mm inner diameter. Nine swirl vanes, aligned in tangential and radial configurations, are used to impart the swirling flow into the airflow. This setup is commonly known as the generic swirl burner and has been widely used in many studies before [23,24].

To simulate the flow dynamics in the rig, the computational domain was constructed based on the assumption that flow is unconfined as the confinement ratio is low at 0.14. This simplification has been implemented to ensure a smooth and efficient mesh generation and reduce allocated computational time. Figure 2 presents a three-dimensional computational domain of the unconfined generic swirl burner with the unstructured mesh built-in. It consists of a plenum chamber feeding ambient air into the burner via two inlets, a swirl burner and an annular fluid volume with three outlets. In Fig. 2(a), X₁ represents the spatial location 5 mm above the burner exit, where CFD data were validated against experimental data.

A grid independency study was conducted with three tetrahedral mesh sizes (ΔX) (3 mm, 2.25 mm, and 1.5 mm), giving a total number of cells ranging from 0.87 × 106 to 3.7 × 106, as shown in Table 1.

To determine the ideal mesh size for the simulations, the axial velocity profile was used as a representative parameter and compared to that obtained from the experimental study at X₁ = 5 mm [8]. The mesh independency test results are provided for the modified law-of-the-wall approach applied to the 8M swirler as shown in Fig. 3.

The mean axial velocity profile of the swirling flow is predicted well by the realizable turbulence model, which matches the peak velocities and recirculation zones of the measured values. Due to the minimum discrepancy between the experimental and simulation results, a grid size of 1.5 mm was selected for performing all other CFD cases.

For this research, a numerical model was used to examine the performance of three swirl inserts that were fabricated using AM and had different surface finishes. These inserts were previously identified as “8R” (raw AM swirler), “8G” (grit blasted AM swirler), and “8M” (traditionally machined swirler) in a study conducted by Runyon et al. [8]. For each swirl burner and its five separate surfaces, R_z values were measured by averaging the ten-point surface roughness, tabulated in Table 2. Detailed information on characterization of surface roughness can be found in Ref. [8]. For all simulation cases, Table 2 values were used for k_s input.

The second approach resolves roughness elements within the computational grid by using enhanced wall treatment. Figure 4 shows the geometrical representation of the surface roughness elements extracted upon the smooth wall surfaces of the swirl inserts.

The roughness elements were aligned uniformly across the wall surfaces within a distance, so-called target distance shown in Table 2. The target distance was set to obtain maximum achievable density of roughness elements, considering the allocated computational capacity.

where u_i and u_i′ are the mean and fluctuating velocity components, respectively, t is time, p is pressure, and μ is the dynamic viscosity.

where μ_t is the turbulent viscosity, k is the turbulent kinetic energy, and δ_ij is the Kronecker delta tensor.

In comparison to other k–ɛ turbulence models, the realizable k–ɛ model uses a variable $Cμ$ proposed by Reynolds.

where$ks$ is the physical roughness height or so-called as sand-grain roughness height.

Various roughness functions have been reviewed in literature [7]. Ansys Fluent adopts Cebeci and Bradshaw formulations [29] based on Nikuradse data, which calculates $fr$ for each of three distinctive roughness regimes: hydraulically smooth, transitionally rough, and fully rough regime.

• For the hydrodynamically smooth regime (⁠$ks+≤2.25)$

  $ΔB=0$

  (8)
• For the transitional regime (⁠$2.25<ks+≤90$⁠)

  $ΔB=1κln[ks+−2.2587.75+Csks+]x sin {0.4258(lnks+−0.811)}$

  (9)

where C_s is a roughness constant. For tightly packed, uniform sand-grain roughness, C_s = 0.5. Higher values imply the departure from the uniform sand-grain roughness. In this study, C_s were set to 1 for the wall boundary conditions of modified law-of-the-wall models.

• For the full rough regime (⁠$ks+>90$⁠)

  $ΔB=1κln(1+Csks+)$

  (10)

The boundary conditions were chosen to match those of the experimental study in Ref. [8], in order to confirm the numerical accuracy of the physical model. For each inlet, a mass flow boundary condition was used with a prescribed flowrate of 0.00805 kg/s and an air temperature of 573 K, which corresponds to an equivalence ratio of 0.55 for a methane-air mixture. At the inlets, turbulence intensity and hydraulic diameter were set to 4.72% and 0.02 m, respectively. Pressure outlet boundary conditions were applied at the outlets, with turbulence intensity set to 10% and hydraulic diameters specified for each outlet. The wall domains were assigned a no-slip wall boundary condition, and the temperature was set to 573 K.

The solution has been calculated using the governing equations of three-dimensional, incompressible flow inside the burner, and realizable k–ɛ turbulence model equations were discretized over the computational cells and iteratively solved by using the software. The pressure-based coupled algorithm for pressure–velocity coupling, second-order upwind scheme for spatial discretization of the governing equations, and Green-Gauss node for evaluation of gradients and derivatives. PRESTO! interpolation scheme was applied to the model for calculating pressure values at the cell faces as it performs well with high Reynolds flows and high swirling flows [25,31]. In the numerical model, the convergence criteria were met by monitoring the axial flow velocity component, especially at locations with significant velocity gradients. Additionally, the residuals of the governing equations were required to have an absolute convergence criterion of 10⁻⁴. For faster convergence, the global time-step formulation for the pseudotime method was used and the time scale factor was set initially to 10⁻⁴. It was gradually increased once the solution stabilized and converged smoothly.

This section focuses on the results from radial locations at a fixed downstream location from the nozzle exit (X = 0, Y = 0, and Z = 0), as shown in Fig. 2(a).

Figures 5 and 6 show the predicted data for the 8M, 8G, and 8R swirl inserts. The predictions are based on the modified law-of-the-wall approach using equivalent sand grain roughness height (ESGR), and roughness resolving (resolved) approaches. The experimental data [8] are also included for comparison (denoted as “Exp”).

In the rough case (8R), both the ESGR and resolved approaches overestimate the experimental peak velocity by 7.2% and 4.2%, respectively. The discrepancy from the experimental data becomes more noticeable in the steepest shear layers, particularly between Y = 9–20 mm and Y = 20–24 mm. In the positive steepest shear layer (Y = 9–20 mm), both methods predict the velocity values, with the resolved model exhibiting an average discrepancy of 9.2%, while the ESGR model shows a larger average discrepancy of 13.8%. Conversely, in the negative steepest shear layer (Y = 20–24 mm), the ESGR method outperforms the resolved model, achieving an average discrepancy of 14.5% from the experimental data. In contrast, the resolved model displays a much higher average discrepancy of 32.8%.

The results indicate that the ESGR based approach struggles to accurately predict the mean velocity shift with relative roughness height in the positive shear layer. This is likely to be due to the low accuracy of the correlation used to estimate the equivalent sand-grain roughness height. On the other hand, the roughness resolving approach predicts the velocity variation with the relative roughness height reasonably well. In terms of computational expense time, both methodologies have a similar average time per iteration for similar mesh size. For the 8R case, the resolved method has an averaged computational time of 34 s per iteration while the ESGR method demands 30 s per iteration.

The dimensionless roughness height, k_s+, was calculated for each surface of the 8M, 8G, and 8R swirlers based on the sand-grain roughness height. The results are displayed in a contour map, given in Fig. 7.

The values range from 0.04 to 5.31, indicating the presence of both smooth and transitional rough regimes for 8M, 8G, and 8R swirlers, as defined by Cebeci and Bradshaw [29]. Note that the lower (k_s+ smooth) and upper bands (k_s+ rough) for the onset of transitionally rough and fully rough regimes were varied in the literature [7]. Additionally, there is no available experimental data to confirm whether the boundary layer remains in a transitionally rough regime on the wall surfaces of the swirl inserts. When it comes to estimating the roughness of sand-grain surfaces, a single correlation parameter like roughness height is not enough to provide accurate results. Studies [16,32] have shown that more complex correlations, such as those that take into account multiple roughness parameters (e.g., skewness function and effective slope), are needed to accurately estimate the roughness height for realistic surfaces. It is important to consider the three-dimensional topology of the rough surfaces in order to get a more precise estimation. The uncertainty in the correlation estimating the sand-grain roughness height could be the reason why the model fails to detect the shift in velocity (Fig. 6(a)) as the surface height changes.

It has been already established [33] that the presence of surface roughness above the admissible level tends to intensify the wall shear stress and the thickness of the turbulent boundary layer. The extent of this impact varies with the scale of the roughness [33]. Figure 8 shows the local variation of the wall shear stress, on the rough surfaces of the 8M, 8G, and 8R swirlers. A comparison was made for the relative roughness height and selected rough surface approaches.

Based on the contour images, the modified law-of-the-wall approach does not indicate any significant changes in τ_w concerning wall roughness height. This approach calculates the maximum values of τ_w at around 5.8 and 6.2 Pa for the 8M and 8R swirlers, respectively, which differ by only 7%. However, the geometrically resolving wall roughness approach predicts that τ_w is almost twice as high for all the swirlers, with a 13% variation in the maximum values.

where $Gz$ is the axial flux of swirl momentum, $Gx$ is the axial flux of axial momentum, R is the radius, $Vx¯$ and $Vz¯$ are axial and tangential velocity component of the flow. In order to study the change in swirl number along the length of the nozzle, the maximum value of swirl number was calculated in the Y direction at every 2 mm interval in the X direction, starting from X = −22.0 mm and ending at X = −2.0 mm.

Figure 9 shows the calculated maximum swirl number variation along the X direction for 8M, 8G, and 8R swirl burners. Within the nozzle, the swirl number deviates locally from the geometric swirl number of 0.8 for both approaches and all cases, varying along the streamwise direction. The modified law-of-the-wall approach overlaps almost entirely the swirl number for all swirler inserts, indicating that surface roughness has no influence on the axial and tangential velocities. On the other hand, the geometrically resolved wall approach clearly predicts the local variation of the swirl number with surface roughness height. At X = −2 mm, the 8R swirler produces lower swirling than the 8M swirler, indicating a change in recirculation zone topology inside the nozzle. This could well explain the inward velocity shift with surface roughness.

The recirculation zone topology was examined for the 8M and 8R swirlers, which have significantly different maximum swirl numbers. The comparison was also made for the selected approaches. The isoprofiles of the axial velocity at zero were drawn to visualize the recirculation zone, as shown in Fig. 10.

As predicted, the modified law-of-the-wall approach demonstrates no variation in the central recirculation zone with the roughness height, in both the X and Y directions. On the other hand, the geometrically resolving roughness approach predicts the shrinking of the recirculation zone with the surface roughness height, resulting in significant shortening in the lengthwise direction. The research study on high-swirl combustion has [35] uncovered a relationship between NOx emissions and the residence time within the recirculation zone. The particle image velocimetry results have indicated that a low-swirl injector has a weaker and smaller recirculation zone, which traps a smaller recirculating mass and has a shorter residence time compared to a high-swirl injector. This results in 60% less NOx produced by the low-swirl injector. The shrinking of the recirculation zone would reduce the residence time and thus NOx emissions. This has been observed in Ref. [8] that an increase in surface roughness leads to a reduction in NOx emissions, even when the adiabatic flame temperature and exhaust gas temperatures are similar.

Overall, the roughness-resolving approach has shown better performance as it uses the enhanced wall treatment, addressing the near-wall zones in swirling flows. This is done by smoothly blending the linear and logarithmic law-of-the-wall, while also accounting for the impact of pressure gradients that are commonly encountered in swirling flows. It is important to note that the predictive capability of the roughness resolving method can be further improved by increasing the number of roughness structures and thus the frequency of the height of roughness to represent the texture of rough surfaces better [35].

This study assessed the predictive capabilities of two common roughness modeling strategies within the RANS approach: the resolved and ESGR. The CFD simulations were carried out for the AM generic swirl burners of different surface textures and validated against published experimental data. Both modeling strategies demand similar computational expense. The results demonstrate that the roughness-resolving model provides better agreement with the experimental data, which predicts the velocity variation with roughness height. Nevertheless, both methods reveal a more noticeable discrepancy from the experimental data in the steepest shear layers. The mean flow field analysis shows that surface roughness shortens the recirculation zone, which can impact flame stability and NOx emissions of fuels, to be appraised in subsequent studies.

This work was supported by the UKRI Industrial Decarbonisation Research and Innovation Centre (IDRIC). Rashed Al-ajmi gratefully acknowledges the receipt of a scholarship from the PAAET of Kuwait during his sabbatical leave and the support given by the Gas Turbine Research Center (GTRC) of Cardiff University.

• UKRI Industrial Decarbonisation Research and Innovation Centre (IDRIC) (Grant No. EP/V027050/1).

The authors attest that all data for this study are included in the paper.

A_noz =

burner exit nozzle area (m²)

A_tan =

swirler tangential inlet area (m²)

C_s =

roughness constant

Exp =

experimental data

k =

turbulent kinetic energy

K_s =

equivalent sand grain roughness height (μm)

p =

pressure

Q_tan =

swirler tangential volumetric flow rate (m³/s)

Q_tot =

burner exit nozzle volumetric flow rate (m³/s)

R_a =

arithmetic average surface roughness (μm)

R_q =

RMS surface roughness (μm)

R_z =

ten-point mean surface roughness (μm)

r_noz =

burner exit nozzle radius (m)

r_tan =

swirler effective radius of tangential inlet (m)

resolved =

roughness resolving

S_g =

geometric swirl number

t =

time

u_i and u_i^′ =

mean and fluctuating velocity components

8G =

grit blasted AM swirler

8M =

machined swirler

8R =

raw AM swirler

=

turbulent dissipation rate

μ =

dynamic viscosity

μ_t =

turbulent viscosity

AM =

additive manufacturing

CFD =

computational fluid dynamics

ESGR =

equivalent sand grain roughness height

GTs =

gas turbines