Comparison of turbulence modelling approaches in simulation of a feather shuttle: a porous conical bluff body

The aerodynamics of a feather shuttlecock, a porous conical bluff body, are unique in comparison to other sports projectiles. Experimental wind tunnel studies have been published that present values of drag coefficient (Cd) for traditional feather shuttles that vary widely (0.48 < Cd < 0.74). It is difficult to compare published experimental data, due to a lack of clarity concerning experimental apparatus. All studies have used traditional sting mounts inserted aft of the shuttle base, and it is believed this has a strong influence on Cd, as significant air movement is known to occur in this region. Flow passes through gaps formed by individual feather shafts, or rachis, inserted into the shuttle base. The use of computational fluid dynamic (CFD) simulation in the analysis of shuttles has great potential as analysis can be performed without the need of an experimental sting. This study presents the first CFD simulations of a geometrically realistic feather shuttle. Careful consideration must be given to applied grids, numeric, and turbulence models (unsteady RANS vs scale resolving) if results obtained are to be reliable. CFD results present detailed insights of shuttle aerodynamics, and the significance of flow passage between the feather rachis and internal to the shuttle. The study raises significant concerns regarding the appropriateness of rear sting mounts in shuttle wind tunnel experiments.


Introduction
Badminton shuttles are incomparable in form and mass to other sports projectiles. The aerodynamics and flight behaviour of the shuttle are unique. Aerodynamically classed as bluff bodies (a porous conical skirt preceded by a solid hemispherical dome), shuttles are categorised as either traditional (feather) or synthetic (plastic). Traditional shuttles are widely considered to exhibit superior performance [1].
A traditional shuttle utilises 16 individual waterfowl feathers, typically goose or duck, to form the conical shuttle skirt. Having been removed from the wing, feathers are first cleaned and graded for quality. Feather barbs are trimmed and stripped from the base of the feather shaft, rachis, to form a desired shape to the vane. This is then inserted into a fabric-covered predrilled cork base, aligned, and glued into position. The stripped portion of the rachis is bound together with twine, and both rachis and twine lacquered. Finally, shuttles are left to fully cure before being graded on flight characteristics (see Fig. 1).
In comparison, a synthetic shuttle skirt is a polymer injection-moulded component (typically a grade of nylon). The form is complex, designed as a porous thin walled structure, to recreate the performance of traditional feather variants. For decades, manufacturers have attempted to create a synthetic shuttle with comparable performance to a feather shuttle; however, this is yet to be achieved. Consequently, the aerodynamic behaviour of shuttles has been of interest to researchers for many years.
Studies have been published detailing wind tunnel investigations of feather shuttles [1][2][3][4][5][6], for which aerodynamic drag coefficients have been collated and are presented in Table 1. The studies report a large range of drag coefficients. Some of this variance is to be expected as feathers are a natural product. It is also known that wind tunnels can have a significant effect on the measured values of drag coefficients; however, many of these studies use different designs of experimental sting (some are designed to permit free rotation of the shuttle in the longitudinal axis) and wind tunnel configuration yet do not provide details such as equipment dimensions or turbulence values.
The influence of wind tunnel arrangement is seen in the work of Alam et al. [2] who compared the performance of five different feather shuttles. All shuttles exhibited fluctuating C d , with a significant reduction in reported values (14-30%) below Re = 96,500, a trait not seen in any other published study. It transpired in a recent paper [5] that the study had been conducted in a large industrial tunnel, and the behaviour was attributed to high tunnel turbulence values (T i > 2%). Using a low turbulence tunnel, Alam et al. [5] presented C d values comparable to Cooke [7].
The work of Cooke [7] presents the lowest published coefficient C d ~ 0.48 for a freely rotating feather shuttle. At the other extreme, Cohen et al. [6] reported C d ~ 0.68 for a freely rotating shuttle, with increased C d ~ 0.74 when statically fixed. This influence of rotation on measured C d was also investigated by Kitta et al. [4]. In contrast to Cohen et al. [6], the study revealed a freely rotating shuttle had a marginally higher C d ~ 0.61, than a statically fixed shuttle C d ~ 0.56. Both studies concluded that rotation had no significant influence on drag [1].  Hoerner [8] presents widely accepted drag coefficients of solid conical bluff bodies. For a conical half angle = 20° (similar to a shuttle), C d ~ 0.43 is reported. The drag coefficients for a porous shuttle as reported in [5,7] are therefore considered to be low; it has been shown through wind tunnel experimentation that the introduction of porosity to a conical bluff body results in a considerable increase in measured drag coefficient [10].
It appears that the key feature of a shuttle with substantial influence on C d , are the gaps between the feather rachis. Kitta et al. [4] confirmed the importance of this region by sealing the rachis cage with tape, and measured values of C d ~ 0.4. A shear layer and vortex roll up at the trailing edge were evident for both shuttles; however, the presence of rachis gaps was demonstrated to increase vortex size, driven and augmented by axial jet flow. The work was extended by Hasegawa et al. [9] who presented an internal PIV visualisation, albeit using a cutaway shuttle where the skirt normal to the viewing direction on either side had been removed, compromising flow structure. This showed the axial jet to increase the overall width and extent of the wake structures. These experimental studies only considered the rachis in a pure open or closed state.
A detailed study of the influence of rachis gap size was conducted by Lin [10], utilising wind tunnels, and computational fluid dynamics (CFD). Using a generic shuttle geometry, with plain smooth conical skirt, a series of prototypes of varying rachis gap size were created. The generic approach was justified as it would not be possible to vary the rachis gap in a feather shuttle; yet this feature will ultimately depend upon the section of rachis used, how well it is stripped, and diameter of twine. Prototypes were wind tunnel tested, and corresponding simulations used to visualise flow structures. An "optimum" gap size was revealed where C d was maximum. The influence of gap on C d was explained through plots of pressure coefficient internal and external to the skirt, demonstrating an increase in pressure differential with the introduction of gaps, thus C d increased with gap size, in agreement with [1,4,7].
What is evident from the preceding studies [1,4,9,10] is the flow through the rachis region of a shuttlecock appears to have significant influence on aerodynamic performance. Flow in this region is influenced by the gap size between rachis [10] and blockage is shown to significantly affect measured drag. Yet all experimental studies detailed in Table 1 used traditional sting mounts inserted into the aft of the shuttle base though this region. It is hypothesised the introduction of a sting may create a blockage, compromising internal flow, and affect the measured drag coefficient. The variation in experimentally reported drag, could therefore be attributable to this method of support.
CFD simulation negates the need for a sting mount and has been used to provide clarity of flow in this region [10,11], albeit with generic or simplified geometries. These studies utilised steady Reynolds-averaged Navier-Stokes (RANS) CFD simulation, which can provide a starting point for understanding the complicated flow structures associated with shuttles, but which do not fully reveal the complex underlying flow structure. This was demonstrated by Hart [12] who presented an analysis comparing simulation of a synthetic shuttle using unsteady RANS, and scale-resolving simulation (SRS). A detailed CFD analysis of flow through a geometrically realistic feather shuttle is yet to be presented. It is therefore the intention of this paper to compare a feather shuttle utilising both unsteady RANS (k-ω SST) and SRS [delayed detached eddy simulation (DDES)]. The shuttle geometry was modelled with an open and closed rachis cage, to explore the significance of this region on flow structure and performance.

