Computational Aerodynamics Analysis of Non-Symmetric Multi-Element Wing in Ground Effect with Humpback Whale Flipper Tubercles

: The humpback whale ﬂipper tubercles have been shown to improve the aerodynamic coefﬁcients of a wing, especially in stall conditions, where the ﬂow is almost fully detached. In this work, these tubercles were implemented on a F1 front-wing geometry, very close to a Tyrrell wing. Numerical simulations were carried out employing the k − ω SST turbulence model and the overall effects of the tubercles on the ﬂow behavior were analyzed. The optimal amplitude and number of tubercles was determined in this study for this front wing where an improvement of 22.6% and 9.4% is achieved, respectively, on the lift and the L / D ratio. On the main element, the stall was delayed by 167.7%. On the ﬂap, the ﬂow is either fully detached, in the large circulation zone, or fully attached. Overall, in stall conditions, tubercles improve the downforce generation but at the cost of increased drag. Furthermore, as the tubercles are case-dependent, an optimal conﬁguration for tubercles implementation also exists for any geometry.


Introduction
In Formula 1 races, a tenth of a second difference in performance can result in huge financial losses, which could have otherwise be used to build the team's infrastructure, resource, and research. The greater is the gain, the more can be invested in research and development of these state-of-the-art machines. That is why every year, engineers and scientists work to provide drivers with the best cars. The performance enhancement is achieved on different parts of the car: the engine, chassis, and aerodynamics. As part of the aerodynamics of the car, the front and rear wings are the most important components, which generate downforce to enhance the grip, channel the air flow around the car, and reduce drag to improve speed. Optimizing such elements has significant impacts on the car's performances. The perfect wing would reduce the drag while improving the lift of the car. In that sense, the humpback whale flipper seems to be a masterpiece as, in stall conditions, these considerably enhance aerodynamic performances.
To optimize the flow behavior, many techniques exist [1,2]. One of them relies on the control of the turbulence that leads to the benefits such as lift increase, drag reduction, less turbulence production, etc. The control of the turbulence can be either passive or active. A passive control of the turbulence could be done by changing the initial geometry or adding devices breaking the turbulence such as the large-eddy-breakup devices that are known to decrease the drag by 30-40% [3,4]. Those devices are steady and require no energy to work. On the other hand, an active device requires energy as it is managed in real time. An active control of the turbulence can inject flow in the boundary layer to keep the flow attached or control the wall of the wing to synchronize with the turbulence production [5,6]. Here, a passive flow control device, which is a leading edge incorporating tubercles, as on the humpback whale flippers, was investigated ( Figure 1). These tubercles allow the humpback whale to be very mobile and get extra ability in turning. It was first investigated by Fish et al. in 1995 [7] and since then widely investigated [8,9]. This device reduces the drag due to lift on a wing [7,10]. The tubercles generate streamwise vortices which increase the momentum exchange within the boundary layer helping to keep the flow attached [11,12]. Using a wing with such a device provides better performance at higher angle of attack (see Figure 4 of [9]), delaying the stall angle by about 40%, and thus increasing the range of use of the wing. Many configurations on many wings have already been tested. Fish et al. [7] showed that a span section of the humpback whale flipper is very close to the NACA 63-021 profile. Other studies have shown that the addition of tubercles on the leading edge is free from this NACA profile. For example, Miklosovic et al. [14] employed the NACA 0020, Hansen et al. [11] employed the NACA 0021, and Van Nierop et al. [15] investigated the NACA 0018, but they arrived at similar conclusions. More recently, Castro and Rana [16] employed S1210 airfoil and analysed the aerodynamic and structural design for the 2022 formula one front wing configuration. In this study, such a passive flow control is used on a non-symmetric wing profile.
The frequency and the amplitude of the tubercles have an impact on the performance on a sinusoidal leading edge wing. In 2008, Van Nierop et al. [15] demonstrated from a mathematical point of view that a big amplitude and a small frequency lead to better performance. They obtained results that were reproduced experimentally by Hansen et al. [11] in 2010, however leading to the different conclusion that the better configuration is the one with the smallest amplitude and the biggest frequency. This shows another big aspect of the humpback whale flipper, the sweep and the taper. In fact, Fish et al. [7] considered the sweep-back and taper of the wing as also influencing the good performance of the humpback whale flipper. This was highlighted by both Hansen and Van Nierop. For wings with taper and sweep, the best configuration is the one with the biggest amplitude and the smallest frequency, as shown by Van Nierop [15]. For wings without taper and sweep, the best configuration is the one with smallest amplitude and the biggest frequency, as shown by Hansen [11].
Besides, a difference in performance is noticeable between full span and semi span wings, with sweep and taper using this passive control device. In fact, for a sinusoidal leading edge, the wing performance is reduced for a full span wing as it triggers an early separation of the flow. For a semi span wing, the spanwise stall progression is inhibited. This enables using a span wise wing in a bigger range of angle of attack, as shown by Miklosovic in 2007 [14]. With no taper and sweep on a wing, there is no difference at all between a full span and a semi span wing, as shown by Hansen in 2010 [11].
However, humpback whale flippers are made of 11 tubercles, which do not have a constant inter-tubercle distance or tubercle chord, as shown in Figure 1 and by the observations of Fish and Srinivras [7,9,18]. During their experiments in 2018, Srinivras demonstrated that a bio-mimicked distribution is more effective than a sinusoidal one, as it produces more lift and less drag. This is because the drag-to-lift enhancement comes from the biggest tubercles amplitude and position, the first and fourth [18]. Such models have only been tested for hydrodynamic purpose and only experimentally [18].
Zerihan's work of 2001 is the pioneer in this field. His motivation was to investigate a Tyrrell wing (026). Two different wings were used for experimental simulations: a single-element wing, the LS(1)-0413 MOD with slight modifications, and a two-element wing. Both wings were with an endplate and, therefore, widely used as a comparative case for both numerical and experimental simulations [19][20][21][22][23][24][25]. Zerihan investigated the impacts of the ground effect and the angle of attack of the wing on the aerodynamic performance and the wake. One part of his work was also to provide numerical data, for the single element wing, using the CFL3D code (CFL3D Version 6 Home Page, https://cfl3d.larc.nasa.gov/).
In terms of numerical simulations, many RANS simulations were carried out with the geometry of Zerihan [19]. Doig et al. [26] in 2011 performed numerical simulation on this geometry to investigate the compressibility in ground effect and compared various turbulence models (k − ω SST, k − , k − realizable). Vogt et al. performed numerical simulations on the Zerihan's geometry in 2007 [27] and a NACA 4412 profile in 2010 [23,28]. They investigated the lift and downforce generation for such wing accordingly to the ground effect phenomenon. He [29] in 2014 performed numerical simulations on the NACA 4412 profile again with a multi-objective genetic algorithm. His work was based on the shape optimization of the wing with and without ground effect, highlighting the difference for the optimized profiles. Mahon et al. in 2005 [30] performed numerical simulations on the Zerihan's geometry. They compared the wake characteristics, sectional forces, and surface pressures to the experimental data for several turbulence model. Therefore, they highlighted the best turbulence model for each flow characteristic. Furthermore, Grabis et al. in 2018 [31] performed numerical simulations on several wing profiles (CH10, E423, FX74, LA5055, snd S1210) on both single-and two-elements wings. They also compared each profile at different angles of attack to highlight the best configuration for a Formula SAE (FSAEOnline.com, https://www.fsaeonline.com/) front wing.
This research is aimed to investigate the impact of tubercles on the leading-edge of a multi-element wing in ground effect. The geometry used is the two-element front wing geometry used by Zerihan. Furthermore, this geometry was optimized through the sensitivity study with the objective to maximize the L/D ratio. The next section focuses on the computational framework employed in this investigation along with the explanation of the meshes developed and flow conditions. Lastly, the results obtained in this research were analyzed in line with the physics of the flow over the geometry to gain detailed understanding.

Computational Domain
For numerical simulations on a Formula 1 front wing geometry, the domain can be characterized by four lengths: L u for the upstream length; L d for the downstream length; L c for the crosswise length; and H for the height. Usually, to carry out numerical simulations on a wing, a C-type in 2D (or a bullet-type in 3D) mesh would be used. However, for Formula 1 front wing, and more generally for race car front wing, the ground effect has a big impact on the wing performance and the flow behavior. This is why an H-type mesh is chosen for this particular problem. This difference between wings with and without ground effect is also noticeable in the experiments as, for race car front wing, the geometry is placed above a moving ground that has the same speed as the free-stream flow. An example domain is shown in Figure 2. All those lengths are recapitulated in Table 1. Big upstream and downstream lengths have been shown to give more accurate results for the drag and lift coefficients, as the wake is better solved. Vogt et al. [23,27] showed a boundary independence where they extended their boundaries by 5c, i.e., setting L u and H equal to 20c and L d equal to 25c; however, the improvement in the aerodynamic performance turned out to be rather small, i.e., only 0.034% increase on C l and 0.051% on C d .

Validation Case
Many cases have been designed for single-element wings in ground effect but not that many for multi-element wings. In that sense, the work of Zerihan accomplished during his PhD thesis was very relevant, providing both experimental and numerical data for wings in ground effect [19]. As mentioned above, his study was carried out on the NASA GA(W) airfoil, type LS(1)-0413 MOD with slight modifications. It has a constant chord of 223.4 mm for the main element and 165.7 mm for the flap. This gives a combined chord of 380 mm. All length scales were normalized by this chord. The span of the front wing was 1100 mm and the endplates used were rectangular of dimensions 400 × 170 × 4 mm. The "true" incidence of the wing is 14.1 • . The position of the flap was optimized to produce the maximum downforce at a constant flap deflection. For the geometry used here, this optimization was done at a height of h/c = 0.263. Finally, the best flap location was found to be a gap of 0.032c and an overlap of 0.024c, corresponding to, respectively, 12 and 9 mm. This geometry was used in a wind tunnel at the University of Southampton. The freestream velocity was 30 m/s and the freestream turbulence was 0.2%. A moving ground at the speed of the freestream velocity was also used. All experimental simulations on the two-element wing were done for Reynolds number in the range 0.735 × 10 6 -0.765 × 10 6 . Here, for numerical simulations, only half of the geometry was used. That corresponds to a span of 550 mm. The geometry used for numerical simulation is shown in Figure 3. All following numerical simulations were done at a ground clearance of h/c = 0.211 and with the initial conditions of Zerihan's experiments for results replication.

Meshing Strategy
Based on the literature, a numerical domain was established, as detailed in Table 1. This domain has 10c in width, 9c in height, 8c upstream of the geometry, and 15c downstream the geometry. Primary numerical simulations demonstrated the size of the wake induced by the front wing and thus a refinement box was used, which was about the width and height of the front wing and 10c in length. The final domain used in numerical simulations is shown in Figure 4. In this investigation, two different y + were used for the wing and the endplate. In addition to this, the boundary layer of the ground was also meshed for wings in ground effect, as identified in the literature [23,26,27,[30][31][32]34,35]. To compute the initial spacing for the wing and the endplate, a primary study had to be conducted. First, using the Reynolds number gives the range of the freestream dynamic viscosity: This range is found to be 1.8255 × 10 −5 ≤ µ ≤ 1.9 × 10 −5 Pa·s. Choosing µ = 1.8375 × 10 −5 Pa·s, which corresponds to the dynamic viscosity of the air at T = 25 • , leads to ν = 1.5 × 10 −5 m 2 /s. Finally, with all those parameters, the Reynolds number for numerical simulations was Re = 0.76 × 10 6 .
Using the y + calculator (Y+ Calculator-Compute Wall Spacing for CFD, https://www.pointwise. com/yplus/) gives then all wall spacings ∆s. Then, the expansion factor λ can be computed with the wanted number of layers. To do so, the boundary layer δ height must be computed first. An analogy with the Blasius relation for turbulent flow through ducts can be done [36,37]: The boundary layer height δ is then found to be δ = 9.67 × 10 −3 . Then, to compute the expansion factor λ, all n layers height are summed: By solving this equation for a given number of layers, the expansion factor can be found. Concerning the meshing of the overall domain, it can be divided into three regions: the Farfield region, the Refinement region, and the Closefield region. The Closefield region has the geometry in it and is a hybrid region, with structured mesh for the boundary layers of the wing and the endplate. The wing has a structured surface mesh, as well as most of the endplate. The exception is for the surface at the junction of the wing and the endplate, where the surface mesh is hybrid. The Refinement region is fully structured while the Farfield region is fully unstructured. All regions have a structured mesh for the boundary layer of the ground. The final domain can be visualized in Figure 5. For each level of meshes, the numbers are summarized in Table 2.

Fluent Simulation
To perform numerical simulations, a commercial software was used, FLUENT 19R2 (ANSYS Fluent Software|CFD Simulation, https://www.ansys.com/products/fluids/ansys-fluent). As the used freestream velocity is u = 30 m/s, the Mach number is Ma = 0.088. Therefore, the flow is incompressible. A pressure-based, coupled solver was used to get a steady-state solution. Incompressible RANS simulations with the k − ω SST turbulence model were performed using the SIMPLE algorithm with second-order discretization schemes. In total, 15,000 iterations were done for each numerical simulation.
For boundary conditions, a velocity inlet was used as Inlet and a pressure outlet was used as Outlet. For the sides and top, a symmetry condition was used. A full moving ground was used for the ground and a wall condition for the wing and the endplate. For the ground, a semi moving ground was tested (with a moving ground starting only 1c upstream the geometry). It provided less accurate results on the drag and the lift coefficient (of about 0.5% for both) and for the pressure coefficient (of about 1%).

Investigation of Tubercles on the Front Wing
The geometry used for the validation case was also used for the implementation of the tubercles, i.e., two elements with a constant chord and no taper and sweep. For meshing issues, the endplate was removed from the geometry. In the literature, it has been proved that small and widely spaced tubercles are the best for wings without taper and sweep [11]. Here, two elements are used and can change those observations, as the ground effect. Then, the tubercles are efficient only when stall occurs. To study how the tubercles can improve the aerodynamic performance of this wing, stall conditions are mandatory. To trigger such conditions, two options are available: (1) increase the angle of attack of the wing for the given freestream velocity (i.e., 30 m/s); and (2) decrease the freestream velocity for the given angle of attack (i.e., 14.1 • ). The retained option was to decrease the freestream velocity, to let the initial geometry as unchanged as possible. Several simulations were run with, for quantity of interest, the lift coefficient. In fact, this coefficient undergoes a high decrease at stall conditions ( Figure 6a). Finally, a freestream velocity of 0.5 m/s was used for all numerical simulations with tubercles. To optimize the geometry with tubercles, a sensitivity study was done in two steps. The first step, was with four different configurations: three geometries with tubercles on the main element (referred as ME Tubercles geometry) and a geometry with tubercles on the flap (referred as F Tubercles geometry). A change in amplitude and number of tubercles was done. Those geometries are shown in Figure 7. At the end of the first step, the best configurations, i.e., the geometries giving the best L/D ratio, were kept and combined. New configurations were then tested in the second step, in order to get the best geometry with tubercles.

Validation Case
Prior to assessing the validation of the numerical method, a grid convergence study was done. Three levels of mesh were used. Their descriptions can be found in Table 2. To investigate the grid sensitivity, the drag coefficient C D and the lift coefficient C L were chosen. The values for the experimental and numerical simulations are shown in Table 3. Numerical values are very close from each other; nonetheless, the drag and lift are overpredicted. The difference between the fine mesh values and the theoretical values are −1.54% for the drag coefficient and 8.54% for the lift coefficient. To assess on the grid dependence, the Richardson's extrapolation was used [38,39]. To do so, the grid convergence index (GCI) was computed between the coarse and medium meshes and between the medium and fine meshes. Then, a ratio between those GCI was calculated to assess on the spatial convergence of the results. Finally, this ratio is equal to 0.998 for the drag coefficient and 1.004 for the lift coefficient. As these ratios are close to 1, the obtained results can be said to be spatially convergent. For computational time concerns, the medium mesh was used for all the following simulations with tubercles.
Then, the aerodynamic coefficients were compared to the experimental data of Zerihan [19]. It was found that the error was of −1.54% for the drag coefficient and 8.54% for the lift coefficient. The results can, nevertheless, be investigated more in depth, to validate the method and the mesh used. To do so, the pressure coefficient and the wake were investigated and compared to the Zerihan's data. The pressure coefficient at the symmetry plane, i.e., at the middle of the front wing in the span direction, is shown in Figure 8a. The wake survey at x/c = 1.066 at the symmetry plane is shown in Figure 8b. The results of all three meshes were plotted together, to show again the mesh convergence of those results. For the pressure coefficient, the experimental solution and the numerical solution have the same behaviors. On the extrados, the numerical solution is nearly matching the experimental solution over the entire main element and flap. However, on the intrados, the pressure coefficient is mostly over-predicted. On the main element, the spike (the first big pressure drop) is under-predicted by about −7.25% and the peak (the second pressure drop) is over-predicted by about 9%. Over the entire main element, the over-prediction remains around 10%. The same observations can be made for the flap. Those differences do not come from the mesh as it is visible that the pressure coefficient is spatially convergent, and the surface mesh was refined for all three levels of mesh. Those surface mesh characteristics can be found in Table 4. Those differences are coming from the limitations of the k − ω SST turbulence model itself. Looking at the wake, here again a clear under-prediction of about 10% for the x-velocity component is noticeable over the whole height. For the two minima measured by Zerihan, there are both an under-prediction of 10% and a upshift of about 0.02c. Here again, the general behavior of the flow is well captured and the differences come from the turbulence model used. All values for those main points are shown in Table 5. Finally, the flow speed under the front-wing can be stated as being under-predicted, causing both an under-prediction on the intrados of the front wing for the pressure coefficient and on the x-velocity profile. This results in the difference between the experimental and numerical aerodynamic coefficients as well. Nevertheless, overall, the results capture the correct behavior, for both the pressure coefficient and the wake profile at the symmetry plane. The measured errors are acceptable and can be explained. The model was thus verified to proceed an implementation of tubercles.

Investigation of Tubercles on the Front Wing
As explained in Section 2.5, prior to finding the best configuration for a front wing with tubercles, a preliminary sensitivity study was done with four different configurations. As a reminder, ME Tubercles geometry (10 tubercles/2A = 0.05c), F Tubercles geometry (10 tubercles/2A = 0.05c), ME Tubercles geometry (10 tubercles/2A = 0.1c), ME Tubercles geometry (5 tubercles/2A = 0.05c) were built. The freestream velocity was 0.5 m/s to guarantee stall conditions for the front wing.
Looking at the aerodynamic performance in Figure 9a and Table 6, i.e., the lift and drag coefficients, it is striking that the tubercles bring benefits to the front wing. The ME Tubercles geometry (10 tubercles/2A = 0.05c) is aerodynamically better than the F Tubercles geometry (10 tubercles/2A = 0.05c), with, respectively, improvements of 6.9% and 3.6% for the lift. Even if the flap is producing most of the reverse flow, the tubercles on the main element have an effect on the downstream flow. This is why the tubercles on the main element have a bigger impact on the downforce generation. Besides, the increase of drag is more important for the F Tubercles geometry (10 tubercles/2A = 0.05c) than for the ME Tubercles one (10 tubercles/2A = 0.05c), with, respectively, increases of 3.6% and 2.3%. This is understandable as the flap has a bigger inclination. The tubercles thus have a bigger projected area, causing a bigger increase in drag. Having a bigger amplitude than the base one also improves the results substantially. Between the ME Tubercles (10 tubercles/2A = 0.05c) and ME Tubercles (10 tubercles/2A = 0.1c) geometries, an improvement of 3.3% on the lift coefficient is noted. Here again, this comes at drag costs, as the drag coefficient undergoes an increase of 3%. However, having fewer tubercles seems to mitigate the benefits. Between the ME Tubercles (10 tubercles/2A = 0.05c) and ME Tubercles (5 tubercles/2A = 0.05c), a deterioration of 4.8% on the lift coefficient is highlighted. Thus far, the tubercles generate more downforce but at the cost of an increase of the drag generated.   Those aerodynamic coefficients can be explained with the negative x-velocity, characterizing the recirculation zones over the front wing, the Q-criterion, and the shear stress lines over the front wing. Those quantities are shown in Figures 10-13. Looking at the baseline geometry, a big recirculation can be highlighted on the main element, and an even bigger one on the flap. This recirculation is on the whole front wing; however, removing the endplate from the original geometry seems to have been beneficial as no circulation is present near the extremity of the front wing. This can be explained as the lack of endplate promoting the vortices creation, which carries a lot of momentum in the boundary layer and thus delays the stall near the extremity of the geometry. That is why the effect of the tubercles is more visible near the symmetry plane. Besides, the vortices, created by the lack of endplate, bend the stall to the inside. As a main difference between geometries, it is worth noticing that, for each geometry with tubercles on the main element, the induced vortices in the flow field affected the recirculation occurring on the flap, as shown in Figure 13a. This can explain the big differences between the ME Tubercles (10 tubercles/2A = 0.05c) and F Tubercles (5 tubercles/2A = 0.05c) geometries in terms of aerodynamic coefficients. Another way to compare each geometry is to characterize the overall recirculation. To do so, the relative number of cells with reverse flow was used to characterize the size, while the total average negative x-velocity was used to characterize the magnitude as in the Table 7. Furthermore, the numbers relative to the stall position at the first intertubercular and tubercular sections are shown in Tables 8 and 9. For example, the ME Tubercles geometry (5 tubercles/2A = 0.05c), compared to the baseline geometry, reduces by about 2.9% the size of the overall recirculation, while it almost does not changed the magnitude (only a decrease of 0.2%). However, looking directly at the negative x-velocity contours in Figure 10a, the area where reverse flow occurs is smaller. The fact that the overall recirculation size is not that much smaller shows the bigger amplitude of the induced vortices in intertubercular sections compared to the baseline geometry. Compared to the ME Tubercles (10 tubercles/2A = 0.05c) geometry, the ME Tubercles (10 tubercles/2A = 0.1c) geometry obtains amplified effects, showing that doubling the amplitude is beneficial. On the main element, at the first intertubercular section, the stall is delayed 52.8% more. Again, this implies an earlier triggered stall in the first tubercular section of about 4.3%. Concerning the overall recirculation, its size is reduced by 8.7%, while its magnitude is reduced by 6.8%. However, reducing the number of tubercles mitigates the improvements brought by the tubercles. Compared to the ME Tubercles (10 tubercles/2A = 0.05c) geometry, the ME Tubercles (5 tubercles/2A = 0.05c) geometry obtains a bigger overall recirculation size by about 1.8%. However, the overall recirculation magnitude is lower by about 1.3%, showing the size of the recirculation plays a big role in the aerodynamic coefficient values. Concerning the stall on the main element, at the first tubercular section, it is delayed by 4.3%. It is worth highlighting that the stall on the flap was only impacted for the F Tubercles (10 tubercles/2A = 0.05c) geometry, with a delay of 17.9%.
Looking at the ME Tubercles, (5 tubercles/2A = 0.05c), the effects of tubercles are fully visible. First, in the intertubercular zones (i.e., between tubercles), the recirculation has a bigger amplitude.
The minimum x-velocity in this area is 4.9% more important than for the baseline case. This velocity characterizes the most important vortex of the detached flow on the front wing. This means the vortices induced by the tubercles have a bigger amplitude and are more pronounced than the ones on the baseline geometry. On the main element, the stall is delayed, but only in tubercular zones (i.e., at a tubercle section); nevertheless, in the intertubercular section, this stall is triggered earlier.
Looking at the flow field, vortices are visible with the negative x-velocity contours in Figure 13a and the Q-criterion in Figure 13b. Surface vortices on the front wing are created due to the vortices creation of the tubercles in the flow field. They are visible with the shear stress lines over the front wing ( Figure 14). Those surface vortices characterize the induced vortices in the flow field in the intertubercular zones. They are responsible for the delay of the stall in the tubercular zones. The vortices in the flow field are carrying momentum in the boundary layer of tubercular zones, which delays the stall. This was already observed by Hansen [12] with, in particular, primary vortices created in intertubercular zones (See Figure 11 of [12]). Talking about the pressure coefficient, here again, a distinction between tubercular and intertubercular has to be made (Figure 15). At the tubercular section, the spike is lower than for the baseline case, by about 23.8%. On the contrary, at intertubercular sections, this spike is higher, by about 30%. It is worth mentioning that between a tubercular and an intertubercular section, the pressure coefficient is very close to the one obtained with the baseline geometry (Table 7). Those differences in the spike pressure coefficient indicate an overpressure at the intertubercular section and an underpressure at the tubercular section. Those pressure drops along the leading edge create surface vortices that induce vortices in the flow field at intertubercular sections. Those last vortices carry momentum in the boundary layer at tubercular sections and delay the stall. This is visible especially with the negative x-velocity contours or the shear stress lines over the front wing. However, the stall is triggered earlier at intertubercular sections.
Moreover, the ME Tubercles (10 tubercles/2A = 0.1c) geometry faces a big recirculation between Tubercles 3 and 4 (Figure 12a). This testifies to the chaotic character of the recirculations created by the tubercles. This phenomenon was already observed by Srinivras [18] (See Figure 25 of [18]). Looking at the shear stress lines in Figure 12c, two secondary surface vortices can be seen, near the trailing edge of the main element. This results from the big recirculation between Tubercles 3 and 4. Now comes the tricky part. Even if increasing the amplitude and the number of tubercles seems to be better, as the overall recirculation is smaller, it is not the case. One other big aspect of the tubercles is the position of the recirculation created. For some geometries, recirculations are near the symmetry plane, improving the performance. Those geometries are favorable configurations. For other geometries, those same recirculations are not near the symmetry plane, altering the results. These are unfavorable configurations.
Finally, it is worth saying that the more upstream the element is, the more benefits the tubercles bring. This was shown by comparing both ME Tubercles (10 tubercles/2A = 0.05c) and F Tubercles (10 tubercles/2A = 0.05c) geometries. This observation seems a little bit counter intuitive at first sight. In fact, the flap is the element with the higher incidence while facing the upstream flow (14.1 • incidence). It is on this element that the recirculation is more important. Nevertheless, the tubercles have more impact on the recirculation of the main element. Talking about their amplitude, bigger tubercles seem to improve the results, as they amplify the creation of big vortices and more delay the stall at intertubercular sections. This amplification can also be seen downstream. Concerning the number of tubercles, the more there are, the better it is. The vortices induced by the tubercles are more numerous, and so they are more beneficial. In addition, those vortices are more pronounced, due to their proximity. If those two last observations are combined, it can be stated that a high amplitude gradient along the span should do the trick. However, a limit exists due to the fact that the recirculation induced by the tubercles is chaotic in both their positioning and behavior. In that sense, using only half of the geometry is not sufficient. As the 3D effects are part of the tubercles, and because they are unpredictable, conducting numerical simulations with a full geometry should be considered. Comparing those observations to the literature, a big difference is highlighted. As a reminder, for wings with no sweep and taper, the best configuration for tubercles is the one with the smallest amplitude and the biggest frequency, as shown by Hansen [11]. His study was done on a single-element wing, at different angles of attack. Here, the best configuration is the one with bigger amplitude and more tubercles, and is was achieved for a given angle of attack. The multi-element wing and the ground effect affect the tubercles behavior. However, the observations should be temperate. In fact, as shown by this study, optimal configurations for both amplitude and number of tubercles exist. In other words, limits exist for the amplitude and the number of tubercles, from which an alteration in aerodynamic performance is highlighted. The tubercles are thus case-dependent. Now that this was stated, the optimal configuration with tubercles could be found. To do so, all benefits obtained with those geometries were combined. Four more geometries with tubercles on the main element and on the flap (referred as MEF Tubercles) were built. These new geometries are shown in Figure 16.       Here again, those results can be explained with the negative x-velocity, Q-criterion, and shear stress lines over the front wing. First, between the baseline and the MEF Tubercles (10 tubercles/2A = 0.1c), high differences can be highlighted. Those differences are even more pronounced than the ones of the preliminary sensitivity study. A big part of the flow is now keep attached on the front wing. There is a big recirculation occurring on the main element bending all other recirculations. Concerning the stall, on the main element, at the first intertubercular section, it is delayed by 230.1% (almost to the trailing edge), and, at first tubercular section, it is triggered earlier by 24.5%. The overall recirculation has its size reduced by 31.8% and its magnitude reduced by 28.4%. As stated with the aerodynamic coefficients, doubling the amplitude greatly affects the flow. Between the MEF Tubercles (10 tubercles/2A = 0.1c) and the MEF Tubercles (10 tubercles/2A = 0.2c), on the main element, the stall is delayed by 5.1% at the first intertubercular section and triggered earlier by 28.6% (almost to the leading edge). Besides, the overall recirculation has its size reduced by 11.5% and its magnitude reduced by 5.2%. Even though those numbers seem to be beneficial, the aerodynamic coefficients are worst for the MEF Tubercles geometry (10 tubercles/2A = 0.2c). The chaotic behavior of the recirculations creation plays a big role in the aerodynamic coefficient. For the MEF Tubercles (10 tubercles/2A = 0.1c), one of the big recirculations is cut by the symmetry plane, while, for the MEF Tubercles geometry (10 tubercles/2A = 0.2c), none of the big recirculations are cut by the symmetry plane. For the whole geometry, the computation of the surface pressure and the aerodynamic coefficients are altered. Increasing the number of tubercles leads to the same conclusions. While a net improvement is shown for the MEF Tubercles geometry (16 tubercles/2A = 0.1c), the results are altered for the MEF Tubercles geometry (20 tubercles/2A = 0.1c). Between the MEF Tubercles (10 tubercles/2A = 0.1c) and MEF Tubercles (16 tubercles/2A = 0.1c) geometries, no more big recirculations are on the main element. The two big recirculations remain on the flap with one of them being cut in half by the symmetry plane. Despite the stall, on the main element, at the first intertubercular section, being triggered earlier by 19%, and, at the first tubercular section, being triggered earlier by 15.1%, the MEF Tubercles geometry (16 tubercles/2A = 0.1c) is more aerodynamic. Besides, the overall recirculation size is slightly decreased while its magnitude is increased by 13.9%. Here again, the fact one recirculation is at the symmetry plane plays a big role on the surface forces. The cancellation of the recirculation on the main element is also influential. Nevertheless, the comparison is simpler for the comparison between the MEF Tubercles (16 tubercles/2A = 0.1c) and the MEF Tubercles (20 tubercles/2A = 0.2c) geometries. The overall recirculation has its magnitude slightly decreased by 2.3% but its size is increased by 6.8%. Concerning the stall on the main element, at the first intertubercular section, it is triggered earlier by 14.5%. This second set of geometries leads to more conclusions. First, using tubercles on both the main element and the flap substantially improves the aerodynamic performance of the front wing. The recirculation over the front wing is mostly removed and only two big recirculations remains. One of them, on the main element, is removed by using 16 tubercles instead of 10. On the main element, the stall is importantly delayed. On the flap, it depends of recirculations. If there is no big recirculation, the stall is delayed to the trailing edge. If there is a big recirculation, the stall is triggered at the leading edge. Finally, the best configuration thus far was the MEF Tubercles geometry (16 tubercles/2A = 0.1c). The lift is improved by 22.6%, the drag is increased by 12.1%, and the L/D ratio is improved by 9.4%. The overall recirculation has its size decreased by 31.9% and its magnitude decreased by 16.8%. On the main element, at the first intertubercular section, the stall is delayed by 167.7%, but, at the first tubercular section, the stall is triggered earlier by 35.9%. On the flap, the stall is delayed to the trailing edge in regions where no big recirculations occurs. It is this geometry that will be used later to be optimized by the adjoint solver. The comparison between the baseline geometry and the geometry with the final configuration of tubercles can be found in Figure 17. Table 7. Comparison of the peak pressure coefficient on the intrados for both baseline and ME Tubercles (5 tubercles/2A = 0.05c) geometries at different section along the span of the front wing.
x/c C p

Concluding Remarks
This paper has as a primary aim to study a passive turbulence flow control, tubercles of the humpback whale flipper, on a F1 front-wing geometry, very close to a Tyrrell wing, in ground effect. To do so, a validation case for multi-element wings was defined and conducted. Based on the Zerihan's work [19], k − ω SST simulations were validated. Then, the tubercles were implemented on the leading edge of the front wing geometry. A sensitivity study was done to understand the behavior of the tubercles and find the best configuration for this geometry. The following points are the key conclusions from this paper:

•
In stall conditions, for all geometries tested, tubercles improved the aerodynamic performance of the front wing, i.e., increased the downforce generation but at the cost of an increase of the drag. Nevertheless, the L/D ratio was increased for each geometry. • By sets of underpressure and overpressure occurring, respectively, at tubercular and intertubercular sections, vortices are created in the flow field, at intertubercular sections. Those vortices carry momentum in the boundary layer that delays the stall in intertubercular sections, at the cost of triggering earlier the stall in tubercular sections. • The more upstream the tubercles are, the more impact on the flow and aerodynamic coefficients they have. • One important aspect of the tubercles is that they alter the creation of the recirculations on the wing that improve the aerodynamic performance. Nevertheless, they can lead to the creation of big recirculations that are chaotic by their positions. For each geometry, the big recirculations highly impact the results.

•
Compared to what was stated by Hansen K. L. [11], it can be said that an optimal configuration can be found for any geometry. Even if having more tubercles with a bigger amplitude seems better, as it reduces the size of the overall recirculation on the front wing, it alters the aerodynamic performance of the wing. This shows the tubercles are case-dependent.
The tubercles were implemented here on a basic front wing in stall condition. Of course, such a condition is favorable for a tubercle implementation. Nevertheless, it could be interesting to implement tubercles on a more realistic geometry and with more realistic initial conditions. With a more realistic geometry, each element of the front wing will have a non-constant chord, a non-constant incidence, and can be bent in the spanwise direction. Moreover, those elements have a higher incidence. In such a case, the legitimacy of implementing tubercles on Formula 1 front wing could be proved. However, more simulations and more parameters should be taken into account: many freestream velocities, many ground clearances, many pitch angles, and many yaw angles. Because, during one lap, all these parameters change, it is worth investigating their effects combined with the tubercles effects to assess the relevance of such a passive turbulence flow control.
Furthermore, a bio-mimicked distribution can be used and compared directly to the sinusoidal distribution. Such a distribution can be found in the the works of Fish F. E. [7] and Srinivras K. S. [18]. The bio-mimicked distribution, closer to the tubercle distribution on the humpback whale flipper, could be aerodynamically more efficient than a sinusoidal distribution. In fact, due to the species' evolution, the humpback whale flipper has done its own CFD on its flippers. This distribution can be legitimately said to be more efficient. However, for the same reasons, this particular distribution is fitted perfectly in a range of conditions (freestream speed, angle of attack, etc). As a front wing uses very different conditions, this distribution cannot be that efficient. Nevertheless, one thing that should be kept from such a distribution is the biggest tubercles. For a bio-mimicked distribution, the Tubercles 1 and 4 are the most prominent. They have a higher impact on the flow field, as highlighted in Figure 17. Srinivras [18] showed that a geometry with only two tubercles was almost equal to the bio-mimicked distribution, in terms of aerodynamic performance. As shown in Section 3.2, only few tubercles have the same behaviors. Testing a distribution with only few-but pronounced-tubercles could also be relevant. This distribution would also be easier to manufacture. Using a different distribution for each element could also be done, which might improve the results significantly.  Acknowledgments: All the simulations were performed using Cranfield High Performance Computing Cranfield HPC) facility. We acknowledge the support received from Mick Knaggs from Cranfield HPC during the simulation process.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: