CFD simulations of cyclist aerodynamics

Computational Fluid Dynamics (CFD) allows researchers and practitioners to analyze cyclist aerodynamics and identify areas for improvement. Despite the numerous CFD simulations of cyclist aerodynamics in the scientific literature, the extent to which user choices in the large number of computational parameters affect the simulation results remains largely unexplored. This paper aims to establish a set of best practice guidelines for CFD simulations of an isolated cyclist in time trial position through a systematic and comprehensive sensitivity analysis. It includes the computational grid in terms of surface, near-wall, and far-field volume grids and the turbulence modeling. The study reveals a high sensitivity of the computed drag area to the surface grid resolution and y + value, while the impact of the growth rate and the grid resolution in the wake is relatively smaller. The results emphasize the importance of complete reporting of grid characteristics and the need for grid-sensitivity analyses, and provide prioritization of the key parameters for such analyses. Satisfactory agreement with wind tunnel measurements is achieved using Scale-adaptive simulations (SAS) and steady RANS with the Transition SST (T-SST) or SST k-ω turbulence model for closure. This work intends to contribute to accurate and reliable CFD simulations of cycling aerodynamics


Introduction
In professional cycling, where every second counts, minimizing aerodynamic drag is crucial for optimal performance.At the high speeds typically observed in professional cycling and on level road, aerodynamic drag is the dominant force acting on the cyclist (Kyle and Burke, 1984;Grappe et al., 1997;Martin et al., 1998;Lukes et al., 2005).Different methods can be employed to optimize rider positions, bicycle design, and equipment choices in terms of aerodynamic drag: wind tunnel testing, on road and track testing, and numerical simulation with Computational Fluid Dynamics (CFD).
Over the past 15 years, CFD has been increasingly applied in cycling aerodynamics research as outlined in several review papers (Crouch et al., 2017;Malizia & Blocken, 2020a, 2021).It enables researchers and practitioners to analyze the aerodynamic forces acting on the cyclist, identify areas for improvement and quantitatively evaluate such improvements.In the scientific literature, CFD has been applied to various aspects of cycling aerodynamics, including rider position optimization, bicycle component analysis, drafting and group riding, crosswind effects, and clothing and equipment optimization.
Rider position optimization involves the use of CFD simulations to analyze the airflow around different rider positions.By investigating the effects of body posture (Defraeye et al., 2010a(Defraeye et al., , 2011;;Griffith et al., 2014;Fintelman et al., 2015;Blocken et al., 2018aBlocken et al., , 2019;;Giljarhus et al., 2020), such as the angle of the back, arms, and legs, on aerodynamic drag, optimal positions can be identified that minimize drag and improve performance.
In addition, CFD simulations help analyze the influence of crosswinds on cycling aerodynamics (Fintelman et al., 2015).By simulating the interaction between the wind and the cyclist's body, bicycle, helmets (Beaumont et al., 2018;Giappino et al., 2019) and other equipment (Malizia et al., 2021), CFD guides the development of aerodynamic designs that minimize the negative impact of (cross)winds, reduce drag and enhance overall performance.CFD simulations can be very sensitive to the choices in the large number of parameters that need to be set by the user.Considering that a 2% variation in drag area can result in approximately a 0.5-s difference per kilometer ridden (Blocken et al., 2021), and recognizing that in a cycling time trial, seconds, and sometimes tenths of a second, can be decisive between winning and losing, it is crucial to assess the extent to which these parameters affect the simulation results.The CFD studies mentioned in the previous paragraphs employed a wide range of computational parameters.An overview of these settings is provided in Section 2 of this paper.These computational parameters include the computational grid, the turbulence modeling approach and the turbulence model.Therefore, CFD simulations should ideally be performed on grids obtained by grid-sensitivity analyses and the simulation results should always be validated with experimental data, which can be wind tunnel measurements or dedicated track and field tests (Franke et al., 2007;Tominaga et al., 2008;Blocken, 2015;Crouch et al., 2017;Malizia & Blocken, 2020a, 2021).Defraeye et al. (2010b) performed a first attempt to evaluate the performance of different turbulence models.They found that the RANS SST k-ω model exhibited the best overall performance among the turbulence models tested.Large-Eddy Simulation (LES) was also found to perform well, albeit with substantially higher computational costs.The study further indicated a preference for low-Reynolds number modeling (LRNM) over wall functions for boundary layer modeling.It is important to note that this study was conducted over a decade ago at a time when computational resources were more limited compared to today.As a result, the study did not include the bicycle, and the grids were generated based on grid generation guidelines rather than grid-sensitivity analyses, and contained fewer than 8 million cells.In more recent times, aided by rapid developments in both software and hardware, more sophisticated turbulence modeling approaches have been applied with a detailed inclusion of the geometry of the bicycle and its components in simulations.Grids can now be constructed based on thorough grid-sensitivity analyses, and the cell counts for a single cyclist commonly range from 30 to 50 million cells (Griffith et al., 2019;Wang et al., 2022;van Druenen and Blocken, 2023).
The cumbersome task of performing detailed grid-sensitivity analyses, testing a range of turbulence models and comparing the results with measurements can be alleviated by best practice guidelines.Such guidelines provide basic advice for the generation of high-quality grids and for the selection of the best performing turbulence modeling approaches and turbulence models.They serve as a starting point for future CFD simulations where they can reduce the number of grids tested in the grid-sensitivity analysis and the number of turbulence models tested in the validation exercise, and substantially reduce the time and computational cost of the simulation exercise.These best practice guidelines themselves need to be based on systematic and comprehensive sensitivity analysis.
Existing guidelines on the impact of computational parameters regarding the use of CFD in cycling aerodynamics are limited to RANS simulations and to the analysis of an isolated cyclist's body (Defraeye et al., 2010b) or bicycle wheels (Malizia et al., 2019;Malizia and Blocken, 2020b).RANS studies focusing on full configurations of one or more cyclists often apply general CFD guidelines (e.g., Casey and Wintergerste, 2000;Menter et al., 2021) and those available for studies in automotive (e.g., Lanfrit, 2005) or urban environments (Franke et al., 2007;Tominaga et al., 2008;Blocken, 2015).For instance for bicycle wheels, it was found that at least six surface cells should be present along the shortest spoke edge, an average y + value below 3.5 should be applied, and the cell growth in the wake should be maximum 1.15 (Malizia et al., 2019).In automotive applications, a maximum growth rate of 1.2 for cell length was advised, with finer cells recommended near sharp edges and regions with high curvature, while larger cells could be used on planar surfaces.The prismatic cells in the surface boundary layer should be aligned with the local flow, and a maximum prismatic growth rate of 1.2 should be applied between consecutive layers.Local mesh refinement is strongly advised in the car's wake region (Lanfrit, 2005).Guidelines for RANS simulations of urban flows (Franke et al., 2007;Tominaga et al., 2008;Blocken, 2015) recommend a domain size of 5H between the modeled area and the inlet, sides, and top of the domain, and a size of minimum 10H between the modeled area and the outlet of the domain, where H is the height of the highest building.The cell length growth rate should be smaller than 1.3, and at least 10 cells should be applied per building side.
To the best of our knowledge, a systematic and comprehensive investigation of the impact of the computational grid on the simulation results in cycling aerodynamics, considering a full cyclist configuration, has not yet been performed, while such an investigation on the impact of the turbulence modeling approach and the turbulence model has not been performed beyond the early work of Defraeye et al. (2010b).This highlights the necessity for such investigations to develop best practice guidelines that in turn should support accurate and reliable CFD simulations of cyclist aerodynamics.

state of the art
Table 1 presents an overview of previous peer-reviewed CFD studies of cycling aerodynamics.This overview only includes 3D CFD simulation studies that have undergone validation.The table specifies the turbulence modelling approach, turbulence model and near-wall treatment used in each study, as well as the grid typology and its main characteristics, including surface cell size, prism layer information, and wake cell size.These parameters will be discussed further in the following subsections.
The turbulence models used to close the 3D RANS equations were the standard k-ε turbulence model, the Shear stress transport (SST) k-ω turbulence model and the Transition Shear Stress Transport (SST) k-ω turbulence model.
In general, two approaches have been used to model the near wall region: wall functions and low-Reynolds number modeling (LRNM).In the latter case, the grid is fine enough to directly resolve the boundary layer including the thin viscous sublayer.This can provide more accurate predictions of boundary layer transition and flow separation (Casey and Wintergerste, 2000;Menter et al., 2021).Most studies listed in Table 1 have used a LRNM approach.Other studies have applied wall functions and sometimes found a fair agreement with measurements (Blocken and Toparlar, 2015;Blocken et al., 2016Blocken et al., , 2018c;;Beaumont et al., 2018).It was however stated that this was expected to be due to the specific shape and position of the cyclist and might not work out well for other cyclists or cyclist positions (Blocken et al., 2018c).

Computational grid
Important grid characteristics for accurate and reliable CFD T. van Druenen and B. Blocken simulations include grid resolution, grid typology, cell type, cell skewness and aspect ratio (Casey and Wintergerste, 2000).A full description of the grid characteristics is important to assess the quality of the grid.However, important grid characteristics are rarely reported in a complete way.

Grid typology
Table 1 shows that different grid typologies have been used in previous studies.Basically, three types can be identified: (i) structured grids; (ii) unstructured grids; and (iii) hybrid grids.Structured grids are regular or Cartesian grids that consist of rectangular and/or cubic hexahedral cells.Structured grids are relatively easy to generate and provide good accuracy for simple geometries.However, they can be difficult to generate for complex geometries and can require a large number of cells to achieve a good resolution.(ii) Unstructured grids are irregular grids that consist of cells with arbitrary shapes, such as tetrahedral, polyhedral or prismatic cells.Unstructured grids are more flexible than structured grids and can be generated for complex geometries more easily.However, they can be more difficult in terms of grid quality control and may require more cells to achieve a good resolution.(ii) Hybrid grids are grids that combine structured and unstructured cells to achieve a balance between accuracy, flexibility and economy, i.e. total cell count.
In cycling aerodynamics research, due to the complex geometry of the cyclist and bicycle, the majority of studies applied hybrid grids, consisting of prismatic cells in the boundary layer on the body surfaces, tetrahedral cells in the vicinity of the cyclist outside the boundary layer and hexahedral cells in the far field.However, more exceptionally, also grids that are mainly polyhedral (Oggiano et al., 2016), hexahedral (Giljarhus et al., 2020) or tetrahedral (Blocken et al., 2013(Blocken et al., , 2018c;;Beaumont et al., 2018) have been employed.

Surface grid
The surface grid resolution should allow for the accurate representation of local flow features, such as flow attachment, separation, reattachment, and areas of high flow curvature.In addition, the grid should be fine enough to capture areas with high pressure gradients on the cyclist.However, only four of the 22 studies listed in Table 1 reported the minimum and maximum surface cell size on the body surfaces, ranging from 0.5 to 5 mm (Giappino et al., 2019;van Druenen & Blocken, 2021, 2023;Wang et al., 2022).

Near-wall volume grid
The grid in the near wall region is important for accurately capturing the near wall flow behavior.The dimensionless wall unit y + is generally used to assess the grid resolution in this region.It represents the dimensionless distance of the wall-adjacent cell centroid to the wall.While some studies report the y* value instead, in the current study the -not provided.a When assessing different turbulence modeling approaches, average y + values were 43 for RANS, 5.2 for DES, and 3.4 for LES.
b Generally lower than 1, locally max 5. c Growth rate defined by last-ratio approach.
T. van Druenen and B. Blocken y + value will be used.It is defined as where ν is the kinematic viscosity (m 2 /s) and u τ the friction velocity (m/ s) τ w ρ √ with τ w the wall shear stress (Pa) and ρ the density (kg/m 3 ).
The required value of y + depends on the type of near-wall treatment.For the wall-function approach, the recommended range of y + is typically 20 to about 300 (Casey and Wintergerste, 2000).A lower y + value, typically less than 5, is recommended for LRNM to ensure that there are at least a few cell layers inside the viscous sublayer and the buffer layer and that the flow in these layers is adequately resolved.
The dimensionless wall distance y + (or y*) is reported in most studies in Table 1.In studies applying LRNM this value was generally lower than 5.In studies with wall functions it ranged from about 21 to 334.In addition, many studies reported the height of the cell center of the walladjacent cell y p .The transition from this low-height high-aspect ratio cell to the larger cells further away from the body is usually performed by a number of cell layers containing prismatic cells with increasing height per layer.This prism-layer growth rate and the number of prism cell layers are important factors considering that the prism cell layers should cover the entire boundary layer region (Casey and Wintergerste, 2000;Menter et al., 2021).Guidelines recommend a maximum growth rate of 1.1 (ANSYS, 2021a) or 1.15 (Gerasimov, 2016) and a minimum of 5-10 cell layers in the boundary layer up to y + = 20, meaning that the total number of layers could go up to 60 (Casey and Wintergerste, 2000).The growth rate between successive layers has been reported in only some studies and ranged from 1.05 to 1.25 (Table 1).The number of prism layers has been reported more often and ranged from 10 to 40 layers.

Far-field volume grid
Sufficient grid resolution in the wake region is essential for accurately capturing flow features such as shear layers, recirculation zones and reattachment (Menter et al., 2021).In many cases the wake region interacts with the boundary layers of nearby objects.Some studies employed a mesh refinement zone in the anticipated wake regions, often with a predetermined maximum cell length (Oggiano et al., 2016;Wang et al., 2022).In contrast, others adopted a gradual increase in cell size from the body in all directions, without any specific refinement in the wake zone (Blocken et al., 2018a(Blocken et al., , 2019;;Giljarhus et al., 2020).The reported cell sizes in the vicinity of the body or within the wake zone seem to indicate a consensus among studies with a maximum cell size of approximately 0.03 m.

CFD simulations: reference case
The starting point for the sensitivity analyses in the present paper is termed the reference case.This section provides the information about the computational geometry and grid, boundary conditions, and solver settings for the reference case.In the subsequent Sections 4 and 5, the reference case is systematically altered in terms of specific computational parameters for further analysis.These parameters are the grid resolution including the surface, near-wall and far-field volume grids, the grid cell type, the turbulence modeling approach (RANS and SAS) and the turbulence model.After each step in the sensitivity analysis, an appropriate setting for the computational parameter under study is chosen, and it becomes the basis for the analyses of the other parameters.Finally, the obtained simulation results are compared with wind tunnel measurements.

Computational geometry
The study focuses on a time-trial setup.A time-trial bicycle with open front and disk rear wheel is selected.The bicycle geometry is simplified by neglecting chain, cables and spokes.The impact of these simplifications on the total drag is considered small.In a similar configuration, Malizia et al. (2021) determined that the non-rotating wheel's spokes contributed approximately 0.7% to the total drag.Furthermore, Malizia et al. (2019) previously assessed the impact of computational settings on the drag of an isolated wheel, emphasizing the need for a minimum of six surface cells along the shortest spoke edge and an average y + value below 3.5.The geometries of wheels, bicycle, cyclist and helmet are obtained by an Eva structured light 3D scanner and processed to yield the computational geometry shown in Fig. 1.The scanning procedure and data processing received rider consent and was approved by the Ethical Review Board of Eindhoven University of Technology with Nr.ERB2020BE_1859456_WT.

Computational domain
The size of the computational domain is determined by following best practices guidelines (Franke et al., 2007;Tominaga et al., 2008;Blocken, 2015) to ensure minimal interference with the flow field around the cyclist (Fig. 2).The dimensions are 16 m (width) x 32 m (length) x 8 m (height).With a cyclist and bike frontal area (A) of 0.338 m 2 , the resulting blockage ratio is 0.27%.This value is below the recommended maximum of 3% (Franke et al., 2007;Tominaga et al., 2008;Blocken, 2015).The cyclist is positioned 6 m away from the inlet of the domain, while maintaining a 0.05 m gap between the wheels and the bottom surface of the domain (Malizia and Blocken, 2020b).

Computational grids
A sensitivity analysis will be performed to investigate the impact of the grid typology, the surface grid, the near-wall volume grid and the remainder of the volume grid on the simulation results.All sensitivity analyses are performed for two grid typologies; (i) poly-hexcore and (ii) tetrahedral-based (tet-based).A poly-hexcore grid typology connects six-sided-base prism cells in the boundary layer to octree hexahedral cells in the volume grid by polyhedral cells.In the tet-based grid typology triangular prism cells with a three-sided base are situated in the boundary layer, while the volume grid comprises tetrahedral cells.A detailed sensitivity analysis concerning the number of surface cells is provided in Subsection 4.1.
The near-wall grid for the bicycle and cyclist body surfaces consists of a series of prism layers down to the viscous sublayer.This is required because the LRNM approach is employed.The sublayer is the near-wall region where the dimensionless normal-to-wall distance (y + ) is less than 5.The impact of y + is evaluated in Subsections 4.2.1 and 4.2.2.
The near-wall grid layers need to encompass the entire boundary layer (Casey and Wintergerste, 2000;Menter et al., 2021) and are constructed using the last-ratio technique (ANSYS, 2021b), employing a specified last-ratio value of 30%.This technique ensures a smooth transition between the prism cells in the final near-wall grid layer and the tetrahedral/polyhedral cells in the volume grid.The impact of the growth rate is examined in Subsections 4.2.3 and 4.2.4.In order to evaluate the influence of the volume grid farther away from the body, different maximum cell lengths are specified in four wake zones surrounding the cyclist.More information on the specification of these wake zones and on the impact of wake cell refinement is provided in Subsection 4.3.

The post-processed scan of the cyclist and bike is imported into
Fluent Meshing and converted into an editable surface mesh.
2. The surface mesh is refined in regions of high curvature and coarsened in regions with less curvature using curvature size functions.
The minimum and maximum cell lengths are varied as explained in Subsection 4.1.3. The boundaries of the computational domain are discretized with a maximum surface cell length of 0.2 m.On the bottom boundary, 0.05 m below the bicycle's wheels, refinement is applied to avoid a large volume cell size growth rate, facilitating a smooth transition from the much smaller wheel surface cells.4. Four mesh refinement zones are defined around the cyclist and in the wake, where grids with maximum cell lengths as described in Subsection 4.3 are applied.5.The near-wall region is meshed with thin layers of prismatic cells.
The properties of these layers are systematically adjusted, as outlined in Subsection 4.2 and 4.3.Grids with different y P (y + ), different growth rates and different numbers of prism layers are used.6.The volume mesh is generated using a maximum cell length of 0.2 m.
In the tet-based grid, the cell lengths are governed by a fixed maximum growth rate of 1.1.For the cubic cell lengths in the polyhexcore typology, a nonconformal 1:2 transition is applied at the interfaces of the grid refinement zones.

Computational settings
The simulations are performed using ANSYS Fluent 21R1 (ANSYS, 2021c).For the evaluation of the impact of the computational grid (Section 4) the 3D steady RANS equations are solved with the shear stress transport (SST) k-ω turbulence model (Menter, 1994).For the evaluation of the impact of turbulence modeling approaches and turbulence models (Section 5) the 3D steady RANS equations are solved with the transition SST (T-SST) (Menter et al., 2006a(Menter et al., , 2006b)), SST k-ω, realizable k-ε (RLZ k-ε) (Shih et al., 1995) and Spalart-Allmaras (SA) (Spalart and Allmaras, 1992) turbulence model and the 3D SAS approach (Menter and Egorov, 2010;Egorov et al., 2010) are coupled with the SST k-ω turbulence model.Though several studies used wall functions to model the flow in the near-wall region, it is known that this approach has difficulties in accurately predicting the location of the separation points.In the assessment of drag a correct prediction of separation regions is crucial.Wall functions assume that the flow near the wall behaves in a logarithmic manner.This assumption is based on the logarithmic law of the wall, which is valid for fully developed turbulent flows over flat walls (Casey and Wintergerste, 2000).However, over complex bodies and curved surfaces, the flow experiences  additional effects, such as adverse pressure gradients and centrifugal forces, which can significantly alter the mean velocity and turbulence profiles near the wall.Wall functions fail to accurately capture these complexities (Casey and Wintergerste, 2000) and are therefore more likely to fail in predicting the point of flow separation.Therefore, in this study, the wall-function approach is not considered and the flow is resolved all the way down to the viscous sublayer by solving the approximate form of the governing equations in every cell of the very fine mesh near the wall to capture the full boundary layer profile.
For all simulations pressure-velocity coupling is performed with the Coupled algorithm and the Least Squares Cell Based method is used to compute gradients.Pressure interpolation is second order.Second order upwind discretization schemes are used for both the convection and the viscous terms of the governing equations in the 3D RANS simulations, while for the SAS approach the momentum equation discretization is performed with central differencing.For the 3D RANS simulations, the pseudo-transient solver is applied with automated time step sizing.Following an initial 1000 iterations used for flow field initialization, the drag area values are averaged over 5000 iterations to achieve a statistically steady average, with a deviation of less than 0.1% compared to a run with twice the number of iterations.For the SAS simulations, temporal discretization is performed by the bounded second-order implicit scheme and the time step size is determined by maintaining a Courant-Friedrichs-Lewy (CFL) number below one in wake region W 1 (the cell zone surrounding the cyclist, as detailed in Subsection 4.3).The CFL number relates the time step size (Δt) to the streamwise cell size (Δx) and the streamwise velocity (U s ), which is equal to the cycling velocity (U c ) in the absence of head-, tail-, or crosswinds, through the equation CFL = U s Δt/Δx.In this study, with a maximum cell length (Δx) of 0.01 m in wake zone W 1 , a cycling velocity of 15 m/s and a conservative CFL choice of 0.75, the SAS simulations are conducted using a time step of Δt = 0.0005 s.

Boundary conditions
At the computational domain's inlet, for the 3D RANS simulations, a uniform mean velocity of 15 m/s and a turbulence intensity of 0.5% are imposed, simulating high-speed cycling in the absence of external atmospheric wind.The decay of the turbulence intensity from the inlet to the location of the cyclist was negligible.For the SAS simulations comparable conditions are applied, with no perturbations imposed on the inlet flow.For all simulations, except those in Subsections 4.2.2 and 4.2.4,and those using the RLZ k-ε model (which doesn't support such modifications, see Subsection 5.2), the cyclist's body surfaces are modeled as no-slip surfaces with an equivalent sand-grain roughness height of k S = 0.1 mm (Blocken et al., 2018b(Blocken et al., , 2020;;van Druenen & Blocken, 2021, 2023).The bicycle and helmet surfaces are modeled as smooth no-slip walls.At the outlet, zero static gage pressure condition is specified.The slip wall boundary condition is implemented at the bottom, top, and lateral boundaries.

Impact of grid resolution
A sensitivity analysis is conducted to examine the influence of the surface grid, the near-wall volume grid (prism layers), and the far-field volume grid on the simulation results.All simulations apply the k-ω SST turbulence model for closure of the 3D RANS equations.The results are evaluated based on the drag area C D A and contours of the velocity ratio U r and the surface pressure coefficient C p are included to highlight differences in the flow field.The drag area is defined by in which C D is the drag coefficient (− ), A the frontal area (m 2 ), F D the drag force (N) and ρ the air density (kg/m 3 ).The velocity ratio U r is defined as U r = U/U c with U the magnitude of the local 3D velocity vector.The pressure coefficient is defined as with p the local static pressure and p 0 the reference static pressure.

Impact of the surface grid
The impact of the surface grid is assessed using 9 poly-hexcore and tet-based grids (Fig. 3).The distribution of surface grid cells was determined using curvature size functions (ANSYS, 2021b) that assign cell sizes based on the local curvature of the surface.In all these functions, the surface cell growth rate is set to 1.1, and the curvature normal angle is set to 14 • .The minimum cell lengths range from 1 to 4 mm, while the maximum cell lengths range from 2 to 16 mm.The near-wall grid layers are constructed using the last-ratio technique, employing a specified last-ratio value of 30%.To maintain controlled growth, the number of near-wall grid layers is determined in a manner that the growth rate between successive layers never exceeds 1.15.For consistency across all grids in this subsection, the y P value is uniformly set to 0.01 mm, resulting in a mean y + value of about 0.5.However, it should be noted that due to flow acceleration, locally higher y + values may be observed.
The selection of different values for the minimum and maximum surface cell size significantly affects the overall surface and volume cell count, as indicated in Table 2. Comparing the grids with hexagonal versus triangular surface cells, the former clearly results in a much lower number of surfaces cells because of the distinctive shape of these surface cells.This difference is also evident in the total volume cell count.
The results from the surface grid refinement study are presented in Fig. 4 in terms of drag area values.Percentage differences are calculated relative to the drag area obtained on the finest grid of the same typology; i.e. grid (1-2) in Fig. 3 and Table 2.Note that results for the two coarsest poly-hexcore grids are not included due to inadequate convergence of the simulations.
Both the poly-hexcore and tet-based grids exhibit a similar trend, with the drag area values finally converging towards the value obtained on the finest grid.This convergence is influenced by both the minimum and maximum cell lengths.Specifically, with a minimum cell length of mm, the differences range from − 2.3% to − 3.8%.The poly-hexcore grids, configuration 4-8 excluded, demonstrate a monotonically converging trend as the maximum cell length decreases, whereas this is not the case for the tet-based grid.However, with a minimum cell length of 1 mm, the results on both grids monotonically converge towards the drag area value on the finest grid as the maximum cell size decreases.
Fig. 5a and h display the contours of C p for the finest grid, for the poly-hexcore and tet-based grids, respectively.In the frontal view, distinct areas of overpressure can be observed on the front-facing surfaces.Additionally, there are regions of under-pressure visible on the sides, with the most prominent ones located on the upper arms and the legs.Some clear differences between both grids appear at the side of the torso and the upper arms.
Fig. 5b-g and Fig. 5i-n illustrate the differences in C p compared to the values obtained on the finest grid.In this representation, blue areas indicate an underestimation of C p for the corresponding grid, while red areas indicate regions where C p is overestimated.An important region with regard to flow separation is the upper arms.It is evident, especially for the tet-based grid typology, that as the grid resolution decreases, the differences in C p become more pronounced in this area.This is likely a consequence of a less accurate prediction of separation points, which, in turn, could lead to a misprediction of the wake flow.This misprediction is also reflected by the presence of red areas on the frontal surface of the upper legs, directly in the wake of the upper arms.The airflow patterns around the cyclist's waist and torso are also influenced by the grid resolution.These complex airflows arise due to pressure interactions between the overpressure below the cyclist's waist and the underpressure above its back.Regarding the bicycle surfaces, the differences in pressure distributions on the front fork become more pronounced as the grid resolution decreases for both grid typologies.
With a deviation within 1% compared to the finest grids (Fig. 4), the grids with a minimum cell length of 1 mm and a maximum cell length of 3 mm are considered to be a good compromise between accuracy and economy and are selected for the remainder of this study.

Impact of the near-wall volume grid
The impact of the near-wall volume grid was evaluated in two steps.First, in Subsections 4.2.1 and 4.2.2, the distance of the cell center point y P (m) to the wall is varied to study the impact of varying y + .Next, the use of varying growth rates in successive near-wall grid layers is investigated in Subsections 4.2.3 and 4.2.4.

Impact of y +
The impact of y + is studied for both the poly-hexcore and tet-based grid topology, for which 12 grids are created in total (Fig. 6).For all grids, surface grid "1-3" from Subsection 4.1 with a minimum and a maximum surface cell length of 1 and 3 mm, respectively, is used.The value of y p is varied as follows: 0.0025, 0.005, 0.01, 0.02, 0.04 mm and 0.08 mm.The corresponding y + values (y + mean ) are derived from the actual simulations, taking into account time-and area-averaging: for the above-mentioned y p values, y + mean is about 0.125, 0.25, 0.5, 1, 2 and 4. Larger y p values are not included in this study, as y + mean values larger than 5 represent y p values that exceed the viscous sublayer.The number of cells increases with reducing y + value, as more prism layers are required to ensure a maximum growth rate between successive layers of 1.15 and a last-ratio value of 30%, as shown in Table 3.
The values of y + mean at different parts of the cyclist and bicycle are listed in Table 4.The drag area results are reported in Fig. 7. Compared to the grid with the smallest y + value (y + mean ≈ 0.125), the drag area values obtained on the poly-hexcore grids exhibit a monotonically converging trend as y p is reduced.There is limited additional value in using mean y + values lower than 1, as the drag area results remain within a 0.3% range.However, when y + mean exceeds 1 the drag area starts to significantly deviate from that of the finest grid, with deviations of − 1.2% (y + mean ≈ 2) and − 2.1% (y + mean ≈ 4).In contrast, the drag area results obtained on the tet-based grids do not exhibit a monotonically converging trend as the y + mean value decreases.An outlier is observed at y + mean ≈ 0.5, where the drag area obtained escapes the expected trend of increasing drag area with decreasing y + mean .Furthermore, the deviations from the finest grid are relatively large.For instance, the drag area obtained on the grid with y + mean ≈ 1 is 2.6% lower than that obtained on the finest grid (y + mean ≈ 0.125).
The drag area convergence observed in Fig. 7 for the poly-hexcore grids is less clearly reflected in the contours of C p -C p-fine in Fig. 8. Multiple dark red and blue areas are present around the upper arms, legs, waist, and lower back.Interestingly, these differences are less pronounced in the tet-based grids, with smaller regions and lighter shades of red and blue.It is possible that the differences in underpressure and overpressure balance each other out in the poly-hexcore grids, resulting in similar total drag area values.Additionally, it could be argued that the application of roughness influences this evaluation.Therefore, additional simulations are performed in which the body surfaces are modeled as smooth no-slip walls.

Impact of y + on smooth surfaces
In this Subsection, the methodology described in Subsection 4.2.1 is employed, but with the cyclist body surfaces being modeled as smooth, no-slip walls.
The drag area results, depicted in Fig. 9, exhibit similar trends for both grid typologies.Compared to the finest grid, maximum deviations of 0.7% (poly-hexcore) and 0.2% are observed across all grids with a y + mean up to 1.As y + mean increases, the deviations from the finest grid become more pronounced, reaching − 2.1% for the poly-hexcore grid   and − 4.9% for the tet-based grids at y + mean ≈ 4. In contrast to the cases with roughness, the drag area results in Fig. 9 now align with the contours of C p -C p-fine in Fig. 10.For both grid typologies, only regions with very light shades of red and blue are observed for y + mean values up to 1, indicating minimal differences in C p compared to the finest grid.However, as y + mean exceeds 1, regions with darker shades of red and blue begin to emerge, which correspond to the increased divergence in drag area observed in Fig. 9.
Based on the above, the values of y P = 0.01 mm and y + mean = 0.5 are used for the grids in the remaining (sub)sections.

Impact of prism-layer growth rate
The impact of the growth rate is studied for both the poly-hexcore and tet-based grid typology, for which 12 grids are created in total (Fig. 11) with growth rates of 1.05, 1.1, 1.15, 1.2, 1.3 and 1.5.Surface grid "1-3" from Subsection 4.1 with a minimum and maximum surface cell size of 1 and 3 mm, respectively, is used for all grids.The value of y P is 0.01 mm (y + mean ≈ 0.5).With a fixed last-ratio value of 30%, the number of required prism layers, and thereby the total cell count, increases with decreasing growth rate as shown in Table 5.
With regard to the variation in growth rate, comparable results are observed across both grid typologies (see Fig. 12).In comparison to the finest grid, which has a growth rate of 1.05, the drag area remains within 1% for grids with a growth rate up to 1.5 for the tet-based grids.For the poly-hexcore grids, deviations within 1% are observed for grids with a growth rate up to 1.3.However, with a growth rate of 1.5, an underestimation of the drag area by 2% is found.

Impact of prism-layer growth rate on smooth surfaces
Considering the demonstrated sensitivity of the impact of y + mean on surface roughness as demonstrated in Subsections 4.2.1 and 4.2.2, the impact of the prism-layer growth rate is further assessed with the cyclist body surfaces modeled as smooth, no-slip walls.Fig. 13 shows the Fig. 6.Volume grid slice in centerplane z = 0 and surface grid for the poly-hexcore and tetrahedral grids with varying y + mean value.results.Similar to Fig. 9, the tet-based grid exhibits higher drag areas in the absence of roughness.For the finest grid, when the cyclist's body surfaces are assigned an equivalent sand-grain roughness height of k S = 0.1 mm, a drag area of 0.1808 m 2 is obtained (Fig. 12).In contrast, for the smooth body surface, a drag area of 0.1883 m 2 is observed (Fig. 13).However, this trend does not hold for the poly-hexcore grid, where similar drag area results are obtained regardless of body surface roughness.For example, on the finest grid, a drag area of 0.1817 m 2 is obtained with roughness, while a drag area of 0.1818 m 2 is found for the case with smooth body surfaces.Fig. 14a and d display contours of mean C p for rough body surfaces, while Fig. 14b and e give the contours of mean C p for smooth body surfaces, for the finest poly-hexcore and tet-based grids (growth rate = 1.05), respectively.Fig. 14c and f present the differences in C p between the cases with and without roughness for the finest grids.Although the poly-hexcore grids yield similar drag area values despite variations in body surface roughness, there are noticeable C p differences between the two cases (Fig. 14c).From a side view, the red areas indicate lower C p for the rough case, specifically at the left shoulder and the inner lower leg on the right side.Additionally, the front view reveals a red color surrounding a large part of the body contour.Consistent with prior studies on cylinders, the impact of surface roughness on C P is most pronounced on the cylindrical-like geometries perpendicular to the wind direction, such as the upper arms and the right lower leg (Kyle et al., 2004;Oggiano et al., 2007).Similar patterns can be observed for the tet-based grids (Fig. 14f), albeit with larger and darker shaded red areas compared to the poly-hexcore grids (Fig. 14c).Furthermore, the front view shows pronounced red areas on the inner side of the upper arms, the right lower leg and the left upper leg.The larger C p values for the rough versus the smooth cases suggest delayed flow separation, while the lower C p values can indicate earlier flow separation.In contrast to the poly-hexcore grids, the cyclist's cheek shows remarkable red areas, and a significant dark-shaded blue area was visible on the outer side of the left upper arm.
Fig. 15 depicts contours of the velocity ratio U r and 2D streamlines for the finest grids (with a growth rate of 1.05) of both grid typologies, for the cases with and without roughness.The cross-section analyzed was a horizontal plane passing through the upper left arm.Dashed lines indicate the tangents on the inside and outside of the arm, which run parallel to the main flow direction.It is expected that the introduction of roughness would result in delayed flow separation, a narrower wake, and reduced drag.
For the poly-hexcore grid in the rough case (Fig. 15a), there is indeed a noticeable delay in separation at the inner (upward) side of the upper arm, leading to a narrower wake where a low-speed recirculation zone (depicted in blue color, U r < 0.8) is located below the upper tangent line.Conversely, in the smooth case (Fig. 15b), the flow separates from the upper arm farther upstream, resulting in a wider wake that extends beyond the upper tangent line.The difference in U r (magnitude) near the separation point on the outer side of the upper arm however is less pronounced.
For the tet-based grid, the impact of roughness on the flow field around the upper left arm is significantly more noticeable.In the presence of roughness (Fig. 15c), there is a distinct delay in separation observed on the inner (upward) side of the upper arm, which is a phenomenon that is more pronounced compared to the poly-hexcore grids.This yields a lower C P for the rough case as shown by the red color in Fig. 14f (front view).In the smooth case (Fig. 15b vs. Fig.15d), the tetgrid shows a flow field similar to that of the poly-hexcore grid, with separation points occurring at comparable locations, resulting in a wake with a similar width.However, in the rough case (Fig. 15a vs. Fig.15c), the tet-based grid has a separation point on the inner (upward) side of the arm that is positioned significantly farther downstream compared to the poly-hexcore grid.As a result, the wake width experiences a substantial reduction.Additionally, the wake pattern undergoes considerable changes and shifts more towards the outer side (downwards) and is characterized by higher U r values.The lower U r values observed near the outer surface of the upper arm for the rough case (Fig. 15c) in comparison to the smooth case (Fig. 15d) contributes to the increase in C p at this location, as depicted in Fig. 14f.
For the grids in the remainder of this study, a growth rate of 1.15 is used.

Impact of the wake volume grid
The impact of the wake volume grid is evaluated by different maximum cell lengths in 4 zones: W 1 , W 2 , W 3 and W 4 , with dimensions as indicated in Fig. 16.For the poly-hexcore grid a 1:2 jump in cell length for the hexahedral grid is present at the interface between the different wake zones (non-conformal grid), while the tet-based grid is given a gradual transition for the tetrahedral cells between these zones (conformal grid).For each zone, a maximum edge length is imposed on the cubic and tetrahedral cells in the poly-hexcore and tet-based grid typologies, respectively.Table 6 provides the maximum volume cell edge lengths and the resulting total cell counts, while Fig. 17 displays the grids.As opposed to the sensitivity analyses in the previous subsections, changes in the cell lengths in the wake have a smaller impact   on the total cell count.From the finest to the coarsest grid, a reduction of 20% was achieved for the poly-hexcore grid, whereas only a 7% reduction was observed for the tet-based grid.
The resulting drag area values are given in Fig. 18.For the polyhexcore grids with a maximum cell length of 0.01 m in W 1 , no significant effect is observed when modifying the cell lengths further downstream.The maximum difference in drag area compared to the finest grid is 0.3%.Only when the maximum cell length in W 1 is increased, the deviations in drag area become substantial: increasing the maximum cell length from 0.01 m to 0.02 m and 0.04 m leads to an approximate 2% increase in drag area.
For the tet-based grids, the use of a coarser mesh in the wake region has only a minimal influence on the cyclist drag area.Similar to the previous case, the cell size in W 1 is most important.When the maximum cell lengths in W 1 range up to 0.02 m and varying cell lengths are implemented in the downstream wake zones, the drag area difference remains within 0.4% compared to the finest mesh.Only with a maximum cell length of 0.04 m in W 1 an increase in drag area of 1% is observed.A specific edge length yields a greater number of tetrahedral cells per volume compared to poly-hexcore cells.Since calculations occur at cell centers, this increased resolution of tetrahedral cells may contribute to the lesser dependency of the evaluated cell lengths W 1 in comparison to poly-hexcore grids.
Fig. 19a and l show contours of C p on the cyclist and bike surfaces and contours of U r in the vertical centerplane, for the finest poly-hexcore and tet-based grids (1112), respectively.Fig. 19b-k and Fig. 19m-v displays the difference of these quantities for the other cases/grids compared to the finest grid.For the poly-hexcore grids, consistent with the drag area results, Fig. 19b-f (W 1 = 0.01 m) shows only minimal flow field variations.Minor differences in U r are mainly observed at the interfaces between the different zones and in the immediate wake of the cyclist.However, these differences are not significantly reflected in C p differences on the cyclist and bicycle surfaces.Fig. 19g-k (0.02 m or 0.04 m in W 1 ) displays larger differences in U r as well as more pronounced variations in C p .
For the tet-based grids, differences in C p on both the cyclist and the bike are not very pronounced in Fig. 19m-v.This observation aligns with the drag area results where a maximum deviation of only 1% is observed.However, there are some clear variations in the wake flow field compared to grid 1112.Notably, these variations occur at the location where the grid is most coarsened.For instance, when a maximum cell length of 0.08 m is implemented in W 4 , a relatively large purple area with slightly lower velocity ratios becomes evident at this location.In addition, as the maximum cell length in W 2 increases, dark shaded regions of green and purple expand behind the rear wheel, originating from the wake zone interface.This implies an insufficient grid resolution to effectively capture gradients in the direct wake structures.

Impact of turbulence modeling approach and turbulence model
The impact of the turbulence modeling approach and the turbulence model is assessed and the calculated drag area is compared with the value from the wind tunnel (WT) tests.The evaluated approaches, listed from high to low computational expense, are: (i) the hybrid LES/URANS approach using Scale-adaptive simulation (SAS) coupled with the SST kω turbulence model, the steady RANS approach with (ii) the Transition SST (T-SST) turbulence model, (iii) the SST k-ω turbulence model, (iv) the Realizable k-ε (RLZ k-ε) turbulence model and (v) the Spalart-Almaras (SA) turbulence model for closure.

Wind tunnel measurements
Wind tunnel measurements were conducted in the Eindhoven opensection WT.A full-scale time-trial bicycle and a cyclist mannequin with geometry as indicated in Fig. 1 were positioned on a two-component force sensor integrated into an elevated platform as shown in Fig. 20.The platform was equipped with a customized support system.The mannequin was created by CNC milling of high-density polyurethane based on the same 3D scan of an elite cyclist as used in the CFD simulations in this study.Afterward, the mannequin received a smooth surface treatment.Also the bicycle and aero-helmet used in the WT were the same as scanned and used in the CFD simulations.
To account for the uncertainties associated with roughness modeling, three skinsuits made from different fabrics with varying levels of surface roughness were tested.The mean drag area of 0.194 m 2 , calculated from measurements on these skinsuits (0.198 m 2 , 0.197 m 2 , and 0.186 m 2 ) was used in the analysis.Fig. 21 illustrates the results where the whiskers denote a standard deviation (σ) of 0.0005 m 2 , based on 30 repeated measurements conducted within a similar experimental setup.The bicycle was equipped with aero-bars, a rear full disc wheel, and a front spoked wheel with sixteen bladed spokes.The combined frontal area of the cyclist and the bicycle was 0.338 m 2 .To determine the aerodynamic drag of the combined cyclist-bicycle system, the forces on the platform, including the support structure, were measured separately and subtracted from the total measured drag of the entire system.The measurements were performed at a reference wind speed of 15 m/s and an angle of 0 • .All results were corrected for differences in temperature, velocity, atmospheric pressure and relative humidity during the measurements.

Turbulence modeling approaches and turbulence models
The influence of turbulence modeling approaches and turbulence models is examined for both the poly-hexcore and tet-based grid configurations.For the grid, surface grid "1-3″ is selected, the chosen y P value is 0.01 mm (resulting in y + mean ≈ 0.5), 28 prism layers are employed with a growth rate of 1.15 and the maximum cell sizes in the wake zones are 0.01 m, 0.02 m, 0.02 m and 0.04 m for W 1 , W 2 , W 3 and W 4 , respectively.
This study evaluates five different turbulence modeling approaches and models, with detailed general computational settings and boundary conditions provided in Subsections 3.4 and 3.5.The application specifics are as follows: (i) For SAS with the SST k-ω model, the schemes and solver settings employed in this study follow the best practice guidelines specific to SAS (Menter, 2015).(ii) For the 3D RANS simulations with T-SST roughness correlation and production limiters (ANSYS, 2021c) are applied.(iii) For 3D RANS with the SST k-ω model the computational settings and boundary conditions are equal to those specified in Subsections 3.4 and 3.5.(iv) The boundary conditions for the RLZ k-ε  turbulence model closely resemble those outlined in Subsection 3.5.Enhanced wall treatment is applied utilizing Wolfshtein's one-equation model (Wolfshtein, 1969) in the viscous sublayer.To address the buffer layer zone towards the wall functions in the log-law region, a blending function is incorporated.However, it should be noted that in this particular combination, the application of rough walls is not supported (Ansys, 2021c).(v) For the 3D RANS simulations with the SA turbulence model vorticity-based production is applied.

Results
Fig. 21 compares the drag areas obtained using various turbulence modeling approaches and turbulence models with the WT results.The difference in drag area values obtained using the same approach and model when comparing both grid typologies is as follows: 2.8% (0.0055   m 2 ) for SAS, 0.7% (0.0013 m 2 ) for T-SST, 0.3% (0.0006 m 2 ) for SST k-ω, 0.9% (0.0013 m 2 ) for RLZ k-ε, and 1.4% (0.0023 m 2 ) for SA.The larger deviation observed with the SAS method could be attributed to its switching between URANS and LES based on the local grid characteristics.
Comparing the CFD results with the mean value of the three WT tests, a good agreement is achieved for the SAS and T-SST models.On the poly-hexcore grid, the SAS model shows a deviation of − 1.8%, while the T-SST model has a deviation of +0.9%.On the tet-based grid, the SAS model exhibits a deviation of 1.1%, while the T-SST model has a deviation of +0.3%.The SST k-ω model yields a less good agreement with deviations of − 6.3% and − 6.6% on the poly-hexcore and tet-based grids, respectively.The performance of the other methods was rather poor.The RLZ k-ε model with enhanced wall treatment underestimated the drag area by approximately 25%.In addition, the SA model underestimated the drag area by about 15%.
Fig. 22 provides further insight into the differences in C p for different turbulence modeling approaches and turbulence models by comparing the obtained values to those obtained by SAS.The drag results from the T-SST and SST k-ω turbulence models align more closely with SAS compared to the RLZ k-ε and SA turbulence models, as reflected in C p differences by lighter shades of red and blue.In the T-SST model  (Fig. 22bg), minor differences are observed in the front view, mainly on the lower part of the body.The side view reveals vertical parallel regions of red and blue on the upper arm, front fork, and lower legs, which indicate local shifts in separation points.C p differences on the rear wheel are similar for both grid typologies, while on the front wheel, these differences are more pronounced for the poly-hexcore grid.The latter observation holds true for the SST k-ω model as well (Fig. 22ch).Conversely, the front view displays larger C p differences compared to T-SST (Fig. 22ch vs. Fig.22bg).The inner sides of the upper arm show a blue shade, suggesting higher velocity and delayed flow separation, which is also visible on the lower legs.In the frontal view of the RLZ k-ε turbulence model (Fig. 22di), a blue shade delineates the upper back, legs, and both inside and outside of the upper arms.This indicates delayed flow separation, smaller wake areas, and reduced drag, supporting observations in Fig. 21.The drag underestimation by RLZ k-ε and SA turbulence models is further evident in Fig. 22deij, with large

Table 6
Grid characteristics of evaluated grids with varying maximum cell edge lengths in the wake zones.

Discussion
This paper aimed to establish a set of guidelines for conducting accurate and reliable CFD simulations of an isolated cyclist in a time-trial position.It systematically evaluated the impact of different computational parameters, namely the grid resolution including the surface, near-wall and far-field volume grids, the grid cell type, the turbulence modeling approach and turbulence model, by comparing the CFD results with each other and with the results from WT experiments.The study had several limitations that are briefly addressed: The computational geometry did not consider the presence of small bicycle components such as spokes and chain.The assumption made was that these components do not significantly impact the overall drag, and any potential effects could fall within the experimental error.However, in specific studies where these small components are of particular importance, such as in the design analysis of these components, it might be necessary to employ a locally higher cell resolution to accurately capture all the details associated with these small components.
The reference case used in this sensitivity study focused on an isolated time-trial cyclist.The authors anticipate that the proposed framework and conclusions of the study would not significantly depend on the specific posture of the cyclist.However, simulating a larger number of cyclists could ask for different computational parameters.For instance, when considering the complex nature of the flow interactions between drafting cyclists (Barry et al., 2015), the refinement of the grid in each cyclist's wake could have a more substantial impact.In such cases the application of scale-resolving turbulence approaches might be necessary due to their improved prediction of wake regions.
Although the present study primarily focused on the RANS approach, it is worth noting that scale-resolving CFD simulations, such as LES, hybrid RANS-LES, and URANS, have the potential to offer additional insights into cycling aerodynamics.Therefore, it is important to investigate whether the proposed guidelines for grid resolution in this study are also applicable and effective for these scale-resolving methods beyond the one tested in this study.
Another limitation of the present study was the absence of crosswinds, as well as the non-rotation of the bicycle wheels.The decision not to include these factors was driven by the desire to mitigate additional uncertainties and computational costs.Furthermore, it was deemed that the time-averaged effect on the wake structure resulting from cyclist pedaling would be minimal, as supported by previous research (Crouch et al., 2016;Griffith et al., 2019;Javadi, 2022).Additionally, the impact of the absence of a moving floor was considered negligible, as the measured boundary-layer height of 40 mm on the elevated platform in the WT fell well below the region of interest.The CFD simulations conducted in this study applied low turbulence level approach flow conditions, and the comparison was made with WT measurements obtained under similar low turbulence level approach flow conditions.This practice is commonly employed in both WT experiments and CFD simulations in cycling aerodynamics.However, it is important to note that low turbulence level approach flow may not always accurately represent real-world conditions (e.g. Brown et al., 2023).Therefore, future studies should investigate the impact of the computational parameters for moderate to high turbulence level approach flow conditions.Furthermore, the extensions of the computational domain were based on guidelines specific to urban wind flows.Future research should focus on investigating the impact of these extensions on CFD simulations pertaining to cycling aerodynamics.
Despite these limitations, this study provided a systematic and comprehensive sensitivity analysis leading to guidelines for accurate and reliable CFD simulations of isolated cyclist aerodynamics.
The results of the grid sensitivity study emphasize the importance of creating a high-resolution surface grid, because the drag area was found to be highly sensitive to this parameter.Using a minimum surface cell length of 2 mm in areas of high curvature, a common approach in scientific literature, led to an underestimation of the drag area by over 2% compared to grids employing a minimum cell length of 1 mm in these regions.This was true for both the poly-hex and tet-based grid typologies.Differences in C p arising from a coarser surface mesh were particularly evident in certain regions of the cyclist body and bicycle, such as the upper arms and front fork.It could be argued that local mesh refinement, involving the application of a finer curvature size function in these specific areas and a coarser curvature size function in others, could enhance accuracy while limiting the total cell count.This aspect should be further explored in future work to identify optimal mesh refinement strategies.
The sensitivity analyses on the impact of the first cell height, and thus the y + value, revealed the crucial role of this parameter.In the case of the tet-based grids, monotonic convergence of the drag area with decreasing y + mean value was not observed when considering the case that included the modeling of roughness to account for the drag-reduction capabilities of cyclist skinsuits.To address the potential influence on the evaluation due to the application of roughness, additional  simulations were performed with the body surfaces modeled as smooth no-slip walls.In these cases, monotonic convergence was achieved with a maximum deviation of 0.2% for grids with a maximum y + mean value of 1.For the poly-hexcore grid, a maximum y + mean value of 1 resulted in good agreement compared to the lowest investigated first cell height (y + mean = 0.25), irrespective of the applied surface roughness.Here, a maximum deviation of 0.7% was observed.However, with a y + mean value of 4, the deviation increased to over 2% for the poly-hexcore grids and to about 5% for the tet-based grids.
The sensitivity analyses conducted on the volume cell size in the near-wall region revealed that the impact of the prism-layer growth rate was considerably smaller than the effect of y + .Although certain guidelines suggest a maximum growth rate of 1.1 (ANSYS, 2021a) based on flow over a flat plate, it was observed that larger growth rates of 1.3 resulted in differences smaller than 1% for both grid typologies in the present study.However, note that for scale-resolving approaches such as SAS and LES, existing best practice guidelines have recommended to use 1.15 or lower growth rates (Gerasimov, 2016).
In terms of the far-field volume grid, it was shown that the cell resolution in the immediate vicinity of the cyclist can be very important.This was particularly pronounced for the poly-hexcore grid, where differences of approximately 2% were observed as a result of coarsening the grid in this specific area.On the other hand, increasing the volumetric cell length for the tet-based grid did not result in significant differences.This phenomenon is likely a consequence of the gradual transition of the tetrahedral cells in the grid and due to the fact that a specific edge length results in a higher density of tetrahedral cells per volume compared to poly-hexcore cells.Additionally, coarsening the grid further downstream did not have a substantial impact on the drag area for both grid typologies.
The sensitivity study performed on the turbulence models demonstrated that specific models achieved good to fair agreements with the drag area values obtained from the WT.In addition, these results were not strongly dependent on the grid typology used in the simulations.The SAS model showed deviations of − 1.8% (poly-hexcore) and +1.1% (tetbased), the T-SST model exhibited deviations of +0.9% (poly-hexcore) and +0.3% (tet-based), and the SST k-ω turbulence model resulted in deviations of − 6.3% (poly-hexcore) and − 6.6% (tet-based).The performance of the RLZ k-ε and SA turbulence models was less satisfactory, showing larger differences of approximately − 25% and − 14%, respectively.
It is important to note that modeling the actual roughness of cycling skinsuits is not straightforward.To take into account this uncertainty, a range of cycling skinsuits made from different textiles was tested in the WT, and the CFD results were compared to the mean WT results.Consistent with prior WT studies on cylinders (Kyle et al., 2004;Oggiano et al., 2007), applying roughness in CFD by increasing the k S value led to a downstream shift of the separation points around the upper arms, resulting in a reduction in local wake size and drag.It is worth noting that the influence of roughness on total drag and the point of flow separation depend on the grid typology, as this effect was more prominent for tetrahedral-based grids compared to poly-hexcore grids.However, accurately representing specific skinsuit roughness in CFD modeling, such as aligning it with k S or other roughness parameters, remains a significant challenge requiring further research.

Conclusions
This study investigated the effects of different computational parameters on the computed drag for an isolated time-trial cyclist.The parameters considered were the computational grid, including the surface, near-wall, far-field volume grids and cell type, the turbulence modeling approach and the turbulence model.The key findings of the sensitivity analysis are summarized as follows: Findings concerning the impact of the computational grid: • The surface cell size significantly influenced the computed drag area.
When compared to the finest grid tested, using a minimum surface cell length of 2 mm resulted in differences of at least 2%.However, with a minimum surface cell length of 1 mm and a maximum cell length of 3 mm, this difference reduced to less than 0.7% for both grid typologies.• The prism-layer growth rate had a relatively minor impact, as long as it remained equal to or lower than 1.3.Differences in drag area were within 1% for such cases.• The application of smooth cyclist surfaces with a mean y + value of maximum 1 resulted in differences of up to 0.7%.However, with a mean y + value of 4, differences increased to over 2% and about 5% for poly-hexcore and tet-based grids, respectively.It is important to note that applying roughness can increase the sensitivity to mean y + value, particularly on tet-based grids.• Regarding the far-field volume grid, close attention to the cell resolution near the cyclist is crucial.Notably, on the poly-hexcore grid, grid coarsening in this area led to approximately 2% differences.The cell sizes in the zones W 2 , W 3 and W 4 further downstream, were found to have a minimal effect on the drag area in the investigated range of 20-80 mm.
Findings concerning the impact of the turbulence modeling approaches and turbulence model: • In terms of the drag area, the SAS approach and T-SST turbulence model provided the closest agreement with the wind tunnel mean results, with deviations of up to 1.8% for both grid typologies in the case of SAS and 0.9% for T-SST.The k-ω SST turbulence model yielded a much lesser agreement, with deviations falling within 6.6%.The RLZ k-ε and SA turbulence models performed less satisfactorily, exhibiting larger deviations of around − 25% and − 14%, respectively.
The findings in this study highlight the importance of selecting appropriate computational parameters in CFD simulations of cycling aerodynamics to obtain accurate drag predictions.They offer insights and enable the development of guidelines to researchers and practitioners towards accurate and reliable CFD simulations of cycling aerodynamics.Furthermore, the study emphasizes the importance of performing grid-sensitivity analyses and prioritizes key parameters for such analyses.Based on the current cyclist and bicycle geometry and application of the SST k-ω turbulence model, the following recommendations are made concerning the impact of the computational grid: • It is recommended to prioritize sensitivity analyses related to surface cell length, with the evaluation going down to a minimum surface cell length of 1 mm.The minimum value should be applied in areas with large curvature and large expected gradients in surface pressure.• A mean y + value equal to or smaller than 1 is advised, which corresponds to a y P value of 0.02 mm for U c = 15 m/s.• The prism-layer growth rate should be kept equal to or lower than 1.3.For scale-resolving approaches, stricter guidelines might be required, such as keeping the growth rate equal to or below 1.15.• The cell size in the region immediately around the cyclist but outside the boundary layer and hence outside the prism layers, including the near wake, i.e. zone W 1 in this study, should be kept equal to or below 10 mm.Further downstream in zones W 2 , W 3 and W 4 , a maximum cell size of 80 mm can be recommended, provided that there is a gradual transition in cell size extending from W 1 to W 2 , from W 2 to W 3 and from W 3 to W 4 .This entails adhering to a maximum volumetric cell growth rate of 1.1 for tetrahedral-based grids.For poly-hexcore grids, it involves maintaining a maximum 1:2 transition in cell length.However, when studying drafting cyclists or cyclists followed by vehicles, the above-mentioned guidelines should be applied to each of the cyclists or following vehicles, which also implies smaller cells between the cyclists and potential vehicles.
Based on the current cyclist and bicycle geometry and evaluated computational settings, the following recommendations are made concerning the turbulence modeling approaches and turbulence model: In conclusion, in view of third-party basic evaluation of the grid and in view of reproducibility of CFD simulations, it is strongly advised to include the following key parameters when reporting CFD simulations of cycling aerodynamics: minimum and maximum surface cell length, y p , y + or y*, prism-layer growth rate, the number of prism-layers, cell type and the resolution in zones W 1 to W 4 , supported by some clear graphics.This is essential information for assessing the grid quality and resolution, thereby ensuring transparency and reproducibility and allowing third-party basic evaluation of CFD simulations in cycling aerodynamics.

Declaration of competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 3 .
Fig. 3. Volume grid slice in centerplane z = 0 and surface grid for the poly-hexcore and tetrahedral grids with varying minimum and maximum surface cell sizes.

Fig. 4 .
Fig. 4. C D A values obtained with varying minimum and maximum surface cell lengths.The percentage differences are with respect to the finest grid (1-2).

Fig. 5 .
Fig. 5. Contours of (a,h) mean C p for the finest grid (1-2), and (b-g,i-n) the difference in C p , obtained on the different (a-g) poly-hexcore and (h-n) tetrahedral-based grids, with respect to the finest grid of each type.

Fig. 7 .
Fig. 7. C D A values obtained with varying y + mean value for rough cyclist surfaces.The percentage differences are with respect to the finest grid (0.125).

Fig. 8 .
Fig.8.Contours of (a,g) mean C p for the finest grid (0.125) and (b-f,h-l) the difference in C p on the grids with different y + mean with respect to the finest grid, obtained on the (a-f) poly-hexcore and (g-l) tetrahedral-based grids for rough cyclist surfaces.

Fig. 9 .
Fig. 9. C D A values obtained with varying y + mean value for smooth cyclist surfaces.The percentage differences are with respect to the finest grid (0.125).

Fig. 10 .
Fig. 10.Contours of (a,g) mean C p for the finest grid (0.125) and (b-f,h-l) the difference in C p , obtained on the different (a-f) poly-hexcore and (g-l) tetrahedralbased grids for smooth cyclist surfaces; with respect to the finest grid of each type.

Fig. 11 .
Fig. 11.Volume grid slice in centerplane z = 0 and surface grid for the poly-hexcore and tetrahedral grids with varying prism-layer growth rate.

Fig. 12 .
Fig. 12. C D A values obtained with varying prism layer growth rate for rough cyclist surfaces.The percentage differences are with respect to the finest grid (1.05).

Fig. 13 .
Fig. 13.C D A values obtained with varying prism layer growth rate for smooth cyclist surfaces.The percentage differences are with respect to the finest grid (1.05).

Fig. 14 .
Fig. 14.For the finest (a-c) poly-hexcore and (d-f) tetrahedral grids (1.05), contours of (a,b,d,e) mean C p for (a,d) rough and (b,e) smooth cyclist surfaces and (c,f) the difference in mean C p between the smooth (k S = 0 mm) and rough (k S = 0.1 mm) body surfaces.The dashed line indicates the location of the horizontal crosssection, shown in Fig. 15.

Fig. 15 .
Fig. 15.Horizontal cross-section at 0.95 m height through the upper left arm (location indicated in Fig. 14) showing 2D streamlines and contours of mean velocity ratio U r for the (a,c) rough (k S = 0.1 mm) and (b,d) smooth (k S = 0.0 mm) body surface for the finest (1.05) (a,b) poly-hexcore and (c,d) tet-based grid.The dashed lines indicate tangents on the inside and outside of the arm, which run parallel to the main flow direction.

Fig. 17 .
Fig. 17.Volume grid slice in centerplane z = 0 and surface grid for the poly-hexcore and tetrahedral grids.

T
. van Druenen and B. Blocken areas of lower C p on the frontal side and large areas of higher C p on the back side of the upper and lower legs.Moreover, substantial areas with dark red shades are visible on the side torso and the bicycle frame at the junction of the top and down tube.In the case of the SA model, dark red shades are also present on the outside of the upper arms.Additionally, dark blue shades are observed on the shoulders, face, inside of the upper arms and the bicycle frame near the head tube for both the RLZ k-ε and SA models.

Fig. 18 .
Fig. 18.C D A values obtained with varying maximum cell edge length in the wake zones.The percentage differences are with respect to the finest grid (1112).

Fig. 19 .
Fig. 19.Contours of (a,b) mean C p on the cyclist and bike surfaces and mean U r in the vertical centerplane z = 0 for the finest grid (1112) and (b-k,m-v) the difference in C p and U r with respect to the finest grid for the (a-k) poly-hexcore and (l-v) tetrahedral grids.

Fig. 21 .
Fig. 21.C D A values obtained by wind tunnel measurements and CFD simulations with varying turbulence modeling approaches/turbulence models for rough (SAS, T-SST, SST k-ω, SA) and smooth (RLZ k-ε) body surfaces.The percentage differences are with respect to the mean drag obtained from the wind tunnel measurements.

Fig. 22 .
Fig.22.Contours of (a,f) mean C p for the SAS approach and (b-e,g-j) the difference in C p , obtained for the different turbulence models for (a-c,e-h,j) rough and (d,i) smooth body surfaces on the (a-e) poly-hexcore and (f-j) tetrahedral-based grids; with respect to the SAS approach of each type.

•
It is not recommended to use steady RANS with the Realizable k-ε turbulence model and enhanced wall treatment, or with the Spalart-Allmaras turbulence model.Instead, the evaluation should consider Scale-Adaptive Simulation, or Steady RANS with the Transition SST turbulence model or the SST k-ω model for closure.

Table 1
CFD studies on cycling aerodynamics published in the peer-reviewed scientific literature.Overview of turbulence modeling, near-wall treatment and computational grid.

Table 2
Grid characteristics of evaluated grids with varying surface cell lengths.

Table 3
Grid characteristics of evaluated grids with varying y + mean value.
T.van Druenen and B. Blocken

Table 4
Mean area-weighted average y + values for rough cyclist surfaces.

Table 5
Grid characteristics of evaluated grids with varying prism-layer growth rate.