Shuttle geometry
The shuttle geometry was based on a Yonex Aerosensa 50, a readily available traditional style cork and goose feather shuttle. Geometry was created in ProE Wildfire 5, using a combination of traditional and laser scanning metrology, as shown in Fig. 2. Shuttle height was H = 86 mm and maximum skirt diameter D = 66 mm. It is accepted practice to use D (characteristic length scale for frontal area drag coefficients [1]).
A single feather was removed from the shuttle skirt and scanned in isolation, (Fig. 2c). Capture of an entire shuttle geometry would have been unfeasible, due to occlusion from vane overlap. Scanned data in a tessellated format provided the basis for construction of a geometrically clean NURBS (non-uniform rational basis spline) CAD representation of the feather. The NURBS feather was inserted into a CAD representation of the cork base, at determined baseline angles = 20° and γ = 12° (defined in Fig. 2) and patterned radially with overlap to form the 16 feathers of the skirt. At γ = 12° a slight gap results, approximately 0.3 mm, between adjacent vanes (internal vane edge, and adjacent vane inner surface). Twine binding the individual rachis, with a cord diameter of 1.4 mm, was then added completing the geometry. A shuttle geometry with a closed rachis was also created, as shown in Fig. 2d. This was achieved by revolving a thin surface through this region.

