OpenFOAM airfoil optimization

Overview

This is a simple set of code that automatically finds airfoils, optimizing for the lift-to-drag ratio $C_l / C_d$. It does this by generating an airfoil shape based on six CST-parameters (courtesy of this repo), using Kulfan's 2007 paper; Universal Parametric Geometry Representation Method. An initial attempt to create a variable meshing code myself with blockMesh turned out to be very painful, so the meshing is handled by curiosityFluids' excellent mesher (blog post)

SciPy's differential evolution is taken as an optimization algorithm. It's far slower than other methods, but using a global optimizer here seems like the better choice, since I want to explore the full space, to see if there are multiple viable solutions. Other derivative-free algorithms like Nelder-Mead also found reasonable airfoils and were much faster, however.¹.

The simulation is then ran. If any issues are encountered with the meshing, blockMesh, or simpleFoam, or convergence, the code returns $+\infty$. For $C_l$, any value is considered feasible, even negative ones - that turns out to generally be the code inventing upside-down airfoils. Negative $C_d$ are penalized, however, since this generally indicates the optimizer cheating, or the case not having converged. The value to minimize is then taken as $-|C_l / C_d|$.

The result is a CSV containing airfoil parameters and their performance. These can be further post-processed with ParaView.

This was a fun Christmas holiday project, and a nice foray into coupling non-trivial simulation problems with (surrogate-based) optimization. If anyone has suggestions on how to improve or adjust things further - I am very much open to them! Overall, I'm surprised at how smoothly this project went. The existing repos helped a lot, especially with meshing. I'm still impressed at how effective differential evolution was - with a previous meshing-template, it was able to find and exploit flaws with ease. I had to adjust the goal function so many times there. I am also very much impressed with how effective random forests were at representing these complex simulations in the surrogate model aspect!

Installation

This assumes you already have OpenFOAM installed - I am using the OpenFOAM.com/ESI version. You can clone the repo. Source the OpenFOAM functions first. For the ESI-version, this can be done with

source /usr/lib/openfoam/openfoam2406/etc/bashrc

This exposes functions like blockMesh and simpleFoam. Dependencies were kept to a minimum: pandas, numpy, matplotlib, scipy, and optionally sklearn. Once these are installed, run python main.py.

Getting started

The code is quite minimal. This is partially on purpose. If you are using OpenFOAM, you are used to editing code files. It's also quite difficult to create a sufficiently flexible way of working with these things without oversimplifying. The intended way of using the code is to look through everything and try to understand it. Everything should be quite straight-forward and readable.

Features

Parallel-processing

The code uses an OpenFOAM template folder, which is repeatedly copied into several directories, depending on the number of workers specified, using UUIDs as names. Relevant parameters are then entered into these templates (such as the adjusted blockMesh, and U) As soon as a case completes (either because of errors, or because it correctly finished), the folder is deleted. Since OpenFOAM requires its folders in a particular structure, differential evolution can easily be parallelized, and because each individual case runs fast, this seemed like a better choice than decomposing the domain across multiple cores. With the current method, there should be almost no overhead from parallelization.

Automatic top-n selection and rendering

The code can post-process the top-n highest scoring airfoils so far, simulating each of them. The results can subsequently be rendered with a ParaView Python macro. First, run

python main.py --custom

This places the top-n runs under custom_runs. Follow this by

python src/post_processing/post_process.py

If all goes well, this places all the results under results/renders.

The latter isn't entirely reliable: ParaView includes its own Python-distribution based on Python 3.10. If you encounter errors here, it's best to either look at the code, or open each .foam under custom_runs individually.

Analysis

The repo includes an analysis-notebook. This tracks performance over time, and allows for selection of the best-performing airfoils.

Updates and adjustments

Turbulence modelling

I started off with the default parameters from a case I found, which used $\tilde \nu = 0.14$. That seemed high, also compared to other OpenFOAM examples I found, but I assumed it was acceptable - up until I started getting $C_l/C_d$ values of over 900. Reducing the layer thickness close to the mesh to fix $y+$ did not improve the situation.

It turns out, $\tilde \nu$ was probably far too high. There are different sources here. CFD-Online argues $\tilde \nu$ should be 0 for the Spalart-Allmaras model (which I'm using), but can also be set to $\tilde \nu \leq \nu / 2$ if that causes problems with solvers. Since $\nu \approx 14.88 \cdot 10^{-6}$ at ambient temperatures and pressures, this is certainly a far lower value than I was using before. A NASA source gives different values still. (Note that this page uses $\hat \nu$ rather than $\tilde \nu$; the latter apparently didn't show up properly on screens when this page was created.) It says to use $\tilde \nu_\text{far-field}$ between $3 \nu_\infty$ and $5 \nu_\infty$, so between $4.5$ and $7.5 \cdot 10^{-5}$.

