1. Introduction
Transport contributed 16 % of the world’s total greenhouse gas emissions in 2021 (ClimateWatch 2024). Despite extensive efforts to promote low-emission vehicles, there remains a critical need to improve the aerodynamic efficiency of road vehicles in order to reduce fuel consumption,
${\rm CO}_2$
emissions and energy consumption in electric vehicles, thereby extending their range. Aerodynamic drag is a major contributor to the resistance force encountered by vehicles when driving at speeds above approximately 70 km h−1 (Hucho Reference Hucho1987; Hucho & Sovran Reference Hucho and Sovran1993; Howell & Gaylard Reference Howell and Gaylard2006; Barnard Reference Barnard2009; Choi, Lee & Park Reference Choi, Lee and Park2014). The drag force is strongly influenced by the vehicle’s wake, which forms as the airflow separates from the rear of the vehicle. This separation leads to the formation of a low-energy, low-pressure wake region behind the body, which comprises the main component of the drag.
In research aimed at improving the vehicle aerodynamics, simplified vehicle models are often employed to reduce the complexity of simulations and focus on the fundamental aerodynamic phenomena. Such models, while stripped of many vehicle-specific details, are designed to preserve the key aerodynamic features that influence the wake dynamics. By increasing the `signal-to-noise’ ratio in these models, they allow researchers to focus on the essential flow phenomena that influence drag, without the interference of vehicle-specific characteristics.
Current understanding of the wakes behind simplified vehicle models identifies three periodic modes, the shear-layer dynamics, vortex shedding and bubble pumping, along with a stochastic mode, typically referred to as bi-modal switching (figure 1). These modes are typically characterised by specific Strouhal numbers: the shear-layer dynamics (
${\sim} 3$
), vortex shedding (
${\sim} 0.2$
), bubble pumping (
${\sim} 0.08$
) and bi-modal switching (
${\sim} 0.006$
) (Grandemange et al. Reference Grandemange, Gohlke and Cadot2013b
), where the Strouhal number for bi-modal switching represents an average `frequency’ as it is a stochastic process.
-
(i) Shear-layer dynamics stems from Kelvin–Helmholtz instabilities at the rear edges of the vehicle, creating distinct shear layers.
-
(ii) Vortex shedding occurs due to interactions between opposing shear layers, generating vertical and spanwise vortices in the wake.
-
(iii) Bubble pumping remains less well understood, with hypotheses ranging from periodically shed vortices at the end of the separation bubble to nonlinear interactions of vortex shedding modes (Berger, Scholz & Schumm Reference Berger, Scholz and Schumm1990; Duell & George Reference Duell and George1999; Khalighi et al. Reference Khalighi, Zhang, Koromilas, Balkanyi, Bernal, Iaccarino and Moin2001; Krajnović & Davidson Reference Krajnović and Davidson2003; Volpe, Devinant & Kourta Reference Volpe, Devinant and Kourta2015; Pavia, Passmore & Sardu Reference Pavia, Passmore and Sardu2017).
-
(iv) Bi-modality is characterised by symmetry breaking, which often manifests as spontaneous switching between two distinct wake states, driven by reflectional symmetry breaking in the wake. The underlying mechanism of the laminar symmetry breaking can be explained through nonlinear analysis of the Navier–Stokes equations (Meliga, Chomaz & Sipp Reference Meliga, Chomaz and Sipp2009; Zampogna & Boujo Reference Zampogna and Boujo2023). The symmetry breaking persists into the turbulent region. However, in this case, turbulent perturbations induce the wake to switch randomly between asymmetric positions (Rigas et al. Reference Rigas, Morgans, Brackston and Morrison2015).
Figure 1. A schematic illustrating different dynamic modes observed in a bluff-body wake: (a) shear-layer dynamics, (b) vortex shedding, (c) bi-modal switching and (d) bubble pumping.
Because these wake dynamics is inherently three-dimensional and complex, wake modifications are often best understood through their impact on the global flow topology rather than through localised changes in velocity or pressure fields. A commonly used framework for describing such modifications is to refer to the global mean topology of each bi-stable state, which is characterised by a toroidal vortex structure that reorients with each bi-modal switch (Lucas et al. Reference Lucas, Cadot, Herbert, Parpais and Délery2017; Dalla Longa, Evstafyeva & Morgans Reference Dalla Longa, Evstafyeva and Morgans2019). This toroidal representation simplifies the characterisation of large-scale wake changes, providing an intuitive means of analysing wake modifications and their aerodynamic consequences. Throughout this paper, wake shape, size and bi-modal properties will be examined primarily through this lens.
Most vehicle wake dynamics studies have employed either the Ahmed body (Ahmed, Ramm & Faltin Reference Ahmed, Ramm and Faltin1984) or the Windsor body (Windsor Reference Windsor1991), which represent simplified square-back vehicles in ground proximity without wheels. However, it has been shown that, while the addition of wheels does not substantially change the dominant dynamic structures in the wake, it does modify the wake’s shape (Pavia & Passmore Reference Pavia, Passmore and Wiedemann2018; Su et al. Reference Su, He, Xu, Gao and Krajnović2023), bringing it closer to the behaviour observed in real vehicles. In this study, we use the Windsor body with wheels as the baseline geometry. This retains the simplified nature of canonical bluff-body configurations, maintaining a focus on fundamental wake dynamics, while introducing enough realism to enhance the relevance and generalisability of the findings to full-scale automotive applications.
Drag on bluff bodies is primarily due to pressure drag, making wake modification a key focus for drag-reduction strategies. This study examines rear roof edge modifications to assess their impact on drag. The design of vehicles to reduce drag is a complex challenge, especially when striving to balance functional design with aerodynamic efficiency. Roof extensions are a promising option in this context because they offer a relatively simple modification that does not significantly alter the overall vehicle geometry. While such modifications have been employed in the automotive industry for decades, fundamental studies guiding their design remain limited. Some investigations have explored roof extensions on industrial geometries, such as Volvo cars (Sterken et al. Reference Sterken, Lofdahl, Sebben and Walker2014), whereas other studies have examined passive-control techniques like tapers (Littlewood & Passmore Reference Littlewood and Passmore2010; Perry, Pavia & Passmore Reference Perry, Pavia and Passmore2016; Bonnavion & Cadot Reference Bonnavion and Cadot2019), flaps (Grandemange et al. Reference Grandemange, Mary, Gohlke and Cadot2013c ; Garcia de la Cruz, Brackston & Morrison Reference Garcia de la Cruz, Brackston and Morrison2017; Urquhart et al. Reference Urquhart, Varney, Sebben and Passmore2020b ; Fan et al. Reference Fan, Parezanović, Fichera and Cadot2024) and cavities (Evrard et al. Reference Evrard, Cadot, Herbert, Ricot, Vigneron and Délery2016; Lucas et al. Reference Lucas, Cadot, Herbert, Parpais and Délery2017; Bonnavion & Cadot Reference Bonnavion and Cadot2018; Lorite-Díez et al. Reference Lorite-Díez, Jiménez-González, Pastur, Cadot and Martínez-Bazán2020; Muñoz-Hervás et al. Reference Muñoz-Hervás, Lorite-Díez, García-Baena and Jiménez-González2024), with some aiming to elucidate the underlying flow mechanisms. Although the physics governing roof extension drag reduction shares similarities with these devices, the exact mechanisms are likely to be a combination of various effects unique to the configuration studied. Nevertheless, it is useful to review the principal findings of prior passive-control studies, as they provide the physical context in which the present optimisation framework should be understood.
For constant-length flaps and tapers, several experimental studies have examined the sensitivity of drag to the angle of incidence, often using two-dimensional Particle Image Velocimetry to infer the underlying flow mechanisms. A particularly influential example is the work of Grandemange et al. (Reference Grandemange, Mary, Gohlke and Cadot2013c ), who investigated the aerodynamic response of top and bottom flaps, referred to here as extensions, on a square-back body. Their analysis decomposed drag variations into distinct contributions: a linear term associated with base-size reduction and quadratic terms linked to lift-induced drag. The latter was further separated into a component due to overall body lift and another attributed to the tip vortices formed along the flap edges. While this framework provided valuable phenomenological insight, it was inherently constrained by the lack of three-dimensional flow-field data, leading to an overemphasis on inviscid lift-induced mechanisms and limited consideration of viscous and wake-interaction effects.
Subsequent studies developed broader interpretations of the same phenomena. The effects that Grandemange et al. (Reference Grandemange, Mary, Gohlke and Cadot2013c ) described as lift-induced drag were later reconsidered in terms of vertical `wake balance’, that is, whether the time-averaged wake is dominated by upwash or downwash, without directly invoking induced-drag analogies. It is important to emphasise that wake balance is an inherently broad concept that has been interpreted in several distinct ways. It can be quantified via the vertical base-pressure gradient at the centre of the base, reflecting the global circulation and the relative strength of the upper and lower recirculation regions. In the present work, we adopt a definition based on the inclination of the central wake streamlines relative to the streamwise direction. One can, for instance, envisage a wake in which the upper recirculation region is larger than the lower one, yet the mean streamlines impinge on the base nearly perpendicularly. In viscous configurations, particularly those including ground clearance and wheels, which modify the momentum fluxes above and below the body, the purely inviscid, lift-induced interpretation becomes less meaningful. In such cases, a definition based on streamline orientation or momentum balance provides a more robust and physically interpretable framework, as it directly captures the effective inflow conditions experienced by the wake and offers a more functional basis for optimisation.
Within this framework, Urquhart et al. (Reference Urquhart, Varney, Sebben and Passmore2020b ) examined flap optimisation using wake balance as an integral descriptor of the global flow state, demonstrating that approaches toward vertical symmetry correlate strongly with reduced drag. Haffner et al. (Reference Haffner, Castelain, Borée and Spohn2021) further isolated the influence of symmetry by introducing controlled spanwise perturbations, showing that vertically balanced wakes consistently exhibit lower drag, even when bi-modal switching persists. Complementary investigations of rear tapers have examined how geometric shaping affects symmetry and switching behaviour (Perry et al. Reference Perry, Pavia and Passmore2016; Bonnavion & Cadot Reference Bonnavion and Cadot2019). In particular, Bonnavion & Cadot (Reference Bonnavion and Cadot2019) demonstrated that the base aspect ratio and pressure-gradient distribution govern both the frequency and persistence of bi-modal transitions.
Furthermore, several studies have examined the performance of flaps and tapers under varied flow conditions, such as yaw (Garcia de la Cruz et al. Reference Garcia de la Cruz, Brackston and Morrison2017; Fan et al. Reference Fan, Parezanović, Fichera and Cadot2024), or explored modified configurations, including flexible flaps (Muñoz-Hervás et al. Reference Muñoz-Hervás, Lorite-Díez, García-Baena and Jiménez-González2024). These works provide valuable insights from a practical implementation perspective; however, the fundamental mechanisms of drag reduction remain consistent with those identified in the aforementioned studies, and are therefore not discussed in detail here.
Taken together, these studies demonstrate that devices acting on the rear edges, such as flaps and tapers, primarily influence drag through base-size reduction and wake reorganisation mechanisms linked to vertical balance and bi-modal switching. A complementary line of research, focusing on cavity-based modifications, was presented by Evrard et al. (Reference Evrard, Cadot, Herbert, Ricot, Vigneron and Délery2016), who showed that incorporating a cavity can suppress bi-modality by stabilising the shear-layer dynamics and revealed the continuous nature of this transition.
Collectively, these studies provide a strong experimental foundation for understanding passive flow control on square-back geometries. Nonetheless, most of these studies rely on two-dimensional or surface-based diagnostics and vary a single geometric parameter at a time, limiting their ability to capture the coupled three-dimensional interactions among wake features, particularly the interplay between taper-induced and cavity-like effects. This limitation highlights the need for an integrated perspective, one that captures how multiple geometric parameters interact within a fully three-dimensional wake. The present study addresses this by optimising the geometric parameters governing both length and angle concurrently and analysing the resulting three-dimensional flow fields in detail.
Optimising drag-reduction devices involves multiple design parameters, and predicting their effects is challenging. Traditional methods such as Computational Fluid Dynamics (CFD) simulations or wind-tunnel testing are time consuming and costly, especially when dealing with large design spaces. Classical optimisation algorithms often struggle with highly nonlinear, multi-peak problems. Evolutionary algorithms are effective but require many function evaluations, making them impractical for CFD-based design processes.
Surrogate models, also known as response surfaces, have been proposed to address these challenges. These models can accurately describe nonlinear variations and be used to guide optimisation. Surrogate model optimisation, in particular, is efficient and provides valuable insights into input dependencies. However, traditional deterministic models, such as polynomial surrogate models, may fail to capture the full design space, often leading to local optima. To overcome this, Bayesian optimisation has been introduced, leveraging probabilistic models to balance exploration and exploitation. Notably, Jones, Schonlau & Welch (Reference Jones, Schonlau and Welch1998) proposed the efficient global optimisation (EGO) algorithm, which integrates Gaussian processes (GPs) with the concept of expected improvement to enhance search efficiency.
The EGO algorithm is particularly useful for optimisation in complex design spaces. By balancing exploration and exploitation, it can efficiently locate optimal designs. Recent applications of this approach, such as in Yan et al. (Reference Yan, Lai, Wang, Hu and Deng2019), have demonstrated its effectiveness in optimising the vehicle aerodynamics, achieving significant drag reductions with fewer design evaluations. This study was based on Reynolds-averaged Navier–Stokes (RANS) simulations and did not seek to link the observed drag-reduction trends to the underlying aerodynamics. In contrast, we employ high-fidelity, wall-resolved large eddy simulations (WRLES), which capture more complex and complete wake behaviours. Notably, Dalla Longa et al. (Reference Dalla Longa, Evstafyeva and Morgans2019) showed that WRLES can effectively capture bi-modal switching, along with other key features of bluff-body flows. Flow fields corresponding to various drag coefficient regimes are analysed, and alterations in the flow field are linked to the observed trends.
It is worth noting that Bayesian optimisation is not the only approach to incorporate uncertainty into the optimisation process. For example, Urquhart et al. (Reference Urquhart, Varney, Sebben and Passmore2020b ) implemented a strategy developed by Goel et al. (Reference Goel, Haftka, Shyy and Queipo2007), which utilised an ensemble of 10 surrogate models based on radial basis function surrogate models. Their method estimated uncertainty by combining the median for prediction and the standard deviation for variability. The algorithm alternated between evaluating points in the design space with the lowest median and those with the highest standard deviation. While effective, these approaches lack the mathematical formalism and adaptive exploration–exploitation balance offered by classical Bayesian optimisation. For this reason, we adopt Bayesian optimisation in this study to efficiently explore the design space while maintaining robustness in optimisation outcomes.
The present study presents a detailed analysis of the aerodynamic effects of roof extensions for road vehicles, focusing on the optimisation of key parameters: taper penetration distance, angle of incidence and extension length. The primary goal of the optimisation study is to minimise the drag coefficient. To provide a comprehensive aerodynamic analysis, a surrogate model is also developed for the lift coefficient. Furthermore, the optimisation function is extended to simultaneously account for both drag and lift, presenting geometry configurations corresponding to minimal drag at specific lift coefficient requirements. The study goes beyond reporting individual mechanisms by showing how they interact when multiple geometric parameters are varied simultaneously, yielding a new system-level understanding and practical design rules. In doing so, this work connects fundamental wake studies with a more practically oriented design framework, and demonstrates how system-level interactions can be exploited to achieve improved aerodynamic performance. These findings establish a robust framework for aerodynamic optimisation and reinforce the effectiveness of Bayesian optimisation in CFD-based design. The remainder of the paper is structured as follows: § 2 details the numerical set-up and data analysis tools, § 3 presents the optimisation results and surrogate models that illustrate the dependence of drag and lift on each parameter, along with partial dependence studies that reveal the physical mechanisms governing drag variations. Section 4 provides the conclusions of the study.
2. Methodology
2.1. Numerical simulation set-up
2.1.1. Geometry and flow domain
The primary geometry investigated in this study is the Windsor body, a modified square-back Ahmed body originally developed by Rover Group (a former UK vehicle manufacturer) for fundamental aerodynamic research (Windsor Reference Windsor1991). Pavia & Passmore (Reference Pavia, Passmore and Wiedemann2018) introduced further modifications by incorporating four stationary wheels, which were simplified in this study to cylindrical shapes with closed inner geometry, in contrast to the detailed design used by Page & Walle (Reference Page and Walle2022). Additionally, mounting pins were removed due to their redundancy in simulation studies.
To examine passive-control strategies, the rear roof edge was modified with a parameterised extension characterised by three parameters: taper penetration distance
$A$
, angle of incidence
$\beta$
and extension length
$C$
. The baseline geometry, along with the parameterised extension, is illustrated in figure 2. The parameter ranges are as follows:
$A/L$
– taper penetration distance,
$0{-}0.07$
;
$\beta$
– angle of incidence,
$-20^\circ$
to
$20^\circ$
; and
$C/L$
– extension length,
$0{-}0.07$
, where
$L$
represents vehicle’s length. The upper bounds of
$A/L$
and
$C/L$
were set to reflect practical design constraints: although larger values can yield further drag reductions, they correspond to extensions that are unrealistic for vehicle applications. A limited number of cases with
$A/L\gt 0.07$
were nevertheless explored outside the optimisation domain, solely to aid interpretation of flow-physics trends. For the angle of incidence: the range was chosen based on previous studies indicating that most drag reduction occurs between
$5^\circ$
and
$15^\circ$
(Littlewood & Passmore Reference Littlewood and Passmore2010; Grandemange et al. Reference Grandemange, Mary, Gohlke and Cadot2013c
; Pavia, Passmore & Gaylard Reference Pavia, Passmore and Gaylard2016; Perry et al. Reference Perry, Pavia and Passmore2016; Bonnavion & Cadot Reference Bonnavion and Cadot2019; Urquhart et al. Reference Urquhart, Varney, Sebben and Passmore2020b
; Fan et al. Reference Fan, Parezanović, Fichera and Cadot2024), while negative angles (upsweep) were retained to capture the full lift response and to test the possibility that a modest upsweep might improve vertical wake balance. Their inclusion had negligible impact on efficiency, since Bayesian optimisation naturally allocates very few samples to unfavourable regions of the design space. For reference, the case
$C/L=0$
corresponds to a pure taper without an added surface, conceptually similar to the slant geometries studied by Littlewood & Passmore (Reference Littlewood and Passmore2010).
Figure 2. Visualisation of the baseline geometry and the parameterised roof extension: (a) front view, (b) side view and (c) zoomed-in view of the roof extension. A zero extension length (
$C/L$
) corresponds to a pure taper configuration (no additional extension surface).
The Reynolds number based on the body height was set at
$Re_H = 3.3 \times 10^4$
, consistent with previous experimentally validated simulation studies on simplified road vehicles (Hesse & Morgans Reference Hesse and Morgans2021; Ahmed & Morgans Reference Ahmed and Morgans2023). Although Windsor-type bodies exhibit some Reynolds-number sensitivity up to
$Re_H \approx 3\times 10^5$
(Varney Reference Varney2020), this dependence is mainly associated with front-body phenomena such as frontal separation bubbles. In contrast, the present study focuses on the wake dynamics, for which the boundary layer is already turbulent at separation and the large-scale flow behaviour remains effectively Reynolds-independent within this range.
This was confirmed by an additional WRLES of the baseline configuration at a higher Reynolds number,
$Re_H = 0.96\times 10^5$
, matching the set-up of Su et al. (Reference Su, He, Xu, Gao and Krajnović2023). The higher-Re case exhibited the expected small decrease in drag (
$5\,\%$
) and an almost unchanged wake topology, with the non-dimensional vertical base-pressure gradient varying only from
$-0.016$
to
$-0.019$
. The higher-Re case is also used later in the paper to support validation of the simulation framework.
The aspect ratio of the body cross-section was set at
$W/H = 1.35$
, where
$W$
and
$H$
represent the body width and height, respectively. The body length was
$L = 3.61H$
, and the distance between the underfloor of the body and the ground was maintained at
$C = 0.174H$
. Additionally, to maintain consistency with prior studies, a gap of
$0.01H$
was kept between the wheels and the floor.
The computational domain dimensions were configured according to ERCOFTAC guidelines, specified as
$(L_x, L_y, L_z) = (8L, 2L, 2L)$
, where
$L_x, L_y, L_z$
denote the domain size in the streamwise, transverse and spanwise directions, respectively. This set-up results in a blockage factor of
$2.5\,\%$
. The distance between the front of the body and the domain inlet was set at
$2L$
. A uniform velocity profile was prescribed at the inlet, and no-slip conditions were applied to both the geometry surfaces and the domain floor. The floor velocity matched the free-stream velocity to represent moving ground conditions, while the wheels remained stationary. A convective boundary condition was applied at the outlet, and slip conditions were used on all other boundaries.
2.1.2. Numerical solver
The incompressible Navier–Stokes equations were solved using OpenFOAM version 10 with the finite volume method. The PIMPLE – a hybrid CFD algorithm combining the PISO (Pressure-Implicit with Splitting of Operators) and SIMPLE (Semi-Implicit Method for Pressure Linked Equations), was employed as the solver.
Sub-grid-scale turbulence was modelled using the wall-adapting local eddy-viscosity model. The mesh was generated using StarCCM+ with the trimmed cell mesher, ensuring high-quality grid resolution. Mesh refinement studies, presented in subsequent sections, were performed to ensure numerical accuracy and solution independence.
Convergence was monitored at both numerical and statistical levels. Within each PIMPLE time step, the momentum and pressure equations were solved until residuals decreased by at least three orders of magnitude. The pressure field was solved using the Preconditioned Conjugate Gradient algorithm with a diagonal incomplete-Cholesky preconditioner (
$\text{tolerance}=10^{-8}$
,
$\text{relTol}=0.01$
), and the velocity field using the Preconditioned Bi-Conjugate Gradient solver with a Diagonal-based Incomplete LU preconditioner (
$\text{tolerance}=10^{-8}$
,
$\text{relTol}=0$
). Final corrector steps were tightened to
$10^{-9}$
, resulting in residuals typically between
$10^{-6}$
and
$10^{-8}$
, which remained stable throughout the simulations.
Each WRLES case was initialised from a converged RANS solution and advanced for approximately
$25\,T_c$
(
$T_c = H/U_\infty$
) to remove initial transients. Statistical averaging was then performed over an additional
$150\,T_c$
, during which the running-average drag and lift coefficients were monitored. The averaging window was progressively extended in increments of
$25\,T_c$
until these mean quantities changed by less than 0.25 % over two successive intervals, confirming statistical convergence. A minimum sampling duration of
$150\,T_c$
was imposed, based on the convergence behaviour of the baseline case, to avoid false stabilisation of short averaging windows and ensure that all configurations captured the dominant wake dynamics. For configurations exhibiting bi-stable wake behaviour, conditional averaging was applied so that statistics were computed within a single stable wake state. In contrast, for cases near the bifurcation point or in nearly symmetric conditions, where the two asymmetric states were dynamically similar and switching occurred on short time scales, global averaging was adopted, as the associated drag fluctuations were small and conditional averaging produced negligible differences. This averaging strategy was designed to capture statistically representative flow states while maintaining computational feasibility across the large set of optimisation cases. Because the characteristic switching frequency of wake transitions is typically very low (
$St_w \lesssim 10^{-2}$
), the flow spends the majority of the time in one of the stable configurations. Consequently, conditional and global averages yield nearly identical mean drag estimates, and the adopted approach provides a physically meaningful representation of the dominant flow dynamics.
2.1.3. Mesh refinement study
The mesh refinement study was initially conducted using the Windsor body without wheels and with a stationary ground to facilitate consistent validation against previous studies. This choice was deliberate, as it enabled direct comparison with prior studies on the same geometry (Varney et al. Reference Varney, Pavia, Passmore and Crickmore2021; Su et al. Reference Su, He, Xu, Gao and Krajnović2023), which would not have been possible for the wheeled case due to geometric differences. Moreover, given the Windsor body’s close similarity to the Ahmed body, it allowed semi-direct comparison with several WRLES studies performed within the same research group (Hesse & Morgans Reference Hesse and Morgans2021; Ahmed & Morgans Reference Ahmed and Morgans2023). These comparisons informed the choice of domain extents, near-wall resolution, and wake-region refinement in independently validated numerical set-ups. Once the baseline mesh resolution was established on the wheel-free configuration, the wheels and moving floor were introduced, and a second refinement study was conducted. Owing to the geometric differences and lack of direct reference data for the wheeled case, this second stage focused exclusively on verifying numerical mesh independence rather than experimental validation.
The refinement targeted both the wake and near-wall regions. Three different meshes were generated, all employing the same refinement regions and growth ratios. Thus, varying the size of the cells in the most refined region, i.e. the near-wake region (
$\Delta L_{wake}$
), resulted in changes in the size of all cells in the domain, along with the number of prism layers (to maintain constant growth ratios). The first-layer height at the wall was kept constant across all meshes to maintain low
$y^+$
values; the domain-averaged value was
$\overline {y^+}\approx 0.58$
(refer to figure 5 for distribution), with majority of the surface having
$y^+\lt 1$
. Additionally, the average normalised near-wall grid size in the streamwise and spanwise directions,
$\Delta x_{1,3}^+ = {\Delta x_{wake} \boldsymbol{\cdot }u_{\tau }}/{\nu }$
, was
${\sim} 30$
. The prism-layer normal growth ratio lay in the range
$1.3$
–
$1.5$
depending on the mesh configuration, and the transition ratio from the last prism layer to the adjacent outer cell was likewise maintained between
$1.3$
and
$1.5$
to avoid abrupt stretching. Table 1 presents the configurations for the three meshes, showing the final cell count, near-wake cell size and the number of prism layers. The near-wake cell size is expressed in terms of the Taylor microscale, approximated as
$\lambda _T \sim 5.5 {Re_H}^{-1/2}H$
(Hesse & Morgans Reference Hesse and Morgans2023), to provide a reference for the extent of resolved and modelled flow fields. The mesh sizes were chosen based on the studies by Hesse & Morgans (Reference Hesse and Morgans2023) and Ahmed (Reference Ahmed2023), as well as the computational resources available.
Table 1. Mesh refinement study for the Windsor body without wheels: conditionally time-averaged drag and lift coefficients (
$C_D$
and
$C_L$
), base-pressure coefficient (
$C_{\textit{PB}}$
) and recirculation region length (
$L_B$
).
The analysis of key conditionally time-averaged simulation results, i.e. time-averaged values for one of the bi-modal wake positions, such as the drag and lift coefficients, shows a high degree of convergence between the intermediate and highly resolved meshes, with errors in all coefficients remaining below 1 %.
Additionally, Q-criterion contours and turbulent eddy-viscosity fields were examined. The Q-criterion contours (figure 3 a, b, c) demonstrate that all three meshes successfully resolve the primary coherent flow structures, including the hairpin vortices shed by the front separation bubbles, and A-pillar vortices. Notably, a higher degree of similarity is observed between the medium and fine meshes. This is evident in the structures at the initiation of the A-pillar vortex, as well as the amount of vortical structures emitted by the side separation bubbles and the shear layers in the wake.
Figure 3. Isosurfaces of
$ Q^*$
with
$ Q^* = 5.3$
(top row) and viscosity ratio
$ \nu _T / \nu = 5$
(bottom row) for the Windsor body without wheels: (a, d) coarse mesh, (b, e) medium mesh and (c, f) fine mesh. The colour bar indicates the magnitude of the conditionally time-averaged streamwise velocity.
Figure 3(d, e, f) presents isosurfaces of the viscosity ratio (
$\nu _t/\nu = 5$
) to highlight regions where the modelled turbulent effects exceed recommended values. It is important to note that the large downstream structures result from artificial grid effects; specifically, large cell sizes amplify even weak strain effects, as the eddy viscosity depends on grid-dependent length scales. These artificial effects can be disregarded in the analysis. Overall, all three meshes perform well in resolving the main flow structures. However, the coarsest mesh exhibits the weakest performance, particularly in the near-wake region. Based on this analysis, the medium-size mesh (approximately 11 million cells) was selected as the base mesh, consistent with resolutions shown to be sufficient for accurate WRLES of similar square-back configurations using unstructured meshes (Hesse & Morgans Reference Hesse and Morgans2021; Ahmed & Morgans Reference Ahmed and Morgans2023).
It is noted that small Q-criterion iso-surfaces appear near the lower side edges of the body (see figure 3). These structures develop within a region of uniform mesh and are therefore not numerical artefacts. They originate from the confined underbody flow: as the velocity beneath the body increases due to the Venturi-like acceleration between the body and the ground, the resulting pressure difference induces an inward curling of the outer flow around the lower side edge. This process generates a coherent vortex along the bottom edge of the body, which interacts with the floor boundary layer and produces secondary counter-rotating vortices. These near-wall interactions are a well-known feature of confined ground-effect flows and are physically consistent with the observed vortex patterns.
Confidence in the wheel-free simulation set-up was further supported by a targeted validation against the simulation-based study of Su et al. (Reference Su, He, Xu, Gao and Krajnović2023). In this comparison, the Reynolds number was adjusted to match their conditions (
$Re_H = 0.96 \times 10^5$
). To account for the Reynolds-number increase, the mesh size was increased to approximately 19 million cells to achieve comparable Taylor-microscale resolution. Under these matched conditions, our simulation yielded a drag coefficient of
$C_D = 0.307$
, differing by only
$1\,\%$
from their reported value of
$C_D = 0.303$
, thereby providing strong evidence for the consistency and reliability of the present numerical configuration. In addition to the close agreement in
$C_D$
, the qualitative wake behaviour, including the vertical wake balance and streamline topology, closely matched that reported by Su et al. (Reference Su, He, Xu, Gao and Krajnović2023), further confirming that the large-scale wake dynamics is faithfully reproduced.
The lowest-Reynolds-number experimental study on the Windsor body identified in the literature was conducted at
$Re_H \approx 2 \times 10^5$
(Varney Reference Varney2020), and reported a corrected drag coefficient of
$C_D^{\textit{corr}} = 0.317$
, where the correction accounts for wind-tunnel blockage effects according to
\begin{equation} C_D^{\textit{corr}} = C_D \left (\frac {A_{\textit{WT}}}{A_{\textit{WT}} - A_{{\kern-1pt}f}}\right )^{-1}, \end{equation}
with
$A_{\textit{WT}}$
representing the cross-sectional area of the wind tunnel and
$A_f$
the frontal area of the body. This value lies within
$3\,\%$
of that obtained in our simulation, which considering the modest Reynolds-number sensitivity discussed in § 2.1.1, constitutes close agreement and further supports the validity of the numerical set-up.
Following the validation and mesh-convergence analysis for the wheel-free configuration, an analogous refinement study was conducted for the complete model including the wheels and moving floor. The same refinement methodology was applied as described above, using the established base wake cell size
$\Delta x_{wake}$
and refinement ratios. The wheel and wake regions shared the same refinement level to ensure consistent resolution of the wheel–wake interaction. Although resolving the smaller wheel geometry would ideally require a locally finer mesh, this was not adopted here to prevent the wheel region from becoming a computational bottleneck due to Courant-number restrictions. Since the present study primarily targets the wake dynamics rather than detailed wheel flow features, this compromise was considered appropriate.
The total cell count in the wheeled configuration ranged from 10 to 20 million across the coarse, medium and fine meshes. The conditionally time-averaged aerodynamic coefficients for these three meshes were as follows: (
$C_D$
,
$C_L$
,
$C_{\textit{PB}}$
,
$L_B/H$
) = (0.396,
$-0.067$
,
$-0.224$
, 1.48), (0.402,
$-0.037$
,
$-0.227$
, 1.46) and (0.404,
$-0.036$
,
$-0.227$
, 1.46), respectively. Variations between the two finest meshes remained within
$3\,\%$
for
$C_D$
,
$C_L$
and
$C_{\textit{PB}}$
, confirming mesh independence at this level of resolution. The medium mesh (approximately 14 million cells) was therefore selected for all subsequent simulations as a balance between numerical accuracy and computational cost.
The qualitative flow structures in the wheeled configuration closely resembled those of the wheel-free case, with the main distinction being the additional vortical structures shed by the wheels themselves. Representative
$Q^*$
iso-surfaces for the medium mesh are shown in figure 4, illustrating the same dominant vortical features as in the wheel-free configuration. For brevity, viscosity-ratio isosurfaces are not shown, as they exhibit an essentially identical distribution. Notably, the small vortical structures previously observed near the lower side edges of the body are no longer present, as the introduction of a moving ground alters the boundary-layer development and thus the structures formed due to the induced circulation. Because the wheels were modelled as simplified cylinders, their wakes remain relatively compact and primarily affect the underbody momentum balance, without significantly modifying the overall wake dynamics.
Figure 4. Isosurfaces of
$ Q^*$
with
$ Q^* = 5.3$
for the Windsor body with wheels: (a) coarse mesh, (b) medium mesh and (c) fine mesh. The colour bar indicates the magnitude of the conditionally time-averaged streamwise velocity.
As noted earlier, refinement of the wheel region represents a computational bottleneck. Nonetheless, an additional targeted refinement test was carried out to verify that the chosen resolution was sufficient. The wheel-region cell size was reduced to
$\Delta x_{{wheels}} = 0.175\lambda _T$
, yielding a drag coefficient of
$C_D = 0.395$
, compared with
$C_D = 0.402$
for
$\Delta x_{{wheels}} = 0.35\lambda _T$
. The finer refinement therefore altered the drag by only approximately
$2\,\%$
, while approximately doubling the computational cost due to Courant-number constraints. Given this small difference and the large number of simulations required for the optimisation campaign, the coarser refinement level (
$\Delta x_{{wheels}} = 0.35\lambda _T$
) was retained for all subsequent analyses.
Lastly, considering all geometric and set-up differences, the results were compared with other studies. The baseline simulation with the present mesh yielded a corrected drag coefficient of
$C_D = 0.382$
(uncorrected:
$0.402$
) at
$Re_H = 3.3 \times 10^4$
. For reference, Varney et al. (Reference Varney, Pavia, Passmore and Crickmore2021) reported a corrected value of
$C_D = 0.372$
at
$Re_H \approx 7.7 \times 10^5$
in experiments, while Su et al. (Reference Su, He, Xu, Gao and Krajnović2023) obtained an uncorrected
$C_D = 0.445$
at
$Re_H = 0.96 \times 10^5$
in simulations. Both reference configurations were tested at higher Reynolds numbers and employed less streamlined wheel geometries than the present model. These differences make direct quantitative comparison difficult, as drag in this Reynolds-number range is known to be sensitive to both wheel geometry and Reynolds-number effects (Varney et al. Reference Varney, Pavia, Passmore and Crickmore2021). Nevertheless, the present results fall between the experimental and numerical values, which is physically consistent given the differing configurations. This level of agreement, together with the validation of the wheel-free set-up and the refinement study for the wheeled case, supports confidence in the overall numerical framework. We acknowledge that minor discrepancies may persist in the fine-scale representation of wheel-shed structures; however, since the focus of this work is on the rear-wake dynamics and global drag trends, this simplification represents a deliberate and appropriate balance between physical fidelity and computational cost. The final mesh, comprising approximately 15 million cells and showing all refinement regions with close ups of key zones, is presented in figure 5.
Figure 5. Visualisation of the intermediate resolution mesh with imposed
$ y^+$
values on the body surface: (a) wheel close up, (b) prism-layer close up, (c) side view and (d) top view. The colour bar on the right indicates the magnitude of
$ y^+$
values.
2.2. Data analysis
This section outlines the mathematical equations for key variables in flow-field analysis, along with explanations of the techniques employed for their calculation in non-trivial scenarios. Additionally, it covers a brief description of the optimisation algorithms used in this study.
2.2.1. Key variables
Flow-field analysis primarily revolves around velocity and pressure plots. Pressure can be non-dimensionalised to form the pressure coefficient,
$C_p$
, defined as
\begin{equation} C_p = \frac {p - p_\infty }{\dfrac {1}{2} \rho U^2} ,\end{equation}
where
$p$
represents the pressure at a specific point,
$p_\infty$
is the free-stream pressure,
$\rho$
is the fluid density and
$U$
denotes the free-stream velocity.
Investigation into the forces acting on the body entails the utilisation of force coefficients, defined as
\begin{equation} C_X = \frac {F_X}{\dfrac {1}{2} \rho U^2 A} ,\end{equation}
where
$F_X$
denotes the analysed force component, and
$ A$
is the projected frontal area. Specifically:
$ X = D$
for the drag coefficient
$ C_D$
, corresponding to the force aligned with the body axis, which in the case of zero yaw coincides with the streamwise direction,
$ X = L$
for the lift coefficient
$ C_L$
, representing the force in the vertical direction,
$ X = S$
for the side force coefficient
$ C_S$
, referring to the lateral force perpendicular to the body’s sides.
Additionally, the non-dimensional base-pressure gradient was computed to quantify the vertical and spanwise wake balance
\begin{equation} \frac {\partial C_p(t)}{\partial n} = \frac {1}{M} \sum _{m=1}^{M} \frac {\sum _M C_P \left(\dfrac {\Delta n}{2}, t\right) - \sum _M C_P \left(-\dfrac {\Delta n}{2} t\right)}{ \Delta n/L} ,\end{equation}
where
$n$
denotes either the vertical (
$y$
) or spanwise (
$z$
) direction,
$\Delta n$
is the probe spacing in the direction of the gradient calculation (equal to
$0.22W$
for the spanwise locations and
$0.7H$
for the vertical ones),
$L$
is the characteristic length scale (the body height
$H$
or width
$W$
) and
$M$
is the number of probe pairs used for averaging along the base surface.
To identify regions of high rotation rate, the
$Q$
-criterion and its non-dimensional version are employed
\begin{align} Q &= \frac {\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{U}}{\|\boldsymbol{\nabla }\boldsymbol{U}\|} ,\end{align}
\begin{align} Q^* &= \frac {Q H^2}{U^2} ,\end{align}
where
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{U}$
represents the divergence of velocity, and
$\|\boldsymbol{\nabla }\boldsymbol{U}\|$
is the magnitude of the velocity gradient.
In studies of steady geometrical modifications, the analysis of time-averaged flow fields is generally the most effective means of identifying the mechanisms responsible for drag changes. The mean pressure and velocity distributions directly reflect variations in separation topology, wake balance and shear-layer organisation, which are the primary drivers of aerodynamic performance. For this reason, the present discussion emphasises statistically converged mean quantities, either conditionally or globally averaged as described in § 2.1.2. Nevertheless, all simulations are fully time resolved, and unsteady flow information was examined whenever questions arose that could not be addressed through mean fields alone.
2.3. Optimisation
The optimisation procedure applied to the three roof extension variables,
$A$
,
$\beta$
and
$C$
, employs the EGO algorithm (Jones et al. Reference Jones, Schonlau and Welch1998), as implemented by Head et al. (Reference Head, Kumar, Nahrstaedt, Louppe and Shcherbatyi2021). This approach, rooted in Bayesian optimisation, utilises GPs to construct a surrogate model that predicts how the drag coefficient depends on
$A$
,
$\beta$
and
$C$
, interpolating between previously evaluated points. The optimisation process is guided by an expected improvement (EI) function, which effectively balances exploration and exploitation to minimise the number of function evaluations required. This methodology is particularly well suited for expensive-to-evaluate functions, where the computational cost of each evaluation is high. A schematic illustrating the optimisation workflow is provided in figure 6.
Figure 6. (a) Schematic illustrating the Bayesian optimisation loop. (b) Expected improvement schematic. Red dots represent collected data points, the black line denotes the GP model predictions, and the blue shaded regions indicate the variance predicted by the model. The purple distribution visualises the probability distribution of a sample at
$ X$
achieving value
$ Y$
.
A GP with mean
$\mu$
and variance
$\sigma ^2$
is used to model the surrogate function. The GPs define a distribution over potential functions, with the mean representing predictions and the variance quantifying uncertainties for each sample point within the design space.
To interpolate between any set of points using this framework, a covariance function is required to quantify the correlation between points and construct the covariance matrix. Assuming a constant mean in this approach, the correlation between points is expressed as the correlation between the errors
$\epsilon (\boldsymbol{x}_i)$
and
$\epsilon (\boldsymbol{x}_{\!j})$
, defined as
\begin{equation} \text{corr}(\epsilon (\boldsymbol{x}_i), \epsilon (\boldsymbol{x}_{\!j})) = \exp \left (-d(\boldsymbol{x}_i, \boldsymbol{x}_{\!j})\right ), \end{equation}
where
$d(\boldsymbol{x}_i, \boldsymbol{x}_{\!j})$
is a distance metric defined as
\begin{equation} d(\boldsymbol{x}_i, \boldsymbol{x}_{\!j}) = \sum _{h=1}^m \theta _h \left | x_{hi} - x_{hj} \right |^{p_h}, \end{equation}
with
$\theta _h$
representing the correlation coefficients and
$p_h$
governing the smoothness of the distance function for each dimension
$h$
. This distance function encodes the spatial correlation between the data points and is integral to inferring the function’s value at unobserved locations.
The model described is commonly referred to as the `design and analysis of computer experiments (DACE)’ stochastic process model. The DACE model includes
$m + 2$
parameters: the mean
$\mu$
, variance
$\sigma ^2$
and the hyperparameters
$\theta _1, \ldots , \theta _m$
. These parameters are estimated by maximising the likelihood of the observed data, with the maximisation process adjusting the parameters to achieve the best fit to the data while accounting for the spatial correlations.
The EI criterion estimates the potential gain in the minimum value of the objective function when sampling a particular point in the design space. The EI function guides the optimisation by selecting the next sampling point where the EI over the current minimum is maximised, balancing exploration and exploitation. The EI at a point
$\boldsymbol{x}$
is given by the expected value of the improvement, expressed as
\begin{align} EI(\boldsymbol{x}) = \int _{-\infty }^{f_{\textit{min}}} \left (f_{\textit{min}} - \hat {y}(\boldsymbol{x})\right ) p(\hat {y}(\boldsymbol{x})) \, {\rm d}\hat {y}(\boldsymbol{x}), \end{align}
where
$ f_{\textit{min}}$
denotes the current optimal value of the surrogate model, and
$ \hat {y}(\boldsymbol{x})$
is the predicted value at the candidate point
$\boldsymbol{x}$
, with
$ p(\hat {y}(\boldsymbol{x}))$
being the probability density function of the Gaussian distribution. This is schematically shown in figure 6.
The expectation is typically evaluated in a closed-form expression, which can be derived based on the properties of the normal distribution. The details of this derivation can be found in Jones et al. (Reference Jones, Schonlau and Welch1998).
The objective of the optimisation was to minimise the drag coefficient
$ C_D$
of the body. An initial 10 samples of the values of
$A$
,
$\beta$
and
$C$
were randomly selected from within the allowed ranges using Latin hypercube sampling (see figure 8). The optimisation process was then carried out over 8 iterations, with each iteration consisting of 10 function evaluations. In the present context, one function evaluation corresponds to a single high-fidelity WRLES simulation of the geometry defined by a given combination of
$(A/L,\,C/L,\,\beta )$
. Thus, 10 function evaluations per iteration meant that 10 geometries were simulated in parallel before the surrogate model was updated. This set-up was chosen to balance the computational cost with the need for a robust optimisation process. The number of function evaluations per iteration was selected to optimise the time required for the simulations (which can all be run in parallel) and the frequency of model updates. Since the true response surface is unknown, formal error metrics for surrogate accuracy are of limited interpretive value, and errors evaluated only at sampled points can be misleading. Instead of enforcing a rigid numerical stopping rule, convergence was assessed adaptively by monitoring the stabilisation of the predicted optimum, the evolution of the surrogate model as additional data were incorporated and the consistency of the emerging trends with physical intuition near the boundaries of the design space.
To further evaluate robustness, additional off-design simulations were performed outside the optimisation loop. These independent results were consistent with the surrogate predictions, reinforcing confidence in both the convergence and the reproducibility of the identified optimum.
3. Results
This section presents the results of the optimisation process, alongside an analysis of the parameter dependencies inferred from the surrogate models, supported by additional sensitivity studies. The section begins with the drag optimisation results, including the final surrogate models that illustrate the dependence of the drag coefficient on the three geometric parameters, and the predicted optimal configuration. Based on trends observed in the surrogate models, further simulations were conducted in regions exhibiting similar behaviour, to investigate the physical mechanisms driving the observed trends. Finally, surrogate models for the lift coefficient are introduced, followed by the Pareto front derived from the combined optimisation of both drag and lift coefficients.
3.1. Drag optimisaiton
3.1.1. General trends
The drag coefficient in the absence of any roof extension or taper is
$C_D = 0.402$
. The optimisation process converges to a drag coefficient of
$C_D = 0.376$
after 60 simulations, without human-guided interventions, as shown in figure 7. This corresponds to a
$6\,\%$
reduction. The optimal parameter set is
$[A/L, \beta , C/L] = [0.07, 6.7^\circ , 0.06]$
, which aligns closely with the minimum predicted by the surrogate model,
$[A/L, \beta , C/L] = [0.07, 7^\circ , 0.07]$
(figure 9), yielding a similar drag coefficient of
$C_D = 0.377$
.
Figure 7. Plot of the convergence of the drag coefficient with the number of iterations, where each iteration corresponds to 10 simulations, and optimal geometry overlaid with the pressure field superimposed on the dot product of the inward-pointing surface normal vector.
Figure 8. Distribution of evaluated points in the design space, coloured by drag coefficient (
$C_D$
) values. The colour scale is capped at 0.4 to emphasise regions where
$C_D$
is lower than that of the baseline geometry. Red points indicate the initial Latin hypercube sampling points. The plots represent a three-dimensional parallel projection of the design space, with each subplot depicting a different two-dimensional projection of the multidimensional optimisation domain.
Figure 9. Top: section plot of the surrogate model for angles of incidence between 0 and
$20^\circ$
at
$\beta =7^\circ$
,
$C/L=0.07$
,
$A/L=0.07$
, respectively. Bottom: section plot of the surrogate model for angles of incidence between
$-20^\circ$
and
$0 ^\circ$
degrees at
$\beta =-20^\circ$
,
$C/L=0.07$
,
$A/L=0.07$
, respectively. Note the different colour bars.
Compared with a conventional parametric sweep, the Bayesian optimisation framework delivered substantial computational savings. For example, a uniform grid study at a resolution of
$\varDelta (A/L)=0.01$
,
$\varDelta (C/L)=0.01$
and
$\Delta \beta =1^\circ$
(restricted to non-negative
$\beta$
for fairness) would require of the order of 1000 simulations, corresponding to approximately 100 000 core hours. By contrast, the present optimisation converged in only 60 simulations, consuming about 6000 core hours – more than an order of magnitude reduction.
Figure 8 shows the distribution of evaluated design points within the domain space. The plot presents parallel projections of the three-dimensional design space when viewed along the extension length
$ C$
axis and the taper penetration distance
$ A$
axis, respectively. Two notable trends emerge: most optimal configurations occur within the
$\beta$
-parameter range of
$5^\circ$
to
$10^\circ$
, with
$C_D$
decreasing as
$A$
and
$C$
increase. Furthermore, no outliers were identified.
The final surrogate model, created using all evaluated design points, is visualised in figure 9. Since the optimisation involves three parameters, a direct visualisation of the full surrogate model is not feasible. Instead, two-dimensional cross-sections are presented, illustrating the dependence of
$C_D$
on two parameters while holding the third constant.
Distinct behaviours arise for positive and negative angles of incidence, leading to separate subfigures. The top three subfigures correspond to positive angles of incidence, with the third parameter set to the optimal configuration for
$C_D$
minimisation. Conversely, for negative angles of incidence, the third parameter corresponds to the configuration yielding maximum
$C_D$
.
For positive angles of incidence, the drag coefficient exhibits high sensitivity to
$\beta$
. A minimum occurs within the range of
$5^\circ$
to
$10^\circ$
, with its exact location depending on
$A$
and
$C$
. The influence of
$\beta$
becomes more pronounced as
$A$
and
$C$
increase, showing a near-linear relationship with
$C_D$
. Among the parameters,
$C$
has a greater effect than
$A$
. The predicted minimum at
$[A/L, \beta , C/L] = [0.07, 7^\circ , 0.07]$
corresponds to
$C_D = 0.377$
. The design space yielding improvements in
$C_D$
spans
$[A/L, \beta , C/L] = [0{-}0.07, 0{-}15^\circ , 0{-}0.07]$
.
For negative angles of incidence,
$C_D$
increases uniformly as the angle magnitude grows. Parameters
$A$
and
$C$
again show a near-linear effect, with
$C$
having a greater influence. The maximum drag coefficient,
$C_D = 0.665$
, occurs at the domain boundary,
$[A/L, \beta , C/L] = [0.07, -20^\circ , 0.07]$
.
3.1.2. Partial dependence studies
To further investigate the individual influence of each parameter, partial dependence studies were conducted. They involved performing additional simulations with equal spacing, varying only one parameter while keeping the other two constant. The values of the fixed parameters were selected based on observations from the general trends discussed in § 3.1.1, ensuring that the partial dependence plots accurately reflect the trends observed in the analysed region.
The surrogate models, overlaid with the results of the new simulations, are presented in this section. These results are further interpreted in terms of the underlying physical mechanisms.
3.1.2.1. Taper penetration distance
$\boldsymbol{A}$
Three distinct behaviours can be seen in the partial drag dependence on the taper penetration distance
$ A$
(figure 9): (i) for negative angles of incidence, (ii) for positive angles of incidence without experiencing separation on the taper/extension and (iii) for positive angles when a separation is present.
In the first region, drag increases with taper penetration distance, primarily due to the effective base area of the body increasing, and enhanced upwash. No distinct physical phenomena associated with the dependence of parameter
$ A$
were identified in this study.
To analyse the pre-separation region, a representative angle of incidence of
$\beta =5^\circ$
and an extension length of
$C/L=0.03$
were selected. Figure 10 shows the surrogate model, along with partial dependence simulation results.
Figure 10. Partial dependence of
$C_D$
on
$ A$
. Left: surrogate model plot at
$\beta = 5^\circ$
and
$C/L = 0.03$
, incorporating the
$1.96\sigma$
confidence interval. Simulation data points are overlaid for (a)
$A/L = 0$
, (b)
$A/L = 0.015$
, (c)
$A/L = 0.03$
, (d)
$A/L = 0.045$
and (e)
$A/L = 0.06$
. Right: drag coefficient variation with roof extension
$A$
for geometries extending outside the optimisation domain. Red points represent data within the domain, while blue crosses indicate points outside.
The surrogate model, which closely aligns with conducted simulations, indicates a nearly linear decrease in the drag coefficient with taper penetration distance within the analysed range of
$A$
, with a total decrease of approximately 2.5 %.
Before discussing the flow-physics trends in detail, it is important to note that, although the visualisations presented below focus on mid-plane and base-surface views for clarity, the underlying analysis was carried out using the complete three-dimensional LES fields. Both instantaneous and time-averaged three-dimensional velocity and vorticity distributions were examined to confirm that these two-dimensional sections capture the dominant physical mechanisms. For the quasi-two-dimensional roof-extension geometry studied here, the essential wake behaviour was found to be well represented in the symmetry plane. Therefore, the mid-plane and surface plots shown in the following sections provide a clear yet physically faithful visualisation of the governing flow structures.
Figure 11 shows the conditionally time-averaged results of the simulations corresponding to figure 10(a–e). The first column presents the pressure field multiplied by the dot product of the surface normal vector (inward pointing) and the streamwise normal vector, which represents drag contributions per unit surface area. The second column displays streamlines at the centreline. Since the bi-modal positions are symmetric (Dalla Longa et al. Reference Dalla Longa, Evstafyeva and Morgans2019), the centreline effectively captures the wake structure and location.
Figure 11. Left: conditionally time-averaged pressure field superimposed on the dot product of the surface normal vector (inward pointing) and the streamwise normal vector. Right: streamlines at the centreline of the conditionally time-averaged velocity, depicting the wake. Both columns correspond to
$\beta = 5^\circ$
,
$C/L = 0.03$
, and: (a)
$A = 0.0$
, (b)
$A/L = 0.015$
, (c)
$A/L = 0.03$
, (d)
$A/L = 0.045$
and (e)
$A/L = 0.06$
.
The pressure contours indicate that the average base pressure increases, i.e. base-pressure recovery improves, with taper penetration distance. As no significant change in the streamlines pattern is observed, this pressure recovery is interpreted as classical diffuser-induced pressure recovery associated with flow deceleration. This interpretation is supported by the average pressure coefficient at the separation point increasing from
$C_P \sim -0.22$
in simulation (a) to
$C_P \sim -0.21$
in simulation (e). This trend, combined with the decreasing base area as taper penetration distance increases, aligns with the classical interpretation of the boat-tailing effect. The analysis of streamlines in the wake’s middle region, together with the location of the saddle point, shows that downwash increases with the taper penetration distance
$A$
. Since the geometry corresponding to
$A/L = 0$
experiences a slightly upwash-dominated wake, i.e. streamlines in the middle region create a negative angle with the horizontal, the net effect is that low-drag geometries correspond to a `balanced’ wake, i.e. one that is neither downwash nor upwash dominated. This finding aligns with the conclusions of Urquhart et al. (Reference Urquhart, Varney, Sebben and Passmore2020b
). Nevertheless, establishing causality or dependence in this case remains challenging.
The analysis of the simulations presented in figure 11(a–e) reveals a trend in the spanwise balance of the wake. Specifically, as
$ A$
increases, the spanwise pressure gradient decreases, strongly suggesting that the wake undergoes a process of re-symmetrisation. Quantitatively, the non-dimensionalised spanwise pressure gradient decreases from approximately 0.15 to 0.10 across the examined range of
$A/L$
, confirming the progressive weakening of lateral asymmetry. This trend aligns with previous studies investigating the sensitivity of bi-modality to passive-control methods. Evrard et al. (Reference Evrard, Cadot, Herbert, Ricot, Vigneron and Délery2016) demonstrated wake re-symmetrisation in cavity flow control, showing that increasing cavity depth progressively reduced the magnitude of the spanwise pressure gradient, eventually leading to full re-symmetrisation. Similarly, Bonnavion & Cadot (Reference Bonnavion and Cadot2019) found that rear base tapers can also promote wake re-symmetrisation. To investigate this effect further, additional simulations were conducted with
$ A$
values extending beyond the optimisation domain (see figure 10).
Figure 12 presents the probability density distributions of the side force coefficient for selected simulations both within and beyond the optimisation domain. This distribution serves as an indicator of the presence and intensity of bi-modality. As
$ A$
increases, the probability density peaks of the side force coefficient magnitudes shift toward zero, achieving full re-symmetrisation at
$ A/L \sim 0.22$
. This confirms not only the presence of re-symmetrisation but also its continuous nature, aligning with previous studies (Evrard et al. Reference Evrard, Cadot, Herbert, Ricot, Vigneron and Délery2016; Bonnavion & Cadot Reference Bonnavion and Cadot2019).
Figure 12. Variation of the probability density histograms of side force coefficients with
$A$
illustrating the re-symmetrisation process for
$\beta = 5^\circ$
and
$C/L = 0.03$
, overlaid with Gaussian mixture model best fit curves.
The presence of the re-symmetrisation process is particularly relevant to this study, as bi-modality has been seen to contribute between 4 % and 9 % of the total pressure drag (Grandemange, Gohlke & Cadot Reference Grandemange, Gohlke and Cadot2014) for geometries of this kind. This implies that mitigating bi-modality through re-symmetrisation serves as an additional, continuous drag-reduction mechanism. As shown by Evrard et al. (Reference Evrard, Cadot, Herbert, Ricot, Vigneron and Délery2016), the benefits of wake symmetrisation plateau once full re-symmetrisation is achieved. A similar trend is observed in figure 10, where the gradient of the drag coefficient diminishes at
$ A/L \sim 0.15{-}0.2$
, coinciding with the onset of wake re-symmetrisation. Although the wake dynamics remains complex, and establishing a direct causal link is challenging, the observed trend suggests that this relationship holds.
A key factor influencing wake re-symmetrisation appears to be the base aspect ratio. Grandemange et al. (Reference Grandemange, Gohlke and Cadot2013a
) found that bi-modality typically disappears entirely when the base aspect ratio falls below
$ H^* = H/W \lt 0.6$
. This threshold suggests that as the wake’s vertical extent is reduced, it becomes more stable, suppressing bi-modality and promoting a symmetric state. This aligns with the results of Bonnavion & Cadot (Reference Bonnavion and Cadot2019), who showed that aspect ratio plays a dominant role in governing re-symmetrisation. As the aspect ratio decreases, the wake transitions gradually toward a more stable symmetric state. This effect is consistent with the influence of increasing the taper penetration distance
$ A$
, which effectively reduces the base aspect ratio. Thus, as indicated in figure 13, reducing the base aspect ratio is likely to promote re-symmetrisation. However, the data presented by Bonnavion & Cadot (Reference Bonnavion and Cadot2019) also suggest that re-symmetrisation is influenced by the base vertical pressure gradient. Specifically, when the vertical pressure gradient is low, achieving re-symmetrisation requires a more substantial reduction in aspect ratio, whereas for larger gradients, re-symmetrisation occurs more rapidly. This observation aligns with the findings of Barros et al. (Reference Barros, Borée, Cadot, Spohn and Noack2017), who conducted a sensitivity study on underbody disturbances that modified vertical pressure distributions. Their results showed that bi-modality persisted when vertical pressure gradients were small but disappeared for larger gradients. Consequently, re-symmetrisation may occur at a higher aspect ratio than the threshold predicted by Grandemange et al. (Reference Grandemange, Gohlke and Cadot2013a
). Additionally, the geometry corresponding to
$ A/L = 0.22$
had the taper slightly extending beyond the base (
$ C/L = 0.3$
), and as demonstrated by Evrard et al. (Reference Evrard, Cadot, Herbert, Ricot, Vigneron and Délery2016), cavity-like control methods can also contribute to re-symmetrisation. As a result of these secondary effects, the geometry that achieved re-symmetrisation (
$[A/L, \beta , C/L] = [0.22, 5^\circ , 0.03]$
) had an aspect ratio of
$ 0.69$
, significantly higher than the threshold of
$ 0.6$
identified by Grandemange et al. (Reference Grandemange, Gohlke and Cadot2013a
).
Figure 13. Figure adapted from Grandemange et al. (Reference Grandemange, Gohlke and Cadot2013a
), depicting the domains of instability development observed in experiments, as a function of aspect ratio (
$H^* = H/W$
) and ground clearance (
$C^* = C/W$
). The star marks the aspect ratios of the base geometry.
The surrogate model depicting post-separation behaviour, along with the conditionally time-averaged results, is shown in figure 14. The model predictions align closely with the simulation results, showing an initial decrease in the drag coefficient followed by an increase. The minimum corresponds to a
$0.75\,\%$
reduction in the drag coefficient.
Figure 14. Partial dependence of
$C_D$
on
$ A$
. Left: surrogate model plot at
$\beta = 15^\circ$
and
$C/L = 0.03$
, incorporating the
$1.96\sigma$
confidence interval. Right: conditionally time-averaged pressure field superimposed on the dot product of the surface normal vector (inward-pointing) and the streamwise normal vector for simulations with (a)
$A/L = 0.015$
, (b)
$A/L = 0.03$
and (c)
$A/L = 0.06$
.
Increasing the taper penetration distance leads to effects similar to those observed in the pre-separation case. However, due to the higher angle of incidence, the induced downwash is more pronounced, resulting in a downwash-dominated wake that is predominantly influenced by the upper recirculation region. Additionally, as the taper penetration distance increases, the angled surface over which an adverse pressure gradient exists grows, leading to more influential separation effects. Initially, the angled surface does not provide enough surface for the shear layer to reattach (
$A/L \lt 0.045$
). The shear layer reattaches at
$A/L \sim 0.045$
, forming a separation bubble and a tip vortex, which result in lower pressure at the angled surface and enhanced downwash. That low-pressure region, induced by the separation effects, outweighs the base-pressure gains at
$A/L \sim 0.045$
. Furthermore, the wake undergoes re-symmetrisation at
$A/L = 0.06$
. The earlier onset of re-symmetrisation can be attributed to two factors: first, the higher angle of incidence leads to a more significant reduction in base aspect ratio; second, it generates a stronger vertical pressure gradient, both of which promote wake re-symmetrisation.
3.1.2.2. Angle of incidence
$\boldsymbol{\beta}$
The dependence on the taper/extension angle of incidence,
$ \beta$
, exhibits two distinct behaviours. For angles below that for which separation occurs on the angled surface, the drag coefficient decreases consistently as
$ \beta$
increases. However, once separation occurs, the low pressure on the angled surface outweighs the gains from base-pressure recovery and this trend reverses, with further increases in
$ \beta$
leading to higher drag coefficient. This dependence is analysed in two stages: first, for positive angles of incidence followed by negative angles. In both cases, surrogate model cuts at
$A/L = 0.07$
and
$C/L = 0.07$
are examined, as they correspond to the conditions yielding the minimum and maximum drag coefficients, respectively.
Moreover, as demonstrated in the previous section, the drag penalties associated with the low-pressure region on the angled surface are considerably lower for shorter geometries. Thus, short tapers at
$A/L = 0.03$
and
$C/L = 0.0$
are also investigated.
Figure 15 illustrates the dependence for positive angles of incidence at
$A/L = 0.07$
and
$C/L = 0.07$
. The simulation results closely align with predictions from the surrogate model, accurately identifying the minimum drag coefficient of
$C_D = 0.375$
at
$\beta = 7^\circ$
, which corresponds to a
$7\,\%$
drag reduction compared with the zero angle of incidence geometry. These results closely align with experimental studies by Perry et al. (Reference Perry, Pavia and Passmore2016), Pavia et al. (Reference Pavia, Passmore and Gaylard2016) and Fan et al. (Reference Fan, Parezanović, Fichera and Cadot2024).
Figure 15. Partial dependence of
$C_D$
on
$ \beta$
. Left: surrogate model plot at
$A/L = 0.07$
,
$C/L = 0.07$
and
$\beta \gt 0^\circ$
, incorporating the
$1.96\sigma$
confidence interval. Right: conditionally time-averaged pressure field superimposed on the dot product of the surface normal vector (inward-pointing) and the streamwise normal vector for simulations with (a)
$\beta = 0^\circ$
, (b)
$\beta = 7^\circ$
and (c)
$\beta = 15^\circ$
.
Interestingly, the zero angle of incidence geometry results in a higher drag coefficient than the baseline geometry. This outcome suggests that, while the extension provides a cavity-like effect by increasing the distance between the top part of the toroidal vortex structure and the base, other factors dominate. This observation underscores the complexity of the interactions and will be further analysed as part of the partial dependence study of parameter
$C$
.
An increase in the angle of incidence,
$\beta$
, induces several effects: diffuser-induced pressure recovery, base-size reduction, downwash generation, spanwise wake re-symmetrisation and separation on the angled surface (only for angles greater than
$7^\circ$
). The interplay of these effects determines the drag coefficient values. Examining the simulations in figure 15(a, b), it is evident that the initially asymmetric, upwash-dominated wake transitions to a symmetric state with significantly higher average base pressure and reduced upwash. The optimal geometry exhibits minimal separation, with the separation region small enough that it does not result in very low pressure on the angled surface.
For larger angles of incidence, the wake symmetry persists, and base pressure further increases. However, the adverse effects of the separation bubble and tip vortices create low-pressure regions that ultimately outweigh the gains in base pressure, leading to a net increase in drag.
Figure 16 shows the dependence of the drag coefficient on the angle of incidence for short tapers without the extension, i.e. geometries corresponding to
$ A/L = 0.03$
and
$ C = 0$
. The drag follows a similar trend to that observed for extensions at
$ A/L = 0.07$
and
$ C/L = 0.07$
, initially decreasing with
$ \beta$
, reaching a minimum at
$ \beta \sim 5^\circ$
corresponding to
$ C_D = 0.396$
and then increasing to
$ C_D = 0.399$
at
$ \beta \sim 10^\circ$
. The increase plateaus for
$ \beta \gt 10^\circ$
. Notably, the increase appears to stem from a different mechanism than that observed in the geometries analysed at
$ A/L = 0.07$
and
$ C/L = 0.07$
. Since the shear layer does not have sufficient space to reattach on the angled surface, no separation bubble is formed. As a result, the low-pressure region that would have been induced by the separation bubble is avoided. Nevertheless, the short taper has a substantial effect on the vertical balance, eventually leading to a wake that is highly dominated by downwash and, thus, an excessive growth of the top recirculation region.
Figure 16. Partial dependence of
$C_D$
on
$ \beta$
. Surrogate model plot at
$A/L = 0.03$
,
$C/L = 0.0$
and
$\beta \gt 0^\circ$
, incorporating the
$1.96\sigma$
confidence interval. Conditionally time-averaged pressure fields are superimposed on the dot product of the surface normal vector (inward pointing) for simulations with (a)
$\beta = 5^\circ$
, (b)
$\beta = 10^\circ$
and (c)
$\beta = 15^\circ$
. Note: as these conditions are far from the optimal design region and correspond to only minor variations in drag compared with other tested configurations, the initial surrogate model failed to capture the observed behaviour. The surrogate model shown in the plot has therefore been updated with the additional simulations.
Interestingly, this suggests that in certain configurations, the optimal angle of incidence may correspond to post-separation geometries, i.e. when the benefits of upwash/downwash effects outweigh the drawbacks associated with separation. This reasoning could explain the findings of Perry et al. (Reference Perry, Pavia and Passmore2016) and Pavia et al. (Reference Pavia, Passmore and Gaylard2016), who investigated short tapers and often concluded that optimal angles of incidence were above
$ 12^\circ$
. However, this was only the case for specific configurations, e.g. symmetric side tapers. For top tapers alone, their optimal incidence angles closely aligned with the findings in this study.
The observed maximum drag reduction of
$ \Delta C_D = 0.006$
, when scaled to account for the difference in the taper size, is consistent with the values reported in earlier studies on short tapers (Littlewood & Passmore Reference Littlewood and Passmore2010; Pavia et al. Reference Pavia, Passmore and Gaylard2016; Perry et al. Reference Perry, Pavia and Passmore2016).
The dependence for negative angles of incidence represents a close-to-linear increase of the drag coefficient with the magnitude of
$\beta$
, with an average rise of approximately
$35\,\%$
per
$10^\circ$
, reaching
$C_D = 0.678$
at
$\beta = -20^\circ$
. Conditionally time-averaged results depicting this behaviour are shown in figure 17. All simulations exhibit symmetrised, highly upwash-dominated wakes, with the wake primarily influenced by the lower part of the toroidal vortex structure. As
$|\beta |$
increases, the strength of the lower recirculation region intensifies, resulting in a decrease in base pressure. This highlights that even for negative angles of incidence, the same underlying mechanisms remain active, including the influence of the vertical pressure gradients on wake re-symmetrisation. At
$\beta = -15^\circ$
, a separation bubble is observed, which persists at lower angles of incidence as well. Even though it contributes to energy losses, no significant changes in the drag coefficient increases were noticed.
Figure 17. Left: rear view of the conditionally time-averaged pressure field overlaid on the dot product of the surface normal vector (inward pointing) and the streamwise normal vector. Middle: front view of the same field. Right: centreline streamlines of the conditionally time-averaged velocity, illustrating the wake. All plots correspond to
$A/L = 0.07$
,
$C/L = 0.07$
, with: (a)
$\beta = -5^\circ$
, (b)
$\beta = -10^\circ$
and (c)
$\beta = -15^\circ$
.
3.1.2.3. Extension length
$\boldsymbol{C}$
The dependence on the extension length,
$ C$
, exhibits behaviour analogous to that observed for parameter
$A$
, encompassing three distinct regimes: (i) negative angles of incidence, (ii) positive angles of incidence without separation on the extension and (iii) positive angles of incidence with a separation bubble present. As with the taper penetration distance, the region corresponding to negative angles of incidence is excluded from further analysis due to the absence of distinct physical phenomena or dependencies.
Following the methodology used for the taper penetration distance, a partial dependence analysis is performed at
$A/L = 0.03$
and
$\beta = 5^\circ$
(pre-separation regime) and
$A/L = 0.03$
and
$\beta = 15^\circ$
(post-separation regime). Figure 18 presents the surrogate model alongside the simulation runs conducted. The analysis demonstrates a continuous, close to linear, decrease in the drag coefficient within the studied region. The total drag reduction observed in this region is
$3\,\%$
, which is
$0.5$
percentage points higher than the drag reduction achieved by varying the taper penetration distance. However, the minimum drag coefficient achieved is larger than in the partial dependence study of the parameter
$A$
, as the geometry at
$C/L = 0.0$
and
$A/L = 0.03$
resulted in higher drag than that at
$A = 0$
and
$C/L = 0.03$
.
Figure 18. Partial dependence of
$C_D$
on
$ C$
. Top: surrogate model plot at
$A/L = 0.03$
and
$\beta = 5^\circ$
, incorporating the
$1.96\sigma$
confidence interval. Bottom: conditionally time-averaged pressure field superimposed on the dot product of the surface normal vector (inward pointing) and the streamwise normal vector, and streamlines at the centreline. Both columns correspond to
$A/L = 0.03$
and
$\beta = 5^\circ$
, and: (a)
$C/L = 0.0$
, (b)
$C/L = 0.03$
and (c)
$C/L = 0.06$
.
The conditionally time-averaged results for three representative simulations from the partial dependence study are shown in figure 18(a–c). Despite the geometric similarity between parameters
$C$
and
$A$
, the underlying physics governing their dependencies differ significantly. Both parameters lead to increased base-pressure recovery as they are increased, but the wake topology is modified in fundamentally distinct ways.
Firstly, as the extension length increases, the centre of the top recirculation region shifts further from the base, causing a vertical tilt in the toroidal vortex structure. The extension also expands the region where the top shear layer remains unaffected by the low pressure in the recirculation zone. Consequently, the total curvature of the top shear layer decreases, shifting the wake’s saddle point upwards. To compensate for this upward shift, the size of the bottom recirculation region increases. This behaviour is seen in both the streamline patterns and the base-pressure contours. The strengthening of the bottom recirculation region is hypothesised to counteract the pressure recovery gains and the displacement of the top recirculation region. This interaction may explain the deviation from the surrogate model trend observed in the simulations corresponding to
$C/L = 0.045$
and
$C/L = 0.06$
, which exhibit a plateau-like behaviour in drag coefficient despite increasing extension length.
This effect is also observed in figure 15(a), where the geometry with an extension at zero angle of incidence exhibits a higher drag coefficient than the base geometry. The lower pressure at the bottom of the base outweighs the pressure gains in the top region, resulting in a net increase in drag. These results emphasise the importance of a `balanced’ wake. While direct correlations remain challenging to establish, unbalanced wakes seem to lead to one recirculation region becoming significantly stronger, negating the benefits of reducing the size of the opposing region. In the context of extension length, this suggests that further improved performance may be possible with a symmetric system consisting of top and bottom extensions.
Secondly, increasing the extension length does not directly reduce the base area. However, the downward movement of the separation point, and consequently the top recirculation centre, leads to a reduction in wake size, resulting in a similar effect to base-size reduction. Additionally, the extension serves as a physical boundary that opposes fluid motion in the corresponding region, thereby increasing pressure in that area.
Spanwise wake re-symmetrisation was not observed within the analysed region. However, as for parameter
$ A$
, additional simulations were conducted to test the hypothesis of wake re-symmetrisation. Simulations at (a)
$ C/L = 0.08$
, (b)
$ C/L = 0.10$
and (c)
$ C/L = 0.12$
resulted in complete wake re-symmetrisation. These results highlight the substantial difference in the influence of parameters
$ A$
and
$ C$
on the wake dynamics. Increasing the extension length induces a higher vertical pressure gradient than increasing the taper penetration length, which positively contributes to wake re-symmetrisation. Moreover, it enhances the cavity-like contribution to re-symmetrisation, as described by Haffner et al. (Reference Haffner, Borée, Spohn and Castelain2020). The combination of these two effects leads to significantly earlier re-symmetrisation compared with that observed for taper penetration length.
All three simulations – (a)
$ C/L = 0.08$
, (b)
$ C/L = 0.10$
and (c)
$ C/L = 0.12$
– resulted in
$ C_D = 0.384$
. The saturation of the drag coefficient values may be attributed to either the saturation of the re-symmetrisation effect or the vertical balance effects described above.
Figure 19 illustrates the partial dependence on
$ C$
in the post-separation regime at
$\beta = 15^\circ$
. Increasing the extension length leads to greater downwash, a downward and streamwise shift in the location of the top recirculation region and higher base pressure. Additionally, as for parameter
$A$
, the shear layer does not experience reattachment when
$C/L$
is lower than a critical value, which in this case is
$0.03$
. Similar to the behaviour observed for parameter
$A$
, an initial decrease in the drag coefficient is followed by an increase, as the effects of the separation bubble begin to outweigh the improvements in base-pressure recovery. In contrast to the results shown in figure 18, the wake corresponding to
$C/L = 0.06$
experiences re-symmetrisation. This further supports the theory that effective height reduction and vertical pressure gradients are the dominant factors leading to wake re-symmetrisation.
Figure 19. Partial dependence of
$C_D$
on
$ C$
. Top: surrogate model plot at
$A/L = 0.03$
and
$\beta = 15^\circ$
, incorporating the
$1.96\sigma$
confidence interval. Bottom: conditionally time-averaged pressure field superimposed on the dot product of the surface normal vector (inward pointing) and the streamwise normal vector, and streamlines at the centreline. Both columns correspond to
$A/L = 0.03$
and
$\beta = 15^\circ$
, and: (a)
$C/L = 0.0$
, (b)
$C/L = 0.045$
and (c)
$C/L = 0.06$
.
In contrast to the pre-separation regime, for the post-separation regime, increasing the extension length does not lead to the growth of the bottom recirculation region. Instead, it results in an expansion of the top recirculation region. This suggests that the downwash induced by the roof extension outweighs the effects associated with the streamwise displacement of the top recirculation region. These findings underscore the complexity of the interactions involved and the necessity of analysing each case individually.
3.2. Surrogate model for the dependence of lift on the roof extension parameters
Figure 20 presents the surrogate model of the lift coefficient, derived from the results obtained during the drag optimisation process. The top two subfigures correspond to the maximum values of extension length and taper penetration distance, respectively, depicting high downforce (negative lift) conditions. The minimum lift coefficient occurs at the edge of the domain at
$[A/L, \beta , C/L] = [0.07, -20^\circ , 0.07]$
, yielding a downforce increase of approximately
$\Delta C_L = 0.48$
compared with the baseline lift coefficient of
$C_L = -0.082$
. While a large downforce can be desirable, especially in motorsport, such an extension would lead to significantly higher drag, which would make it less desirable in real-world applications. To provide a more feasible perspective, the bottom subplots illustrate the lift dependencies for short roof extensions, demonstrating that a substantial downforce increase of approximately
$\Delta C_L = 0.15$
can be achieved at
$[A/L, \beta , C/L] = [0.01, -20^\circ , 0.01]$
.
Figure 20. Section plot of the surrogate model of the lift coefficient for short and long extensions.
3.3. Drag–lift Pareto frontier
To explore the trade-off between low-drag and low lift conditions, a Pareto frontier was constructed, showcasing the optimal solutions of a modified optimisation function
\begin{equation} f = \alpha C_L + (1-\alpha ) C_D ,\end{equation}
as shown in figure 21. This Pareto frontier is based on data collected during the drag optimisation study. In the context of road vehicles, lift often serves as a constraint rather than a parameter to be optimised. That is, the goal is to achieve the lowest possible drag without exceeding specified lift limits. To determine the optimal geometry at a given lift coefficient constraint, it is possible to identify the value of
$ \alpha$
at the limiting lift coefficient and then extract the associated geometric parameter values.
Figure 21. Pareto frontier illustrating the trade-off between drag and lift coefficients in the multi-objective optimisation study. Left: drag versus lift coefficient as the cost function weight
$\alpha$
increases, highlighting how the optimisation shifts from drag minimisation to lift (downforce) prioritisation. Each point corresponds to an incremental increase in
$\alpha$
by 0.1, starting from
$\alpha = 0$
(lower-left point in the plot). Right: corresponding geometric parameter values
$(A, \beta , C)$
that yield the optimal solutions for each
$\alpha$
.
It is worth noting that constructing surrogate models directly from an optimisation process targeting the entire Pareto frontier, e.g. via hypervolume-based acquisition strategies, would likely improve the robustness and coverage of the surrogate model across the design space. However, this approach requires a significantly higher computational cost and is therefore suggested as a direction for future work.
Figure 21 shows that, as
$\alpha$
increases, the optimisation prioritises downforce over drag reduction, saturating at
$\alpha \sim 0.6$
. At
$\alpha = 0$
, the parameters correspond to the values determined by the drag optimisation (§ 3.1). With increasing
$\alpha$
, the angle of incidence increases, reducing downwash and thus increasing downforce, while parameters
$A$
and
$C$
decrease to prevent significant drag increases. Both
$A$
and
$C$
reach values close to zero for
$\alpha \sim 0.2$
, representing a minimum. Once lift starts dominating the optimisation function,
$A$
and
$C$
begin to increase.
Notably, the weighting parameter
$\alpha$
can also be determined based on specific performance objectives, such as mathematical modelling of race time for a given track, incorporating the effects of drag and lift.
Figure 22 illustrates the surrogate model of the cost function at
$\alpha \sim 0.3$
, representing an optimisation that balances downforce with drag. The optimal design corresponds to negative angles of attack at low values of parameters
$A$
and
$C$
. Physically, negative angles of attack promote increased upwash and, consequently, greater downforce, while small extension dimensions minimise the impact on drag. Importantly, even small geometry add-ons, i.e.
$A/L = C/L \sim 0.01$
, correspond to significant reductions in the cost function, from 0.257 (baseline) to 0.241.
Figure 22. Section plot of the surrogate model representing the lift–drag optimisation.
4. Conclusions
In this study, a drag optimisation analysis of a roof extension parameterised by three variables – taper penetration distance, angle of incidence and extension length – was conducted. The square-back Windsor body with wheels and zero yaw served as the baseline geometry. High-fidelity, WRLES were utilised in combination with Bayesian optimisation based on GPs (Kriging), guided by an EI criterion. The optimisation results, including trends captured by the surrogate models and the simulation outcomes, were analysed to establish links between parameter dependencies and the underlying physics. To further facilitate this analysis, partial dependence studies were performed, where two parameters were held constant while the third was varied.
The optimisation process converged after six iterations, corresponding to a total of 60 simulations, achieving a drag reduction of
$6.5\,\%$
. The surrogate model analysis revealed three distinct behavioural regions: negative angles of incidence, positive angles experiencing separation on the angled surface and positive angles without separation. The optimal design corresponded to an angle of incidence at the onset of separation (
$\beta = 7^\circ$
) and the other two parameters at their maximum values within the analysed domain,
$A/L = C/L = 0.07$
, where
$L$
represents the length of the baseline geometry.
The analysis of the conditionally time-averaged results revealed six distinct physical mechanisms whose interplay determined the performance benefits (see figure 23):
-
(i) Diffuser-induced pressure recovery – the angled extension acts as a diffuser, slowing the air down, and enhancing pressure recovery.
-
(ii) Base-size reduction – increasing either the taper penetration distance or the angle of incidence leads to a physical reduction of the base size.
-
(iii) Vertical wake asymmetry (‘wake balance’) – variations in the extension parameters alter the wake balance, influencing the dominance of the top or bottom recirculation region. It was observed that asymmetric wake structures, where one recirculation region dominates the other, often led to higher drag.
-
(iv) Separation effects – if there is flow separation on the angled surface of the taper/extension, a low-pressure region forms. For sufficiently large angled surface length, this leads to the formation of a separation bubble and tip vortices, further decreasing the pressure in this region. When this effect becomes sufficiently pronounced, it outweighs any gains from base-pressure recovery, resulting in increased drag.
-
(v) Physical relocation of the recirculation region cores – the extension shifts the position of the wake’s recirculation regions both in the streamwise direction (cavity-like effect), lowering the base pressure in that region, and in the vertical direction, influencing wake size and interactions between recirculation zones.
-
(vi) Spanwise wake re-symmetrisation – in certain configurations, the wake’s symmetry-breaking behaviour is reduced or disappears entirely, altering the base-pressure distribution and contributing to drag reduction. This mechanism is hypothesised to result mainly from the reduction in the effective body height, alongside the influence of vertical pressure gradients and the physical boundary imposed by the extension.
Figure 23. A schematic summarising the mechanisms of drag reduction: (a) diffuser-induced pressure recovery, (b) base-size reduction, (c) vertical wake asymmetry (‘wake balance’), (d) separation effects, (e) physical relocation of the recirculation region cores, (f) spanwise wake re-symmetrisation.
In general, increasing any of the studied parameters enhanced diffuser-induced pressure recovery and reduced the effective base area, consistent with classical boat-tailing effects. These changes led to increased base pressure and a corresponding reduction in drag. Furthermore, the reduction in base height lowered the base aspect ratio, which promoted spanwise wake re-symmetrisation. This effect was more sensitive to variations in the extension length than in the taper penetration distance, owing to the stronger influence of the former on vertical pressure gradients and the physical boundary imposed by the extension surface.
All three parameters – taper penetration (
$A$
), angle of incidence (
$\beta$
) and extension length (
$C$
) – typically improved aerodynamic performance up to the onset of flow separation on the angled surface. Beyond this point, the adverse pressure distribution associated with the separated shear layer began to outweigh the benefits of increased base pressure. This trade-off generally emerged just after separation was initiated. However, for configurations with short angled surfaces, the limited extent of the separation region reduced the associated drag penalty, allowing certain post-separation designs to remain aerodynamically effective. Such configurations may be advantageous in applications where the induced upwash or downwash is more important than the drag penalty arising from surface separation.
The most significant differences between the effects of taper penetration
$A$
and extension length
$C$
stemmed from how they altered the wake topology. Surfaces extending beyond the original base shifted the separation point downstream, displacing the top recirculation region and producing a classical cavity-like effect. This downstream displacement also modified the effective curvature of the shear layer, which in turn influenced the vertical balance of the wake. In some cases, this led to counterintuitive outcomes, such as the generation of net upwash in geometries that would nominally be expected to produce downwash. These findings ultimately suggest that, to fully exploit the cavity-like effect, both the top and bottom edges should be extended in order to ensure the vertical balance of the wake.
Surrogate models of the lift coefficient were also presented. Additionally, a Pareto frontier was constructed to identify optimal solutions from the multi-parameter optimisation involving both drag and lift coefficients. This frontier was derived based on the drag optimisation results and could be further improved by performing a dedicated Pareto front optimisation in future work. The analysis demonstrated how the optimal design parameters evolved with changes in the lift-to-drag trade-off and how they could be adapted to satisfy specific performance objectives, such as minimising drag under prescribed downforce constraints.
This work provides the first simultaneous optimisation of taper penetration, extension length and incidence angle for roof extensions, establishing how these parameters interact to shape wake topology and drag. By identifying their distinct and coupled roles, the study moves beyond cataloguing individual mechanisms to deliver system-level understanding and practical design rules, enabling more informed aerodynamic design. In parallel, the findings demonstrate the effectiveness of Bayesian optimisation in CFD-based aerodynamic optimisation, reinforcing its potential for future aerodynamic design applications.
Acknowledgements
This research was supported through an EPSRC iCASE studentship in collaboration with Jaguar Land Rover and UK Research and Innovation (UKRI). The simulations were performed using the ARCHER2 UK National Supercomputing Service (Beckett et al. Reference Beckett, Beech-Brandt, Leach, Payne, Simpson, Smith, Turner and Whiting2024). The authors acknowledge the use of these resources and thank the ARCHER2 team for their support.
Competing interests
Author 2 is employed by Jaguar Land Rover. The authors declare that they have no other competing interests.
Artificial intelligence
AI tools were used for language editing and clarity improvement only.
Open access
