2.1 Boundary conditions

To complete the specification of the problem, boundary conditions must be defined over the entire boundary of the computational domain. In the case of simulating the full BLOODHOUND SSC vehicle aerodynamics, this includes farfield boundaries, viscous walls, jet engine/rocket exhausts and jet intake^{(Reference Evans and Rose11)}.

2.1.2 Viscous walls

For viscous flow, the no-slip condition:

(1) $$\begin{equation} \vec{u} = \vec{u}^w \end{equation}$$

is strongly applied, where $\vec{u}^w$ is a specified wall velocity. When simulating a vehicle with rotating wheels, the vector $\vec{u}^w$, at a point with position vector $\vec{r}$ on the wheel relative to a user-specified position on the axis of rotation, is computed as:

(2) $$\begin{equation} \vec{u}^{\,w} = \vec{\omega }\times \vec{d}, \end{equation}$$

where

(3) $$\begin{equation} \vec{d} = \vec{r} - \frac{(\vec{\omega }.\vec{r})}{\left|\vec{\omega }\right|}{\vec{\omega }} \end{equation}$$

is the position vector of the point relative to the centre of the wheel and $\vec{\omega }$ is the user-specified wheel rotation vector. All the simulations are performed with the vehicle at rest, and this requires that the ground moves with the speed of the vehicle but in the opposite direction. This boundary condition definition for the ground plane assumes that it behaves as an impervious surface at the same temperature as the local flow. The validity of this assumption will be explored during early vehicle testing. Note that the wheel intersects with the rolling ground surface, assuming a penetration depth of 15 mm. This assumption is based on observations of wheel–ground interaction behaviour during the Thrust SSC LSR project^{(Reference Noble1)}. This is another assumption that will require validation during vehicle testing. Figure 6 shows a slice through one of the rear wheels and fairing. Note that the flow is fully simulated in the space between the wheel and fairing.

Figure 6.

A longitudinal (constant y) mesh cut through the wheel centreline indicating the resolution and boundary layer mesh. Note that the flow is being modelling in the gap between the rotating wheel and static fairing.

2.2 Improving computational performance

To ensure a reasonable turnaround time for each simulation, the computational performance of the solution algorithm is improved by the use of parallel processing. The parallel implementation implemented in FLITE3D is based on physical domain decomposition. In order to optimise performance and achieve scalability on a large number of processors, the chosen domain decomposition strategy must produce sub-domains of a balanced size and with a minimum number of cut edges. This ensures that each processor has to perform an equal amount of work and that the amount of communication between processors is minimised. This has been achieved by utilising the METIS family of partitioning algorithms^{(Reference Karypis and Kumar27)}.

During this study, FLITE3D’s solver and ancillary functions were executed using the computational power of Swansea University’s Zienkiewicz Centre for Computational Engineering parallel computer cluster. Each CFD study was conducted using a cluster of 36 Opteron 240 dual-core 1.4 GHz processors. Each individual steady state CFD simulation took approximately 6 days of wall clock time on meshes of approximately 25 m elements.

2.3 Aerodynamic coefficients

The FLITE3D system post-processor integrates pressure and skin friction loading across the vehicle surface and outputs aerodynamic forces as coefficients normalised by the dynamic pressure as:

(4) $$\begin{equation} C_L = L/q_{\infty }, \\ \end{equation}$$

(5) $$\begin{equation} C_D = D/q_{\infty }, \\ \end{equation}$$

(6) $$\begin{equation} C_T = T/q_{\infty }, \end{equation}$$

3.1 Creation of the parametric geometry

In order to reduce the full vehicle modelling complexity, simplifications were introduced: the intake was blocked, jet and rocket exhausts were omitted, and the winglets, vertical fin and front wheels were removed. The vehicle was also halved in the longitudinal axis, which would later be accounted for in the FLITE3D system with the use of a symmetry plane to mirror the model. These simplifications are illustrated in Fig. 7, where the half-car parametric model is presented.

Figure 7.

An annotated diagram of the simplified parametric geometry with the parametric fairing.