I tested some of these for one of the airfoils that previously messed up, and it appears there is almost no difference between the smaller values.

Case	$\tilde \nu$	$C_l$	$C_d$	$C_l/C_d$
Original	$0.14$	0.516667	-0.000555	-931.013514
OpenFOAM default	$4.0 \cdot 10^{-5}$	1.13645	-0.0742927	-15.295453
Lower NASA-bound	$4.5 \cdot 10^{-5}$	1.1369	-0.0743574	-15.287780
Far lower value	$7.0 \cdot 10^{-6}$	1.13665	-0.0743269	-15.290234

I will continue to use $4.5 \cdot 10^{-5}$ for now. If anyone has better suggestions, I would love to hear them!

Convergence and populations

Another interesting issue; some airfoils never converged with SIMPLE, instead oscillating at different $C_l/C_d$ values. The intermediate SIMPLE steps looked a bit similar to a Kármán vortex street, and more iterations did not appear to help. I resolved this simply by specifying a bound on $\sigma_{C_l/C_d}$ and let differential evolution handle it.

With that, we get an interesting population. Three out of the top-four are all very different; the best performer is a fairly standard airfoil, albeit a bit thick. The next-best is almost bird-like, and the third is high-camber instead. It's surprising to see such variation even after a fairly long run.

53.156

52.751

51.072

50.497

$C_l/C_d$ for different archetypes in the top-4 airfoils at 5° AoA.

$C_l/C_d$ curves

I added code to evaluate the lift-drag ratio as a function of the angle-of-attack (AoA). The fact that the airfoil performs best at 5° isn't very surprising; it's optimized for that point. For higher angles, performance rapidly decreases. I considered multi-objective optimization to create an airfoil that performs well over a wider range, but that would be very slow to run.

One interesting aspect here: I ran this up to 45°, but did not obtain sufficiently converged solutions there. Examining the forces for the highest AoA where simpleFOAM did not simply crash, we observe an oscillatory solution - I think this is basically the air detaching from the airfoil, resulting in a kind of Kármán vortex street, with SIMPLE unable to converge to a steady-state solution. Examining ParaView here, we indeed see oscillations in the flow field.

This behavior previously caused a lot of issues with very small or negative $C_d$ values, until I simply specified that any solution where $\sigma_{C_l/C_d} \geq 1$ over the last 500 steps is disqualified, to let differential evolution take care of it.

Model reduction

I am curious about potential model reduction: by predicting performance based on the six inputs, a lot of time could be saved. If a rough prediction on which airfoils perform best is accurate, a simple machine learning model like random forests could be used for an initial optimization stage. I doubt a simple model like this would be sufficient, but it's an interesting avenue to explore.

After some attempts, it seems surprisingly good. I get MAEs of 1.5 - 5, for a small training set, even for the best-performing airfoils, where there is correspondingly less data available. This is technically a surrogate model-approach.

After running it a bit longer, it gets better and better; I'm very surprised. We do have a fair amount of data, but this is spread out in six dimensions; the curse of dimensionality should be kicking in here, yet somehow, even with quite sparse data, it's doing well. However, I am not using a randomly sampled set; the data is all from an optimizer, so it's likely to be clustered around certain regions, effectively reducing dimensionality.

Optimizing over the surrogate models

I added a grid-search with 5-fold cross-validation to optimize a classification and regression model, then optimized those with the same differential evolution code. This is a two-step process; we first predict whether we will have any result at all (i.e. no failures in overlapping airfoils, blockMesh, simpleFoam, or convergence issues), and if the random forest predicts there aren't, we regress our vector to obtain $C_l/C_d$.

Oddly, this gets stuck below the best-performers that we previously found using the regular optimization method. Overall, though, it's very similar in shape and design to the optimal version, and it's close; the best airfoil I have found thus far reached 59.68698, and this one is at 59.31823; it's very close, and it only took a few minutes to run, compared with 72 hours for the full model.² After some thought, the reason for this is obvious; a random forest partitions space with constant values, but cannot extrapolate. A switch to an SVM for regression gives better results; it now hits $C_l/C_d = 60.84197$ as predicted by the SVM, with simulations estimating it at $C_l/C_d = 60.38784$ - noticeably better than the regular optimization was able to find in three days, after running for only fifteen minutes.

Footnotes

1. We can't use derivative-based optimization methods very easily, because if we have a crash in blockMesh, we have no measure of how 'badly' things messed up, so steering gradients away from there is difficult. ↩

2. The figure name has 58.57 in there, but I messed up the naming; that was for a previous, worse airfoil. ↩