Certainly not

Verification: the act of waiting for the next simulation result, desperately hoping it's not some crazy value that invalidates the last 3 days of work

This is a bit of a different post that focuses on a much less practical and instead more numerical issue. The concept itself has nothing to do with aerodynamics, but is an important aspect of any performance prediction method or experimental testing. 

I once heard someone say that verification is the process of making sure you're solving the problem right, and validation is the process of making sure you're solving the right problem. This is both the simplest and most accurate description I've heard, so I've adopted this as my own explanation whenever someone asks. What I'm covering in this post is verification of the EV24 CFD setup featuring the as-manufactured aero surfaces; in particular the process of conducting an uncertainty analysis to determine the effective numerical precision/repeatability of our simulation results. We ran verification tests on EV23 back when we first put together the original version of the current CFD setup but that was much less rigorous due to limited time, no correlation data (and thus mostly default physics and boundary/wall settings), very low force magnitudes due to having no aerodynamic components, and consequentially a very simple flow field. It was conducted inefficiently with some poor decisions made, which at the end we "resolved" by slapping a huge safety factor on all our results that completely dominated the actual uncertainty calculations. We then somewhat ignored the existence of any uncertainty for the rest of the year until expressing the final performance estimates for the 2024 aerodynamic floor. That's why I've waiting until now to post about this process.

Throughout 2024, small improvements were constantly being made to the CFD setup to address stability/repeatability issues and to increase mesh efficiency. This resulted in far nicer and more expected results during the first few tests from this new uncertainty study, so we decided we'd put in the effort and do it properly this time out, as it's probably a process that the team will only do rigorously and in its entirety every second year at most. This meant setting up dedicated field probes, reading up on uncertainty procedures and formulas, and creating a program to automate calculations. 

Parameters and monitors

There were two key results we wanted for each group of verification tests. These were the ideal parameter values for mesh and domain sizes, and an uncertainty estimate for each result. Both of these results rely on assessing solution convergence between iterative simulations. The first order of business was to set up new monitors so we had a solid selection of results from which to assess convergence. Lift, drag, and balance were the three results we wanted uncertainty estimates for, so obviously these value were also used to assess convergence. Additional to these were lift-drag ratio, six total pressure monitors, and two static pressure monitors. The total pressure monitors were placed in critical areas that were expected to be influential on total force estimates (either locally or downstream) and/or sensitive to small input changes. The static pressure monitors sampled the maximum positive and negative pressure magnitudes on the floor part. 

Total pressure converge monitor locations:

The input parameters considered for uncertainty estimates were mesh and domain sizes, as mentioned. AMR cut-off values were also assessed, but not rigorously. Initially, uncertainty due to iteration-level oscillations in results within an individual simulation were considered, but the moving average used to provide the final value for each result suppressed these oscillations such that the associated uncertainty was an order of magnitude smaller than that from the mesh and domain size changes, and was therefore deemed negligible. 

General method

To minimise superimposing of convergence or scatter associated with each individual input parameter with one another, it was decided that the "baseline" setup would be conservative and adopt large domain sizes and a small mesh base size. From here, each parameter would be iteratively coarsened until either divergence clearly dominated over any underlying scatter, or the simulation itself became unstable. All tests used a refinement ratio of 1.2, except for changes to mesh size in the skidpan case, which used a ratio of 1.1 to increase available data points without blowing out the mesh count. It was then decided if values finer than those of the baseline setup should be investigated. This process was repeated for both the straight ride-height case and the skidpan case, from which the uncertainty estimates would then also be applied to the braking case and median corner case respectively, as percentages. 

The straight ride-height case was tested first, in which the domain sizes were the first investigated parameter. Domain width and height were tested simultaneously to reduce the total number of simulations required (this may be less suitable if a rear wing is present, given the increased up-wash). Figure 1 shows the table of results for these first tests as an example of how results were recorded.  Figure 2 shows the convergence plots for all domain and mesh size tests for the straight ride-height case. Pressure monitor convergence was assessed only qualitatively, hence their inclusion as "spark lines" in Figure 1 rather than detailed plots. All domain sizes for the straight airflow ride-height case are measured from the origin to the respective boundary. 