The fairing itself needed to be defined by several constrained design parameters that could be varied to control the resulting geometry for the DoE optimisation process. To achieve this, the fairing was split into three components: the leading spike, the trailing spike, and the Static Central Hub (SCH). The SCH made up the central part of the fairing to which the leading and trailing spikes were anchored. The SCH was also attached to the suspension struts connecting the wheel to the main body of the vehicle, as well as being the point of attachment for the wheel axle. Due to the degree to which the SCH hugged the profile of the wheel (around 16 mm of clearance between the wheel and the fairing in the parametric model), opportunity to vary its geometry was limited; this led to the SCH remaining static throughout all of the DoE optimisation permutations.

The scope for variation in the fairing geometry was therefore limited to adjusting the positions of the leading and trailing spike tips in three-dimensional space. To enable both leading and trailing tips to be adjustable in 3 degrees of freedom but also independently of one another to allow for asymmetric configurations three parameters were required for each tip. This gave the parameters depicted in Fig. 8. The leading and trailing spikes are pyramids, with their bases attached to the rectangular faces on the SCH and tips defined by the variable design parameters.

Figure 8.

The variable dimensions (ϕ 1-6) are shown. The variable dimensions define the position of the leading and trailing tip, which in turn define the shape of the leading and trailing spikes. The SCH forms the central part of the fairing and remains static throughout.

3.2 Validation of the parametric model

The parametric model needed to closely mimic the aerodynamic behaviour of the full-fidelity model to be suitable for basing design decisions on, particularly in the region of the fairing. To assess the validity of the parametric model, its performance in the flow was compared to the full-fidelity model with the aforementioned baseline wheel fairing attached. This baseline fairing model featured smooth, continuous surfaces and far more accurate modelling of the struts and connections between the fairing and the fuselage. The parametric model, in contrast, features rough simplifications and sharp edges. These differences are illustrated in Fig. 9.

Figure 9.

Comparison of the parametric fairing geometry (left) with the full-fidelity fairing geometry (right).

After conducting a CFD study on both models at a Mach value of 1.3, the lift and drag coefficients (C_L and C_D) were compared with the full-fidelity geometry CFD predictions. The differences in lift and drag characteristics between the two models are shown in Table 1.

Table 1

Comparison of full-fidelity and parametric model CFD force coefficient predictions for baseline model at M_∞ = 1.3

It is clear from Table 1 that the there is a discrepancy between the predicted C_L and C_D of the parametric and full-fidelity models. However, when the aerodynamic characteristics of the two fairings were compared, it was found that the parametric model was the characteristics of the flow phenomena predicted were similar. This is illustrated by the flow visualisation depicted in Fig. 10. A plane parallel to the ground raised to a height approximately equal to the vertical midpoint of each fairing was cut through the flow field. The coefficient of pressure was then mapped across it. Both models predicted a detached bow shock of similar strength and shock angle, secondary shock inboard of the fairing close to the ‘throat’ and downstream inboard and external expansion fans (albeit stronger and larger in extent on the parametric model). Due to the similarities in key flow features around the fairings, it was deemed reasonable to assume that optimising the parametric model would result in improvements that would also translate back to a full-fidelity model.

Figure 10.

A cut was made parallel to the ground plane through (a) the parametric fairing, and (b) the full-fidelity fairing, to show the similarities in the pressure coefficient field surrounding them at M _∞ = 1.3. The cut planes were made at the same height from the ground, and as a result the parametric fairing appears to not come to a point on each end. Note that neither the jet nor rocket efflux is simulated in the parametric model.

5.1 Selection of the design parameters and constraints

In any optimisation problem, the optimum design parameters are those which minimise the objective function under consideration. In the case of this study, the design parameters were the dimensions $\phi _1\text{ to }\phi _6$ that were defined earlier in the paper (Section 3.1). A set of these design parameters defines a full, unique fairing geometry that could be investigated with CFD simulations.

The constraints of the design variables were dictated by multi-disciplinary engineering considerations. For example, ϕ₂ was limited to prevent the trailing point from extending into the region of the jet efflux at the rear of the vehicle; likewise, ϕ₁ was limited to avoid over-extending the leading point to the extent that lack of rigidity would translate to undesirable structural vibration or aeroelastic effects. The full set of limits for the dimensions is detailed in Table 3 with reference to Fig. 8.

Table 3

Constraints placed on the parametric design space

5.2 Objective function

The primary focus of this study was to minimise the aerodynamic drag caused by the rear wheel fairing. The most direct approach to achieve this was to optimise considering only the total integrated aerodynamic drag. An objective function based solely on aerodynamic drag is represented as,

(7) $$\begin{equation} f(\phi _1,\phi _2,\phi _3,\phi _4,\phi _5,\phi _6) = C_D = D/q_{\infty } \end{equation}$$

More involved approaches, however, were conceivable; objective functions that accounted for the lift and lateral aerodynamic forces produced by the fairing would give even more insight into how changes in the geometry affect the aerodynamic characteristics of the fairing. This would allow for even greater tailoring of the fairing to the needs of BLOODHOUND SSC. For example,

(8) $$\begin{equation} \begin{array}{ c c l}f(\phi _1,\phi _2,\phi _3,\phi _4,\phi _5,\phi _6) &=& C_D + \alpha C_L \\ & = & D/q_{\infty } + \alpha L/q_{\infty } \end{array} \end{equation}$$

Prior to this study, BLOODHOUND underwent optimisation based on CFD experiments conducted with the baseline wheel fairing^(6-8). As a result, it was decided that it would be valuable to measure how greatly the lift force incurred by the optimised parametric fairing deviated from the lift incurred by the baseline fairing. An example of an objective function that would represent this is,

(9) $$\begin{equation} f(\phi _1,\phi _2,\phi _3,\phi _4,\phi _5,\phi _6) = C_D + \alpha \left| C_L - {C_L}_{baseline} \right|, \end{equation}$$

where α again represents an arbitrary weighting factor, and the normalised lift term with the subscript of ‘baseline’ denotes the total integrated lift produced by the vehicle model with the baseline fairing configuration.

It was concluded that multiple objective functions would be investigated: one based purely on drag, as described in Equation (7), and one that considered both drag and lift, as described in Equation (9).

5.3 Sampling of the design space

The design space is the six-dimensional Cartesian space defined by the six design parameters and their associated constraints. This design space was sampled by conducting CFD experiments on fairings defined by sets of parameter values.

Based on previous experience with simulations conducted on similar-sized meshes with FLITE3D, a goal of conducting 60 CFD simulations was set, thereby sampling 60 points in the design space. Based on this restriction, several experiment designs were considered and rejected. For example, a classical factorial experiment design would require 3⁶ = 729 sampling points for an experiment with six design parameters^{(Reference Narayanan and Iyengar35)}. Ultimately, the Latin Hypercube Sampling (LHS) technique first described in 1979 by McKay^{(Reference McKay, Beckman and Conover36)} was employed. The LHS technique is capable of generating sets of parameter values for multi-dimensional design spaces that are more effectively spread than if they were selected at random^{(Reference Stein39)}.

The sets of parameter values, normalised between zero and one, are shown in the Appendix. The resulting total integrated aerodynamic force coefficients are given for each fairing. Note that some of the design space positions resulted in solutions for which a sufficiently converged steady-state solution was not possible. These positions were therefore left out of the sampling used to construct the DoE response surface.

5.4 Response surface modelling

Response Surface Modelling (RSM) is a statistical technique first described by Box and Wilson^{(Reference Box and Wilson37)} to define relationships between multiple input variables and one or more response variable. Response surfaces were fitted to the data collected from the CFD experiments using the rstool in the MATLAB statistics toolbox⁽³⁸⁾. The resulting Root Mean Square Errors (RMSEs) of the C_D fit for each of the available models is presented in Table 4 along with an indication of the number of free parameters used in each model.

Table 4

Comparison of the errors in the response surface models for f = C_D

The linear model was selected as the basis for the optimisation since, despite having a slightly higher RMSE, the predicted error at the point of minimum drag was lower than in the pure quadratic model. The generated linear response surface model function was

(10) $$\begin{equation} \begin{array}{c c l}C_D & = & f(\phi _1,\phi _2,\phi _3,\phi _4,\phi _5,\phi _6) \\ & = & 1.0009-0.0334\phi _1-0.0102\phi _2 \\ & &+\, 0.0024\phi _3-0.0164\phi _4 \\ & &+\, 0.0254\phi _5+0.0217\phi _6 \end{array} \end{equation}$$