Boundary conditions
The shuttle was placed in a cylindrical computational domain of diameter 9D, with the rotational axis of shuttle and domain aligned. Velocity flow inlet was placed 4H upstream of the shuttle and pressure flow outlet 12H downstream. Domain outer walls were modelled with a slip condition. The shuttle was modelled in a non-rotating condition, and porosity of the feather vane was not considered. The inlet flow velocity was varied between 10 m/s < U ∞ < 60 m/s (45,000 < Re < 271,000). A turbulence intensity T i = 1%, with eddy viscosity ratio T VR = 2, was specified at all flow boundaries.

Optimisation of computational approach
The study used ANSYS Fluent v15. Two approaches to modelling turbulence phenomena were investigated. Unsteady RANS modelling using a k-ω SST model (with inclusion of the Wilcox low Reynolds number terms [13]), and SRS using DDES. It should be noted that although DDES provides scaled content, it is still coupled with an underlying k-ω SST model with low Re terms for near wall modelling. Governing equations and turbulence models were solved using non-iterative time advancement. This was implemented using a fractional step pressure velocity coupling, and a second-order temporal discretisation. Comparable second-order finite volume spatial discretisation was used for all simulations. A time step of T s = l min /U where l min is minimum cell length, and U is free stream velocity, was used for advancing solutions.

Computational grids
A series of computational grids were constructed to ascertain the sensitivity of the numerical solution to grid resolution. Reynolds number was fixed at Re = 90,000. All grids were hybrid in construction and used full geometry. A tetrahedral surface grid was constructed over the shuttle and flow domain boundaries. Two alternative surface resolutions for the shuttle were tested: 213 k elements (element size range 0.4-1.0 mm), 838 k elements (element size range 0.2-0.5 mm). The largest elements were placed over the shuttle base, the minimum size of element was required to capture local rachis and vane curvature. Hanging node hexahedral meshing was used in the construction of volumetric grid. A series of volumetric grids were created that investigated: the systematic refinement of the immediate region around and behind the shuttle, the inclusion of prismatic elements upon wall boundaries. Refinement of a volume grid was achieved through specification of a maximum cell length through and behind the shuttle. The influence of near wall prismatic elements, on the resolution of surface boundary layers, and the viscous sub-layer in simulation, was observed. Ten layers of prismatic boundary layer elements were constructed over the entire shuttle surface. Simulations were run until flow structures had established in full, prior to activation of data sampling to obtain time averaged values. Details of constructed grids, and associated time averaged data, are provided in Table 2.
At a fixed surface resolution, what appears as a "grid independent" solution can be obtained through systematic volumetric mesh refinement. For example, refinement of Grid A (2.75 million elements) to Grid D (10.55 million elements) reduces the predicted C d by 4.3%. This is a comparable mesh error to that accepted by Verma [11], and Lin [10], who had claimed mesh independency using only 3-4 million unstructured elements. Both, however, had reached this conclusion through the comparison of only two grids [10,11], and after conducting only 100 iterations [10].
Data presented are time averaged as the shuttle form was found to produce transient flow structures with even the coarsest mesh, C d fluctuations of ~ 1%. Such fluctuations were evident in the figures of Verma [11] yet could not have been fully resolved due to the implementation of a steady state RANS solution. Although volumetric mesh refinement increased the clarity of the predicted flow structure, this had little influence on instantaneous C d fluctuations (Grids A-H).
Introduction of near wall prismatic elements (Grids I-L) had a significant influence on resolved flow features, transient behaviour, and time-averaged C d . This reduced values of wall y + (normalised distance between wall surface and the first adjacent computational cell centre), and increased the number of elements spanning the wall adjacent viscous sub-layer, improving capture of the adjacent boundary flow gradient. Flow separations and subsequent reattachments on both base and vane surfaces were sharply captured. As a result, instantaneous fluctuations of C d increased to 4%.
Based upon convergence of C d and resolved flow structure, parameters as applied in the construction of Grid K (16.9 million cells) were chosen as the basis for construction of all subsequent grids used in simulation. This is comparable to the grid size required (21 million cells) to ensure sufficient scale resolution in simulation of synthetic shuttles [12]. Slightly higher overall volume cell count results in simulations of synthetic skirts. This is due to the preservation of the smaller scales of a perforated hole pattern, in comparison to the gaps between rachis and vanes in a feather shuttle form.  Fig. 3. These values were obtained by averaging the instantaneous drag over the normalised period 10 = T s U/D. Average C d values obtained from wind tunnel investigations [6,7,9] are provided for comparison, as are values obtained in simulation by Verma [11].
The baseline shuttle geometry had a time-averaged drag coefficient of C d ≈ 0.72. Instantaneous values fluctuated by 4%, due to the prediction of large vortex structures, that will be shown through flow visualisation. These fluctuations are shown in Fig. 4, which provides an example of the temporal history of the drag force coefficient. C d is normalised by C dRMS (C d root mean square) for comparison purposes. These fluctuations are sufficient in magnitude to necessitate the use of a time-stepped solution, to ensure a correct solution convergence. The obtained C d  is comparable in magnitude to Cohen [6]; however, it is higher than data reported by other researchers.
Closing the rachis cage reduces time averaged drag by 37% to C d = 0.45, which is comparable in magnitude to Hasegawa [9] who reported a 32% reduction, C d = 0.38. As previously stated, interpolated drag of a solid conical bluff body of = 20°, which is comparable to a shuttle, is C d ≈ 0.43 [8]. The obtained results for a shuttle with a closed rachis are therefore comparable to those of a conical bluff body. What is also notable is the reduction in the rate of vortex shedding when the rachis is closed, Fig. 4. This disappears completely when simulations are conducted using a k-ω SST turbulence model.

Component force contribution
A breakdown of the force contribution for individual shuttle components is shown in Fig. 5. Pressure drag is dominant, accounting for approximately 93% of measured C d . The measured contribution of each component was found to be independent of Re. The feather vanes are the largest source of drag, 57%, followed by the rachis cage 22%, base 13%, and the twine used in binding the rachis 8%. In contrast closing the rachis cage is seen to have a significant influence on individual contribution of components. A reduction in the measured drag is observed on all components of the shuttle, with the exception of the rachis. Drag on the rachis rises slightly which is attributable to the increase in surface area caused by closing this region.
The individual force contributions can be explained through the investigation of variables, such as mean pressure coefficient (C p ) and shear stress (τ ), acting upon the surface of the geometry. Figure 6 provides visualisation of surface contours of C p and τ x (shear stress aligned with axial flow direction) acting on the shuttle base, with circumferentially averaged values plotted in Fig. 7. As C d and statistical surface values were found to be independent of Re, all results plotted herein are for Re = 90,000.
Regardless of whether the rachis is open or closed, flow stagnation occurs on the tip of the hemispherical section as expected. However, a subsequent boundary layer separation, that occurs as fluid passes over this surface, is evidently influenced by the downstream rachis region.
A lower C p and greater τ x is evident on the open rachis base in comparison to the closed rachis. This suggests a faster moving flow, that is drawn in through the rachis, delaying separation and inducing reattachment. Boundary layer separation occurs at 89° on the dome, with a subsequent reattachment positioned 9 mm downstream on the cylindrical portion of the base. In contrast, closure of the rachis causes an earlier separation, 82.9°, comparable to that of a laminar smooth sphere separation, 82.5° [14], and flow remains detached. Comparing C p acting centred on the aft base surface, C p = − 0.4 for the closed rachis, and C p = − 1.2 for the open rachis, a result of no fluid passage within the shuttle for the closed rachis. Measured drag on the shuttle base is therefore higher for the open rachis condition.
Differences in mean C p between the shuttle skirt internal and external surfaces are shown in Figs. 8 and 9. For clarity, a single feather is shown due to the overlap between adjacent vanes. Both turbulence models produce comparable mean pressure distributions. A lower differential between the inner and outer vane surface is evident when the rachis is closed, C p values on the skirt internally and externally are lower. A positive C p is visible on the external surfaces of each feather vane and rachis. Towards the outer trailing edge of each vane C p ≤ − 0.45 attributable to flow separation and wake formation behind the  shuttle. A pressure recovery can be seen to occur between sections of overlapping vanes, with a localised increase in C p . The open rachis shuttle exhibits a slightly higher external C p in comparison to the closed rachis, however, a significant region of difference is seen in the vanes' leading edges. External edge C p = 1 at this point, with an opposing region of C p ≈ − 3.3 on the internal vane surface. This is indicative of flow passing through the rachis cage and separating internally at the vane's leading edge. In contrast C p ≈ − 0.4 throughout the closed rachis shuttle inner, as was measured on the aft base surface.