Figure 1: Results table for changes in domain width and height between 17.3 m to 4.02 m (cross sectional areas of 300 m2 and 16.2 m2 respectively). In this case, the "baseline" was 10m, but larger sizes were tested to increase confidence in convergence behaviour and uncertainty calculations. The green cells represent the final selected value for that parameter. 

Figure 2: All main convergence plots for the straight ride-height case (the "coarsest" value for both the domain width and the domain aft parameters were excluded from uncertainty calculations based on apparent divergence beyond the asymptotic range - discussed below)

Discussion of results

For the straight airflow ride-height case results in Figure 2, most exhibited good convergence curves with varying degrees of scatter/noise. Lift and drag had the lowest scatter relative to convergence across all parameters, followed by aero balance and then lift/drag ratio. There are a few interesting point to mention:

There is a degree of sanity-checking that can be done for the domain size tests, particularly the width/height and fore length parameters. By decreasing cross-sectional area, the free-stream assumption at the boundaries becomes less valid. This introduces a blockage effect that promotes upstream pressure build-up and potentially increased velocity over the vehicle, both of which would increase downforce and drag, as was observed. Similarly, bringing the velocity-inlet closer the vehicle forces air that would previously be slowed down through conversion to increase static pressure to know retain its freestream velocity. This also boosts static pressure build-up and increases velocity over the vehicle. 