Figure 13 shows predictor plots for each of the six design space parameters centred on the middle of the design space. The red lines indicate upper and lower 95% confidence bands. Note that gradients of these plots give a measure of how important each of the six parameters is in determining the supersonic fairing drag. It suggests that ϕ₁, the leading spike length, ϕ₃ the height of the leading spike tip, and ϕ₄, the height of the trailing spike tip from the ground, are the most important in controlling the fairing drag, whereas drag is insensitive to ϕ₂ the trailing spike length.

Figure 13.

Predictor plots for each design parameter centred on the middle of the design space showing predicted D/q against normalised parameter position. Red lines indicate upper and lower 95% confidence bands.

Despite these useful insights, it was decided, that the errors in the response surface were too large, particularly at the edges of the design space. Linear response surfaces will always have a minimum at the edge of the design space. An alternative approach for interpreting the sampled data that would require minimal CFD runtime was sought.

6.1 RBF interpolation

A Radial Basis Function (RBF) is a real-valued function whose value is dependent solely on the scalar Euclidean distance of the argument from a defined centre such that:

(11) $$\begin{equation} \psi (r)=\psi (\left|x-x_i\right|), \end{equation}$$

where r is the Euclidean distance between x_i and x, which are the position of the centre, and the position for which the interpolated value is to be found, respectively. RBFs are particularly well suited for interpolating in cases of high dimensionality with scattered sampling points^{(Reference Buhmann40)}. This makes RBFs particularly suitable for this study, which had generated data at sampling points determined by LHS across a six-dimensional design space.

There are a range of RBF types and the following are popular choices:

• • Gaussian: ψ(r) = e ^{−θr 2}

• • Polyharmonic spline: ψ(r) = r^k, k = 1, 3, 5. . . or ψ(r) = r^kln(r), k = 2, 4, 6. . .

• • Thin Plate Spline: $\psi (r)=r^2 {ln}(r)$

• • Multi-quadratic: $\psi (r)=\sqrt{1+\frac{r^2}{\theta ^2}}$

The positive constant θ that features in the Gaussian and multi-quadric RBFs is the shaping parameter. The optimal value for this parameter is estimated by the average distance between the interpolation nodes. A function approximation is then constructed typically of the form:

(12) $$\begin{equation} f(x)=p(x)+\sum ^{n}_{i=1}\lambda _i\psi (\left|x-x_i\right|), \end{equation}$$

where n is the total number of sampled points, x_i is the unique centre about which each RBF is associated and λ_i is the weighting coefficient. The approximation function is augmented by the polynomial p(x), which is of 1 degree less than the RBF itself to ensure a unique solution for the given weighting vector so that p(x) is typically of the form:

(13) $$\begin{equation} p(x)=c_0+c_1x, \end{equation}$$

where c ₀ and c ₁ are constants that are chosen in conjunction with the values of the weighting vector λ_i so that the approximation of the function matches the value of the function at the nodes.

The previously listed range of RBF function options were tested, and it was determined that the multi-quadratic option was best at interpolating the CFD drag values at the sample points with minimum error. The tests involved systematically removing sampling point data and then interpolating the objective function at this point using all other sample points.

6.2 Nelder–Mead optimisation

MATLAB’s Nelder–Mead optimisation function, in conjunction with RBF interpolation of the known sample point data, was used to seek an optimum fairing geometry. The Nelder–Mead algorithm requires a starting ‘guess’, and since it suffers from the problem of getting trapped in local minima^{(Reference Nelder and Mead8)}, a range of starting points was trialled to establish how sensitive the final ‘optima’ found was to the starting position. It is a well-known property of RBF interpolation that accuracy degradation occurs near boundaries of sampled data. This limited the perceived validity of solutions with ϕ values at extrema^{(Reference Fornberg, Driscoll, Wright and Charles41)}. To counteract this, the searchable sample space was restricted to 0.1 < ϕ_i < 0.9. It was found that the same ‘optimum’ was discovered regardless of start position (implying that the simple C_D objective function, Equation (7), did not contain multiple local minima). The convergence criterion used to stop the algorithm was that convergence was assumed when neither the objective function nor design parameter values changed for 100 iterations. The evolution of the drag coefficient and design parameters as a function of iteration number for a starting point at the centre of the design space (i.e. the baseline case) is shown in Fig. 14.

Figure 14.

Evolution of C_D (top) and design space parameter values (bottom) with Nelder–Mead iteration number for a starting position at the centre of the design space.

The supersonic drag coefficient reduces from the baseline value of 1.04 to 0.979, a reduction of 6%. It is interesting to note that these results are broadly in agreement with the indications given by the simple DoE study. When compared with Fig. 13, the parameters identified as being likely to have a strong influence on drag are indicating the sense in which they influence drag correctly. For example, the DoE indicated that the a large leading spike length, ϕ₁, is likely to have positive influence on drag minimisation and the Nelder–Mead optimisation validates this. Also, the DoE indicated that both low leading- and trailing-edge offset values were likely to result in minimum drag, and this, too, is corroborated in the results of Fig. 14. The increased complexity that is captured in the RBF inter-polation allows the optimiser to find the exact position of the optimum in the design space rather than just being pushed to the design space limits, as was the case with the linear response surface in the DoE study.

Having constructed the geometry resulting from the predicted optimum design space position and run a CFD simulation at this new point in the design space, it was verified that the RBF interpolation of the CFD prediction supersonic drag coefficient was within 1% of the RBF interpolation value. The optimum geometry and Mach 1.3 pressure distribution across it is shown in Fig. 15 compared with the baseline case.

Figure 15.

Comparison of the (a) baseline and (b) minimum drag optimum pressure distribution at Mach 1.3.

Figure 16 shows pressure coefficent plots at constant z (horizontal) cuts through the flow domain in the vicinity of the rear wheel fairing for the baseline fairing geometry and minimum drag optimum. Note that on the ground plane, the optimised fairing has resulted in the focus of compression downstream of the fairing’s bow shock has been shifted in board (towards the car). Above the ground plane, the strength of the bow shock has been reduced, which is one of the contributing factors towards the drag reduction. Also, it is evident that the extension of the rear spike in the optimised configuration gas reduced the strength and extent of both the inboard and outboard expansion fans which will also contribute to the drag reduction.

Figure 16.

Comparison of baseline (a,b,c,d) and optimum (e,f,g,h) pressure coefficient plots at constant z (horizontal) slices through the domain in the vicinity of the rear wheel fairing. Note that z = 0.0 represents the ground plane and that the wheel contact patch can be observed in plots (a) and (e).

The optimum minimum drag geometry had a lift coefficient of 0.172, 50% lower than the baseline value. This is due to a more pronounced high-pressure zone sitting above the ‘delta’ suspension strut, which is evident in Fig. 15. Although this study was focussed on drag minimisation, it was also decided to be beneficial to not adversely affect the lifting behaviour of the rear wheel and suspension system. Therefore, the optimisation was reconsidered with an objective function formulated such that there was a penalty incurred in changing the supersonic lift coefficient from the baseline model lift coefficient (Equation (9)). The alpha parameter in Equation (9) controls the extent to which the design case is penalised for having a supersonic lift coefficient that is different to the baseline case. When alpha is zero, this is not considered and it is simply drag that is being minimised, and as it is increased to unity, equal weighting is given to drag and lift variation. The effect of taking lift variation into account on the optimum design space position is shown in Fig. 17.

Figure 17.

Variation in predicted optimum parameter values as a function of α in the objective function.

It is evident from Fig. 17 that if the requirement to keep the supersonic lift response of the rear fairing geometry similar to the baseline characteristics is strongly enforced, then this requires shifting from the simple minimum drag optimum geometry. In particular, the results in Fig. 17 suggest that in order to minimise drag with minimum influence on lift, then the trailing-edge spike length should be minimised rather than maximised and the trailing spike offset distance increased rather than decreased. This could be because these parameters have a much stronger influence on the lifting characteristics of the geometry than the drag characteristics of the geometry.

Ultimately, it was decided that drag minimisation was more critical than the lifting behaviour (since the vehicle has winglet trimmers to account for lift response variations), so the final design for installation on the BLOODHOUND SSC vehicle was that in optimum shown in Fig. 15.