Flow visualisation
The complexities of the flow field and vortex structure surrounding, and internal to the shuttle, are clearly visualised in CFD. These are unaffected by the presence of a wind tunnel sting, or by the modification of shuttle geometry to visualise internal structure [9]. The influence of closing the rachis is further visible, as is the influence of the applied turbulence model approach in capture of the flow structure.
The mean flow velocity field is shown in Fig. 10. Both turbulence models produce comparable averaged flow fields. When the rachis region is closed, the flow internal to the shuttle is predominantly stationary. In contrast, localised flow accelerations can be observed through the open rachis, with an axial core of circulating flow located aft of the base, and a large ring of circulation aft of the vanes. Further small circulating pockets of flow reside behind the twine regions. Passage of flow through the rachis, and the lesser fluid movement between the vanes, shape these circulating flows, Fig. 11. This shows streamwise cross-sections of time-averaged normalised swirl (normalised by maximum swirl velocity).
A series of counter rotating annular flows can be seen throughout the shuttle. The rotational directions are driven by the cross-sectional shape of the rachis, vane curvature, and insertion angle γ. Significant and complex swirling flow therefore exists within the open rachis shuttle. Of note, aft of the base the core of circulating flow exhibits a strong anticlockwise rotation, that is seen to extend to the vanes. It is directly into this region that a wind tunnel sting would be inserted. The full complexity of the internal flow is revealed in Fig. 12 by visualising the instantaneous flow structures.
Centre plane contours of instantaneous vorticity (applicate Z direction) are shown in Fig. 12a-d. The formation of shear layers and flow separations across the shuttle base, as indicated in Figs. 6 and 7 are clearly visible, and this flow structure is drawn through the open rachis. Transient vortex structures, driven by the swirling flow field, form within the skirt as discussed. The resolved strength and number of these are dependent on the applied turbulence approach. The k-ω SST simulations reveal large-scale single mode vortex structure only. To obtain detail across a range of turbulence scales requires a scale-resolving model such as DDES. In contrast with no flow passing through the closed rachis, whereas a weak internal structure is predicted by the DDES model (attributable to passage of fluid between adjacent feather vanes), k-ω SST fails to resolve any internal structure.
Differences between the capabilities of the applied turbulence approach is further emphasised in the prediction of external wake structures. Three-dimensional structures are visualised using the Q-criterion method [15] Fig. 12e-h. Considering the open rachis geometry, large vortex structures are observed that separate at the vane trailing edge. Formation of these structures begins at the vane leading edge as thin streamwise cores that roll up between adjacent vane surfaces. These thin cores remain distinct a short distance after separation from the trailing edge, before merging to form larger-scale vortex rings that roll up in a Kelvin-Helmholtz pattern. DDES clearly predicts greater scale content, however, k-ω SST can be seen capable of comparable prediction of these largerscale structures.
In contrast, the k-ω SST model again fails to predict any significant transient structure for the closed rachis. The trailing edge wake forms as a steady vortex sheet, explaining the constant value of measured C d , Fig. 4. DDES however predicts significant vortex structure that varies both spatially and temporally over the entire shuttle surface, however, no clear Kelvin-Helmholtz behaviour is observed due to the lack of internal flow. This is in agreement with the PIV observations of Hasegawa [9] who noted reduced core flow as delaying vortex roll up at the trailing edge.

Summary and conclusions
Computational fluid dynamic simulations were presented of a geometrically accurate feather shuttle. This study has demonstrated the importance of fluid flow between the shuttle rachis, on measured variables including coefficient of drag, through the formation of complex flow structures internal to the shuttle. Careful consideration in the application of computational grids, and of turbulence model through the consideration of unsteady RANS (k-ω SST) and SRS (DDES), was demonstrated. Sensitivity of obtained results to the applied computational grid has been shown, as has the importance of time-dependent simulation. This is necessary to capture the transient behaviours and flow structures that are to be expected with such a geometry, a notable failure of previously published studies. The following main conclusions were obtained: • The detail of resolved flow structure is dependent on applied turbulence model; fully transient behaviours were only captured with the application of scale-resolving simulation, particularly when flow was restricted between the rachis. This has important implications for the simulation of conical bluff bodies with restricted flow. • Fluid movement between the feather rachis create complex flow structures internal to the shuttle with direct influence on measured variables. This raises questions regarding the insertion of traditional sting mounts aft of the shuttle base in experimental wind tunnel studies. Insertion of a sting in this region will cause interference compromising measured drag; restriction of flow in the rachis region is known to reduce measured drag. This interference of the flow in an important region may explain some of the large variation in measured drag reported in literature. • The use of CFD in the analysis of shuttles therefore has great potential due to the omission of an experimental sting. However, as stated previously simulations must be carefully assembled if reliable results are to be obtained. • It is recommended that future wind tunnel studies should dramatically reduce the size of the experimental sting or consider different mount locations. Alternatively, direct trajectory analysis methods should be explored.