A surprising result was the insensitivity of the pressure-outlet distance to the rear of the vehicle. Only minor changes were observed up to as close as 4 m, which was about 30% of the required inlet distance to achieve similarly low sensitivity. Some additional tests were run with these short aft lengths while varying the other domain size parameters, mesh base size, and AMR settings, but the results were consistent and thus accepted (the low effect on mesh count led to 8 m being chosen as the final value to add a buffer to the unstable results observed with aft lengths lower than 4 m. Observing the total and static pressure distributions behind he vehicle, the best explanation for this insensitivity is that the wake is dominated by low energy recirculation close to the vehicle which slowly dissipates downstream. At this point, gentle rotation of the wake cloud is the most complex remaining behaviour; all coherent vortices are lost quickly to bluff wake. Combined with the lack of rear wing or diffuser to create more significant low static pressure and expand the flow aft of the vehicle, a pressure-outlet seems quite adept at managing the situation without affecting upstream flow, providing it remains aft of the recirculation that occupies the 3 m directly behind the vehicle. 

Calculating uncertainty

I've mentioned both scattered monotonic and oscillatory convergence so far. These are two potential types of convergence that can occur, and the uncertainty for each is calculated differently. Purely monotonic convergence is also an option, but this would likely only occur in a very simple problem and where mesh size is structured with perfect subdivisions to avoid the shifting of primary nodes as the mesh is refined. I expect this to be the cause of the oscillatory convergence seen for the mesh base size tests as seen in Figure 2, where the scaling of mesh size means nodes gradually shift closer to the origin as the mesh is refined. For small flow features like the diffuser vortices this can mean that unless the local mesh size is excessively small, details like the exact location of the vortex core can be dependant on the node location and thus the results will appear oscillatory as the nodes move in and out of the solver's "preferred" vortex core location (Figure 3). 

The method we used to obtain uncertainty estimates was generalised Richardson extrapolation, which at its core provides an estimate for the "true" solution, from which the difference to the accepted solution can be determined. Additional complexities exist to account for scatter and for adverse order of accuracy estimates. It was actually quite difficult to track down all the generally accepted methods and in the end I had to reference two documents (one of which referenced the other but had contradicting or at best confusing interpretations of the same method) as well as use a bit of logic.

For scattered monotonic convergence, a best-fitting curve of the form 

is generated based on a least-squares method. Here, δ is the error estimate for the ith value of the currently examined input parameter h, Φ₀ is the "true" solution estimate, Φ₁ is the solution for the ith value of h, a is a simple multiplier to be found, and p is the estimate for the order of accuracy, also to be found. In other words, it's a simple polynomial centred at h = 0. If p is between 0.95 and 2.05, this curve is kept and the uncertainty for the ith value of h is the sum of the error estimate and the standard deviation Uₛ of all points about the fitted curve, all multiplied by a recommended safety factor of 1.25. If p is below 0.95, a new curve is fitted of the form 

that forces an order of accuracy of 2 to avoid underfitting. The uncertainty estimate then depends on the best-fitting and re-fitted error estimates relative to the order of accuracy estimate for the best fit:

A similar process is followed (with different curve and uncertainty definitions) for values of p greater than 2.05, which can be found in the references above. These methods for monotonic convergence were clear enough, but the process for oscillatory convergence, while simpler, was not clear and recommendations differed. Oscillatory convergence is first identified by the number of times the solution values cross over the fitted curve, with the threshold being at least every third point on average. I already don't like this diagnosis method, but I'll come back to that later. As for the calculation, some suggest to take the difference between the largest and smallest solution values as a representative error estimate, while others specified that it should be the largest positive and negative errors between any two points and a best-fitting curve. Perhaps though these are practically similar, as the oscillation magnitude, if present, is typically more significant than any underlying monotonic convergence trend (hence there is no use of Richardson extrapolation in determine oscillatory uncertainty). A recommended safety factor of 3 is then applied, presumably to account for the lack of any consideration of such convergence. I'm not convinced by the method of taking a window between the largest positive and negative values though; Richardson extrapolation and standard deviations are both a magnitude of potential error, so should be applied as such with a +/- denotation. To me, better consistency between monotonic and oscillatory uncertainty interpretation would be achieved by taking the average of the largest positive and negative errors to provide a magnitude instead of the entire envelope. A safety factor of 3 also seems excessive, particularly in our case given the well-defined nature of the oscillatory convergence observed for mesh base size in Figure 2. I therefore chose to halve the recommended safety factor to 1.5. 

There is a lot of curve fitting and "if" statements going on in this method, and we have many uncertainties to calculate, so I wrote a Matlab app to take inputs of parameter values and results, and return the calculated uncertainty (Figure 4). To correctly assess convergence for domain sizes using this method, the lengths were raised to the power of -1 so greater distances to boundaries became smaller parameter values. Figure 5 shows an example output plot (mostly for debugging/sanity checks) for a convergence curve where the best-fitting order of accuracy was within the 0.95 to 2.05 window. Figure 5 shows an example output plot where the best-fitting order of accuracy exceeded 2.05, and an alternative curve is automatically generated to avoid overfitting. 

Figure 3: Diagram showing how incremental changes to mesh size can result in cells moving in and out of alignment with a specific location on a small flow feature, potentially creating oscillations in the results. 

Figure 4: Matlab app UI for automated simulation uncertainty calculations

Figure 5: Output plot from the Matlab script for Cd*A across various values for domain height/width, as an example of a case where the order of accuracy estimate provided by the best fit is deemed appropriate. This is also an example of a case where, strictly, convergence is oscillatory, but I chose to intervene and force the use of scattered monotonic calculations. 

Figure 6: Output plot from the Matlab script for Cd*A across various values for domain aft length, as an example of a case where the order of accuracy estimate provided by the best fit is deemed too high and at risk of overfitting

The consequence of the rather crude oscillatory convergence assessment is that more often than not, the behaviour is diagnosed as oscillatory when logically and visually it is better explained by scatter/noise. For example, in a case where the "oscillations" are inconsistent in magnitude and period, and where there is no reason to believe that the input parameter being adjusted would create cause oscillatory behaviour (such as the domain size parameters), it makes more sense to attribute this behaviour to noise (Figure 7). True oscillatory convergence should be sinusoidal in nature; magnitude and period may still change, but are not random. There were a number of cases where I intervened with the automated oscillation detection and forced the uncertainty to be calculated based on monotonic convergence. In the end I decided to just remove this functionality from the app and have the user always judge for themselves the convergence type.

Adjustments for the skidpan case

The skidpan tests were conducted the same way for the most part. To save time, the uncertainty and best value of domain height, fore length, and aft length were carried over from the straight airflow ride-height case, which was assumed to be valid and conservative given the lower freestream velocity in the skidpan case. The fore and aft lengths were specified as the arc length at the corner radius (9 m for skidpan). The inner and outer domain radii were rigorously tested (based on radial distance to the corner radius arc), as was mesh base size. Logically, the inner radius would have negligible effect, however there are detrimental effects associated with having the velocity inlet and pressure outlet too close to each other, resulting in some simulations with small inner radii diverging. 

In general, the convergence trends were similar, but standard deviations of scatter were increased across the board whenever monotonic convergence was present. This was not unexpected given the more complex flow interactions created by the cross flow and freestream velocity gradient. 

Figure 7: Comparison of (top) scattered monotonic vs (bottom) oscillatory convergence. Both are assessed as oscillatory by the simple crossing-the-curve method discussed above, but the inconsistent magnitude and period of "oscillations" in the top plot make it clear that it is just scatter/noise.

Sensitivity tests

A few separate simulations were run to assess sensitivity of both initial conditions and AMR cut-off values. 

Three initial condtions were tested using a setup with the final selected domain and mesh size parameters:

All three produced identical results for all forces and monitors within 0.1%

Different AMR functions have been tested throughout last year, but the original four as mentioned in CFD from the ground up remain the most efficient and effective. Some minor changes to cut-off values were made though:

Final values

The uncertainties for each parameter were combined through root mean squares to give the total uncertainty for each of lift, drag, and balance, as listed below:

Straight airflow ride-height

ClA: +/- 0.030 | 2.8%
 CdA: +/- 0.007 | 1.1%
 Bal: +/- 0.013

Skidpan

ClA: +/- 0.031 | 2.7%
 CdA: +/- 0.010 | 1.60%
 Bal: +/- 0.016

It makes sense that drag would have a relatively lower uncertainty compared to lift as its nearly always only acting in one direction (-Z), so its total magnitude is proportional to the individual static surface pressure magnitudes acting in the Z axis. With lift, some parts of the car will be making lift while others make downforce. The total static pressure magnitudes acting in the Y axis therefore add up to much greater than the net lift force, giving more room for sensitivities to creep in behind the façade of just one overall lift value. I'm reluctant to use the percentage values for this reason, which also applies to aero balance (particularly since a value of 0 is not a neutral balance; far from it in fact). We'll instead take a conservative approach of using the maximum of either magnitude or percentage when applying these uncertainties to future results. Significant changes to surface pressures or flow field behaviour will warrant a new uncertainty analysis, which at this point will likely be the planned addition of new aerodynamic components for 2026. 

Also interesting is the similar lift coefficient uncertainties for both cases; it was expected that the skidpan case would be more sensitive. Generally though, while noise/scatter was higher in the skidpan case, oscillation amplitudes (when present) were lower than those in the straight airflow tests. 

The final selected parameter values produced the following mesh counts (initial count, and final count after 4 AMR updates):

Straight airflow ride-height

13.5M | 19.0M

Skidpan

26.5M | 36.0M

Given the similar uncertainties between both cases, it seems appropriate to stick with the original plan of applying the straigh airflow ride-height results to the braking case, and the skidpan results to the median corner case, without conducting a rigorous analysis on either. Updated AMR values for the median corner case will be tested separately though, given it has a unique freestream velocity. 

Implications

There's no point going to so much effort if we can't justify it. Uncertainties are pointless if they are not given context or application. We consider the effort justified through the application of these rigorous uncertainty estimates to three aspects that we deem important: