Conducting a convergence and uncertainty analysis
What I'm covering in this post is verification of the EV24 CFD setup featuring the side-floors. In particular, the process of conducting an uncertainty analysis to determine the effective numerical precision of our simulation results based on sensitivity and convergence. 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, and the car having a very simple and thus naturally less sensitive 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 waited 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 whenever there is a major change to the flowfield around the car. Doing it properly meant setting up dedicated field probes, reading up on uncertainty procedures and formulas, and creating a program to automate calculations.
For convergence tests, there were two key parameters we were intested in: mesh and domain sizes. Each would have an ideal value and an uncertainty estimate associated with it. 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 used to assess convergence. Additional to these were lift-drag ratio, six total pressure monitors at critical and/or sensitive locations, and two static pressure reports to sample the maximum positive and negative pressure magnitudes anywhere on the vehicle.
Total pressure monitor locations:
Brow wing inboard vortex core, 300mm aft of trailing edge
Brow wing outboard vortex core, 300mm aft of trailing edge
Diffuser inboard vortex core, 300mm aft of shedding edge
Diffuser outboard vortex core, 300mm aft of shedding edge
Inboard "edge" of front tyre wake, 300mm aft of brow wing trailing edge and 300mm above the ground
In the wake of the vehicle, 2m behind the rear axle and 0.5m above the ground at the centreline
As well as domain and mesh sizing, AMR cut-off values were assessed for sensitivity. A full AMR convergence test was deemed unecessary since the sensitivities could be reduced far below those of the global mesh size with only a tiny increase in total mesh count, such is the efficiency of AMR functions compared to basic geometric refinements.
A "baseline" setup was decided on for both the straight airflow and skidpan cases. These baselines used mesh and domain sizes that were good enough to minimise cross-influence between their effects, but while still keeping an acceptable solve time. From here, each parameter was iteratively coarsened until divergence dominated over 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. Finer-than-baseline values were used if the baseline results did not appear converged.
Other notes on the method:
Domain width and height were tested simultaneously to reduce the total number of verification simulations required.
All domain sizes for the straight airflow case are measured from the origin to the respective boundary.
Inner and outer wall distances for the skidpan case were measured from the vehicle centreline, extending radially from/toward the corner centre.
Pressure monitor convergence was assessed visually only.
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). The green cells represent the final selected value for this parameter.
All main convergence plots for the straight airflow case
Most results for the straight airflow case exhibited good convergence curves with varying degrees of scatter/noise. There are a few interesting point to mention:
The sudden divergence of lift and drag for variations in domain width/height suggest a high order of accuracy. In isolation, this is not unexpected given the simultaneous adjustment of width and height changing area rather than just length. However, the behaviour of the aero balance and lift-drag ratio plots at this right-most point clearly depart from the trend up to that point, suggesting that the results have instead diverged beyond the asymptotic range, where it is no longer valid to compare the results. 12 m was selected as the final value.
The sudden divergence of lift and drag for variations in domain aft length suggest divergence beyond the asymptotic range. 8 m was selected as the final value.
There is oscillatory convergence across all results for the mesh base size parameter. The shape or even presence of any underlying convergence curves is unclear. 8.3 mm was selected as the final value, which in this case was based more on the sensitivity and convergence of the pressure monitors.
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, a blockage effect is introduced 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 to the vehicle forces free-stream velocity where it would previously have been reduced by the pressure buildup infront of the vehicle, thus having much the same effect as the width/height parameter.
A surprising result was the insensitivity of the pressure-outlet ("Domain Aft") distance behind the vehicle. Only minor effects were observed up to as close as 4 m, which was about 30% of the distance found ideal for the inlet. A value of 8 m was selected to give a bit of a buffer, since domain size has a minimal effect on mesh count and thus solve time. A potential explanation for the insensitivity is the generally streamlined shape of the car, resulting in quite gentle pressure gradients in the wake region that are insensitive to a nearby imposed zero-pressure boundary. I expect the inclusion of a rear wing would require a longer aft distance.
The method we used to obtain uncertainty estimates was generalised Richardson extrapolation, which at its core provides an estimate for a fully-converged solution. Additional functions exist to account for scatter and correct 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).
ITTC, 2017. ITTC Quality System Manual, Recommended procedures and guidelines 7.5-03-01-01: Uncertainty analysis in CFD verification and validation, methodology and procedures
Eca, L., and Hoekstra, M., 2008. Testing Uncertainty Estimation and Validation Procedures in the Flow Around a Backward Facing Step. 3rd Workshop on CFD Uncertainty Analysis, Lisbon.
For scattered monotonic convergence (as observed/assumed for all domain size tests), 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, Φi is the solution for the ith value of h, a is a simple multiplier to be solved for, and p is the estimate for the order of accuracy, also to be solved for. 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 R.E. error estimate and the standard deviation 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 both the best-fitting and re-fitted error estimates:
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.
For oscillatory convergence, some suggest to take the difference between the largest and smallest solution values as a representative error estimate, while others specify that it should be the largest positive and negative errors between any two points and a best-fitting curve. The underlying convergence curve is assumed negligible in both cases and thus there is no use of Richardson extrapolation. Since the convergence curve is assumed negligible, the two different methods of acquiring an uncertainty estimate are unlikely to yield significantly different values. A recommended safety factor of 3 is then applied.
I decided to only count the single largest error magnitude, as the uncertainty was to be expressed as a deviation (i.e., "+/-") rather than a range, so this method seemed more appropriate.
To save time calculating all the uncertainties, I wrote a Matlab app to take inputs of parameter values and results, and return the calculated uncertainty.
Matlab app UI for automated simulation uncertainty calculations
The skidpan tests were conducted the same way with some different definitions for domain sizes.
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 specified as radial distance from the corner radius arc.
In general, the convergence trends were similar, but standard deviations of scatter were typically higher than the straight airflow case. This was not unexpected given the more complex flow interactions created by the cross flow and freestream velocity gradient. Another interesting note was that while setting inner domain radius to 0 would be acceptable on a cell-count basis, having the radius too low put parts of the inlet and outlet boundaries in close proximity which led to solver divergence.
A few separate simulations were run to assess the sensitivity of results to 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:
uniform velocity of 0
uniform freestream/target velocity
solution field from a prior simulation
Results from all sensitivity tests returned differences well within the uncertainties already calculated.
Different AMR functions were 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 threshold values were made though:
The turbulent kinetic energy refinement was relaxed behind the front wheels, and disabled entirely behind the vehicle.
Results showed a low sensitivity to the velocity gradient refinement threshold but also the least mesh-heavy, so was left as-is.
Results were insensitive to the total pressure refinement threshold until it increased by over 20%, at which point there was excessive loss of downstream front tyre, rocker, and brow wing wake definition, extending to lift and drag estimates. It was therefore kept at its previous value.
The vorticity refinement was maintained. An additional finer refinement level was applied, limited to within 120mm of the ground to target the vortex cores under the floor. This gave more consistent results and a reduction in the oscillation amplitudes of the mesh base size convergence tests when reconducted.
The uncertainties for each parameter were combined through root mean squares to give the total uncertainty for 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 lower uncertainty than lift as the associated forces are almost exclsively acting in one direction. 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, allowing more sensitivities to hide behind the lift estimate.
The uncertainties for the straight and skidpan cases will also be applied to the braking and median corner cases respectively, using the percentage values.
The final selected parameter values produced the following mesh counts (initial count, and final count after the usual 3 AMR updates):
Straight airflow ride-height
13.5M | 19.0M
Skidpan
26.5M | 36.0M
There's no point going to so much effort to calculate all this if we don't actually use it. Uncertainties are pointless if they are not given context or application. We consider the effort justified through the application of the uncertainty estimates to three aspects that we deem important:
The process allows numerical comparisons of the uncertainty introduced by adjustments to mesh and domain sizes, allowing us to informedly select the most time-efficient setup with respect to results repeatability.
It provides a "required improvement" threshold to confidently say if a CFD design iteration is better than a previous result. Our lack of wind tunnel access for finer optimisation means having a low uncertainty is particularly valuable to us, as we want to be able to compare similar setups in CFD and fine tune things normally done in a wind tunnel.
When combined with experimental uncertainty, the uncertainty estimates act as a window to qualify experimental results as being accurately or inaccurately predicted by CFD