---
jupyter:
  description: 'The Kutta-Joukowski theorem with planex2d and NGSolve: a Joukowski airfoil drawn exactly from its
    conformal map, potential flow around it with the Kutta condition at the trailing edge, and the circulation and the
    lift coefficient against the exact values.'
---

# The Joukowski airfoil

How much lift does a wing section give? In potential flow, the Kutta-Joukowski theorem gives the lift per unit span as
$L' = \rho U \Gamma$, the density times the speed times the circulation $\Gamma$ around the section. The Kutta
condition fixes the circulation: the flow leaves the sharp trailing edge smoothly instead of turning around it. For a
Joukowski airfoil, the image of a circle under the conformal map $z = \zeta + b^2/\zeta$, the circulation is exact:

$$\Gamma = 4\pi U a \sin(\alpha + \beta),$$

with $a$ the circle's radius, $\alpha$ the angle of attack and $\beta$ the angle at which the circle passes through
$\zeta = b$, the point that becomes the trailing edge. The lift coefficient is $C_L = 2\Gamma/(U c)$ for the chord $c$.

Here, `px.Curve` draws the airfoil exactly from its map, NGSolve solves the potential flow around it, and the
circulation is checked against the exact one.

![Streamlines around a Joukowski airfoil at 8° angle of attack](../_static/gallery/joukowski-airfoil.png){.result}

```python
import math

import ngsolve as ngs
import numpy as np
import planex as px
```

```python nbsphinx="hidden"
from docs_figures import Row
```

## The airfoil

The circle has its center at $(-0.1, 0.1)$ and passes through $\zeta = b = 1$, so its radius is $a = \sqrt{1.1^2 +
0.1^2}$ and $\beta = \arctan(0.1/1.1)$. The real part of its center makes the airfoil thick and the imaginary part
cambers it. With $\zeta = (-0.1 + a\cos t) + i\,(0.1 + a\sin t)$, the map gives the airfoil as two expression strings in
$t$, starting at the trailing edge, $t = -\beta$ (Figure 1):

```python
b = 1.0
mx = 0.1
my = 0.1
a = math.hypot(b + mx, my)
beta = math.atan2(my, b + mx)
u = f"(-{mx} + {a}*cos(t))"
v = f"({my} + {a}*sin(t))"
x = f"{u} + {b * b}*{u}/({u}^2 + {v}^2)"
y = f"{v} - {b * b}*{v}/({u}^2 + {v}^2)"
start = -math.degrees(beta)
foil = px.Curve((x, y), (start, start + 360))
(x0, _), (x1, _) = foil.bounding_box
chord = x1 - x0
print("chord", chord, "| trailing edge at", x1)
```

```python nbsphinx="hide-input"
Row((px.Face([foil]), "The Joukowski airfoil, its cusp at the trailing edge"), caption="A cambered Joukowski airfoil")
```

## The flow

The stream function $\psi$ of a potential flow satisfies $\nabla^2 \psi = 0$, and the airfoil, which the flow does not
cross, has a constant $\psi$ on it. Far away, on a circle of radius 50, the flow is uniform at the angle of attack.
The solution is the uniform flow around the airfoil, $\psi_s$ with $\psi = 0$ on the airfoil, plus a multiple of the
pure circulation, $\psi_c$ with $\psi = 1$ on the airfoil and 0 far away:

```python
body = px.Face([foil])
body.edges.name = "foil"
body.edges.maxh = 0.01
air = px.Circle((0, 0), 50).Face()
air.edges.name = "far"
mesh = ngs.Mesh(px.OCCGeometry((air - body).Face()).GenerateMesh(maxh=5, grading=0.2))
mesh.Curve(3)
V = ngs.H1(mesh, order=3, dirichlet="far|foil")
p, q = V.TnT()
K = ngs.BilinearForm(ngs.grad(p) * ngs.grad(q) * ngs.dx).Assemble()
inverse = K.mat.Inverse(V.FreeDofs(), inverse="umfpack")


def potential(far, on_foil):
    psi = ngs.GridFunction(V)
    psi.Set(mesh.BoundaryCF({"far": far, "foil": on_foil}), ngs.BND)
    psi.vec.data += inverse * (-K.mat * psi.vec)
    return psi
```

## The Kutta condition

A flow that turned around the sharp trailing edge would do so at an infinite speed. The Kutta condition picks the
multiple $k$ of $\psi_c$ for which the flow just behind the cusp runs along the trailing edge's direction, the angle
$-2\beta$, instead of across it. The circulation is then the flux of $\nabla\psi$ through the airfoil's outline:

```python
edge = np.array([x1, 0.0])
along = np.array([math.cos(-2 * beta), math.sin(-2 * beta)])
across = np.array([-along[1], along[0]])
behind = edge + 3e-4 * along
normal = ngs.specialcf.normal(2)


def velocity(psi):
    gx, gy = ngs.grad(psi)(mesh(*behind))
    return np.array([gy, -gx])


for alpha in (0, 4, 8):
    rad = math.radians(alpha)
    stream, circulation = potential(ngs.cos(rad) * ngs.y - ngs.sin(rad) * ngs.x, 0), potential(0, 1)
    k = -(velocity(stream) @ across) / (velocity(circulation) @ across)
    gamma = abs(ngs.Integrate(ngs.BoundaryFromVolumeCF(ngs.grad(stream) + k * ngs.grad(circulation)) * normal, mesh,
                              ngs.BND, definedon=mesh.Boundaries("foil")))
    exact = 4 * math.pi * a * math.sin(rad + beta)
    print("α", alpha, "° | Γ", round(gamma, 5), "exact", round(exact, 5), "| C_L", round(2 * gamma / chord, 4),
          "exact", round(2 * exact / chord, 4))
```

The circulation and the lift come within 0.2 % of the exact values, and the lift grows with the angle of attack as
$\sin(\alpha + \beta)$. At zero angle of attack, the camber alone gives lift. The streamlines at 8° leave the trailing
edge smoothly:

```python nbsphinx="hide-input"
import matplotlib.pyplot as plt
from docs_figures import FACE, INK

psi = stream + k * circulation
gh, gv = np.meshgrid(np.linspace(-4, 4, 401), np.linspace(-2.2, 2.2, 221))
outside = (air - body).Face().Distance(np.c_[gh.ravel(), gv.ravel()], signed=True) < -1e-4
values = np.full(gh.size, np.nan)
values[outside] = np.asarray(psi(mesh(gh.ravel()[outside], gv.ravel()[outside]))).ravel()
fig, ax = plt.subplots(figsize=(7, 3.6))
ax.contour(gh, gv, values.reshape(gh.shape), levels=np.linspace(-2.4, 2.4, 33), colors=INK, linewidths=0.7,
           linestyles="solid")
body.Plot(ax=ax, face_col=FACE, line_col=INK, spline_col=INK, legend=False)
ax.set_aspect("equal")
ax.axis("off")
plt.show()
```

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
from gallery import save_figure
save_figure(fig, "joukowski-airfoil", bare=True)
```