Content provided under a Creative Commons Attribution license CC-BY 4.0; code under BSD 3-Clause license. (c)2015 L.A. Barba, Pi-Yueh Chuang.

2D Multi-element airfoil

This is your AeroPython assignment for the fourth course module, titled "Vortex-panel method for lifting bodies." You will investigate a 2D multi-element airfoil, or wing section with extended flap, adapting the source-vortex panel method of Lesson 11 to compute the flow around two airfoils.

The key concepts for using the vortex-source panel method with a multi-element airfoil are:

  1. each airfoil has its own vortex strength, and
  2. the trailing edge at each airfoil must follow the Kutta condition.

In the first part of the assignment, you will use a test for which we have a theoretical solution, due to Williams (1973). Have a look at the reference, and get a feel for its mathematical genius! In the second part, you will use a more popular airfoil (the NACA 23012) for both the main wing section and the flap. For this case, there is an old classic NACA report by Wenzinger (1938) giving data from experiments in a wind tunnel. This exercise will get you wondering about some very interesting aerodynamical questions!

Part 1: Williams' test case

Your task is to calculate the pressure coefficient on a multi-element wing section, consisting of a main airfoil and an external-airfoil flap. This test is from Williams (1973), who obtained a beautiful theoretical solution using the techniques of conformal mapping.

The profile of the wing section looks like this, with flap extended:

We provide CSV files with the $x$ and $y$ coordinates for the end nodes of each panel; these files are in the resources directory of our GitHub repository. The files have the following naming convention, where *** represents the value of N, the number of panels: MainFoil_N=***.csv, FlapFoil_N=***.csv.

You can use fewer panels for debugging while developing your codes and then use more panels for more serious calculations when answering the questions in this assignment.

There are also two files named Cp_Main_theoretical.csv and Cp_Flap_theoretical.csv containing the theoretical values of pressure coefficient versus $x$-coordinates, for the flap configuration we give you in the data files (a $30$º flap deflection, with the main airfoil at $0$º angle of attack).

The theoretical lift and drag force for potential flow over this 2D wing section with zero angle of attack are $3.7386$ and $0$, respectively (non-dimensionalized by dynamic pressure).


Recall the exercise in Lesson 11, which asks you to derive the mathematical formulation of the vortex-source panel method.

Now, for the 2-element airfoil, you need to derive the required mathematical expressions, following the same process as in the exercise, starting from the following expression for the potential:

$$ \begin{split} \phi(x, y) &= U_{\infty}x\cos\alpha + U_{\infty}y\sin\alpha \\ & + \int_{main} \frac{1}{2\pi} \sigma(s) \ln \sqrt{(x-\xi(s))^2+(y-\eta(s))^2} ds \\ & + \int_{flap} \frac{1}{2\pi} \sigma(s) \ln \sqrt{(x-\xi(s))^2+(y-\eta(s))^2} ds \\ & - \int_{main} \frac{1}{2\pi} \gamma(s) \tan^{-1} \frac{y-\eta(s)}{x-\xi(s)} ds \\ & - \int_{flap} \frac{1}{2\pi} \gamma(s) \tan^{-1} \frac{y-\eta(s)}{x-\xi(s)} ds \end{split} $$

Assume the following:

  1. $\sigma(s)$ is constant on each panel,
  2. $\gamma(s)$ is constant on each airfoil, and
  3. the 1st to $N_a$th panels are on the main airfoil and the $(N_a+1)$th to $N$th panels are on the flap (where $N=N_a+N_b$, and $N_b$ is the number of panels on the flap).

There will be $N+2$ unknowns, that is, $\sigma_1\cdots\sigma_N$, $\gamma_a$, and $\gamma_b$. The values of $\gamma_a$, and $\gamma_b$ represent the vortex strengths on the main and flap airfoils, respectively.

You should be able to obtain the following matrix forms of normal and tangential velocity on the $i$th panel:

$$ U^n_i = b^n_i + \left[\begin{smallmatrix}A^n_{i1}\cdots A^n_{iN}\end{smallmatrix}, \sum_{j=1}^{N_a}B^n_{ij}, \sum_{j=N_a+1}^{N}B^n_{ij}\right] \left[\begin{smallmatrix}\sigma_1 \\ \vdots \\ \sigma_N \\ \gamma_a \\ \gamma_b\end{smallmatrix}\right] $$$$ U^t_i = b^t_i + \left[\begin{smallmatrix}A^t_{i1}\cdots A^t_{iN}\end{smallmatrix}, \sum_{j=1}^{N_a}B^t_{ij}, \sum_{j=N_a+1}^{N}B^t_{ij}\right] \left[\begin{smallmatrix}\sigma_1 \\ \vdots \\ \sigma_N \\ \gamma_a \\ \gamma_b\end{smallmatrix}\right] $$

Using the non-penetration conditions on the $1$st to $N$th panels, you will have $N$ linear equations. To solve for the $N+2$ unknowns, the extra two linear equations you need are the Kutta conditions on the two airfoils:

$$ U^t_{1} = U^t_{N_a} \\ U^t_{N_a+1} = U^t_{N} $$

And you are ready to solve for potential flow around a 2D multi-component wing!


Run your codes using 100 panels on both main and flap airfoils and answer the following numeric-value questions.


  • What is the pressure coefficient on the $3$rd panel on the main airfoil?


  • What is the pressure coefficient on the $10$th panel on the flap airfoil?


  • What is the value of lift?

Hint: $~L = - \oint_{main}p\vec{n}\cdot\vec{j}dl- \oint_{flap}p\vec{n}\cdot\vec{j}dl$


  • What is the value of drag?


  • Run the program again with 200 panels on both main and flap airfoils. What is the percentage of improvement of the absolute error in the value of lift?

Part 2: Wenzinger's test case

You should now have a usable panel solver for 2-component airfoils. Next, you'll change the target airfoil to a more realistic profile. There are experimental results available (Wenzinger, 1938) for a NACA 23012 airfoil with a NACA23012 external-airfoil flap, which we'll compare with.

The following figure shows the profile of this wing section with 0º flap angle:

We provide CSV files for the endpoints of a panel discretization of this wing section, using 150 panels on each airfoil: NACA23012_MainFoil.csv and NACA23012_FlapFoil.csv in the resources folder.

The flap has a rotating center (hinge) at the location $(1.03, -0.054)$ (see the configuration on Figure 1 of the reference). Using coordinate rotation, you can obtain different configurations with different flap deflection angles.

Try to use different flap deflection angles and different angles of attack (of the main airfoil) and compare to Wenzinger's experimental results. Note that the definition of the total lift coefficient is $L/(l_{main}+l_{flap})$, where $L$ is the lift force (per unit span) like in the previous problem, and $l_{main}$ and $l_{flap}$ are the chord lengths of the two airfoils. Think about what could be the sources of difference between your results and the experimental data.



  • Use the bisection method to determine the angle of attack, $\alpha$, that will cause zero lift when the flap angle is $20$º. Round your answer to the second digit after the floating point.

HINT: you can start the bisection method within the range $-14\le\alpha\le14$.


  • Use the bisection method to determine the flap angle that will cause zero lift when the angle of attack, $\alpha$, is $0$º. Round your answer to the second digit after the floating point.


  • Suppose that the main airfoil is attached to the fuselage of an airplane at an incidence angle of $4$º (that is, the built-in angle of attack of the wing when the airplane is in level flight is $4$º). Compute the lift coefficient with the flap at $0$º deflection, and at $5$, $10$ and $15$ degrees. How much slower (in percentage) can the airplane fly (in level flight) with the flaps extended at these settings, compared with $0$º?

Think about the meaning of what you observed in Q3, above. What is the effect of using the flap?


  • B. R. Williams (1973), An Exact Test Case for the Plane Potential Flow About Two Adjacent Lifting Aerofoils, Reports & Memoranda No. 3717, Aeronautical Research Council of the United Kingdom // PDF at Cranfield University

  • C. J. Wenzinger (1938), Pressure distribution over an NACA 23012 airfoil with an N.A.C.A. 23012 external-airfoil flap, NACA Technical Report No.614 // PDF at NASA

In [1]:
from IPython.core.display import HTML
def css_styling(filepath):
    styles = open(filepath, 'r').read()
    return HTML(styles)