Potential-flow theory has been used to study the flow over 2D airfoils for a long time. Using complex analysis, we can obtain analytical solutions for many potential flows, especially the flow over various 2D airfoils. When computers became available, people could also solve potential flows numerically using panel methods, finite-difference methods, boundary-element methods, etc.

Though there are many unrealistic assumptions in potential flow, the theory provides a good beginning for further analyses. For example, some CFD codes can solve potential flow first, in order to provide a good initial guess of viscous flow solver.

In 2D potential flow, we can solve for either the potential function or stream function through their governing equtaions: the Laplace equation. The governing equation for stream function is $$ \nabla^2\Psi=\frac{\partial^2 \Psi}{\partial x^2}+\frac{\partial^2 \Psi}{\partial y^2}=0 $$ where $\Psi$ is the stream function. If the velocity vector $\vec{V}=u\vec{i}+v\vec{j}$ represents the flow velocity, the relationship between flow velocity and the stream function is: $$ u=\frac{\partial \Psi}{\partial y}\\ v=-\frac{\partial \Psi}{\partial x} $$

In this notebook, we show how to solve potential flow over 2D airfoils using the finite difference method.

We use NACA 4-digit airfoils for our example in this notebook. For an introduction to NACA 4-digit airfoils, refer to reference [1].

Let's first define a Python function that can be used to generate the thickness profile of a NACA airfoil, given by：

$$ y_t = 5t\left[ 0.2969\sqrt{\frac{x}{c}} -0.126\left(\frac{x}{c}\right) -0.3516\left(\frac{x}{c}\right)^2 +0.2843\left(\frac{x}{c}\right)^3 -0.1036\left(\frac{x}{c}\right)^4 \right] $$where $y_t$, $t$, $x$, and $c$ represent the thickness, maximum thickness, $x$-coordinate and chord length, respectively. （Note: the last coefficient is modified from $-0.1015$ to $-0.1036$ in order to "close" the trailing edge. If we used $-0.1015$, there would be a thickness at the trailing edge, just like airfoils in the real world. But in the computational world, we prefer to use zero-thickness trailing edges.)

```
In [1]:
```import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

```
In [2]:
```def AFyt(x, t, c):
'''
input:
x: x coordinate of center line, float
t: maximum thickness, in fraction of chord length, float
c: chord lrngth, float
output:
half thickness of airfoil at corresponding x coordinate
'''
return 5. * t * (0.2969 * ((x/c)**0.5) -
0.126 * (x/c) -
0.3516 * ((x/c)**2) +
0.2843 * ((x/c)**3) -
0.1036 * ((x/c)**4))

** half thickness** of the airfoil, so the upper and lower profiles must be generated separately. The lower part takes the negetive values of those corresponding to the upper part. Using this function, we can generate the profiles for airfoils with the first and the second digits being zero (i.e., symmetrical airfoils). Here we show the profile of the NACA-0014 airfoil as an example:

```
In [3]:
```x = np.linspace(0., 1., 101)
plt.plot(x, AFyt(x, 0.14, 1.0), 'k-', lw=2) # upper surface
plt.plot(x, - AFyt(x, 0.14, 1.0), 'k-', lw=2) # lower surface
plt.axis('equal'); plt.ylim((-0.5, 0.5))

```
Out[3]:
```

For cambered 4-digit arifoils (i.e., the first and the second digit are not zero), we need to define Python functions for center lines (mean camber lines) and the angles, $\theta$, between center lines and the horizontal line.

$$ y_c = \left\{\begin{array}{lr} m\frac{x}{p^2}\left(2p-\frac{x}{c}\right)\text{,} && 0\le x\le p\times c\\ m\frac{c-x}{(1-p)^2}\left(1+\frac{x}{c}-2p\right)\text{,} && p\times c \le x \le c \end{array}\right. $$and

$$ \theta=\tan^{-1}{\frac{dy_c}{dx}} $$where

$$ \frac{dy_c}{dx}=\left\{\begin{array}{lr} \frac{2m}{p^2}\left(p-\frac{x}{c}\right)\text{,} && 0\le x\le p\times c\\ \frac{2m}{(1-p)^2}\left(p-\frac{x}{c}\right)\text{,} && p\times c \le x \le c \end{array}\right. $$```
In [4]:
```def AFyc(x, m, p, c):
'''
input:
x: x coordinate of center line, float
m: the maximum camber (100 m is the first of the four digits), float
p: location of maximum camber (10 p is the second digit), float
c: chord lrngth, float
output:
y coordinate of center line at corresponding x coordinate
'''
if (x >= 0.0) and (x <= p*c):
return m * x * (2. * p - (x/c)) / (p**2.)
elif (x > p*c) and (x <= c):
return m * (c - x) * (1. + (x/c) - 2. * p) / ((1. - p)**2)
else:
raise ValueError

```
In [5]:
```def AFth(x, m, p, c):
'''
input:
x: x coordinate of center line, float
m: the maximum camber (100 m is the first of the four digits), float
p: location of maximum camber (10 p is the second digit), float
c: chord lrngth, float
output:
angle between center and horizontal line at corresponding x coordinate
'''
if (x >= 0.0) and (x <= p*c):
return np.arctan(2.0 * m * (p - (x/c)) / (p**2))
elif (x > p*c) and (x <= c):
return np.arctan(2.0 * m * (p - (x/c)) / ((1. - p)**2))
else:
raise ValueError

Finally, the relationship between the coordinates on the upper and lower surfaces of cambered airfoils ($x_U$, $y_U$, $x_L$, and $y_L$), the center line coordinate ($x$ and $y_c$) and the thickness ($y_t$) are $$ \begin{array}{ll} x_U=x-y_t\sin{\theta}\text{,} && y_U=y_c+y_t\cos{\theta}\\ x_L=x+y_t\sin{\theta}\text{,} && y_L=y_c-y_t\cos{\theta} \end{array} $$

We wrote a function that can generate both symmetrical and cambered NACA 4-digit airfoils and use a sign to represent upper or lower profile:

```
In [6]:
```def AF(x, t, sign, m, p, c):
'''
input:
x: x coordinate of center line, float
t: maximum thickness, in fraction of chord length, float
sign: indicate upper (1) or lower (-1) surface of airfoil
m: the maximum camber (100 m is the first of the four digits), float
p: location of maximum camber (10 p is the second digit), float
c: chord lrngth, float
output:
x, y coordinates on airfoil surface at corresponding
center line x coordinate
'''
if (m == 0.) or (p == 0):
return x, sign * AFyt(x, t, c)
else:
return np.array([x[i] -
sign * AFyt(x[i], t, c) * np.sin(AFth(x[i], m, p, c))
for i in range(np.size(x))]), \
np.array([AFyc(x[i], m, p, c) +
sign * AFyt(x[i], t, c) * np.cos(AFth(x[i], m, p, c))
for i in range(np.size(x))])

The following shows the profile of a NACA 2412 airfoil as an example.

```
In [7]:
```x = np.linspace(0., 1., 101)
xU, yU = AF(x, 0.12, 1, 0.02, 0.4, 1.0)
xL, yL = AF(x, 0.12, -1, 0.02, 0.4, 1.0)
plt.plot(xU, yU, 'k-', lw=2)
plt.plot(xL, yL, 'k-', lw=2)
plt.axis('equal'); plt.ylim((-0.5, 0.5))

```
Out[7]:
```

Here we introduce the concepts of physical and computational domains, following chapter 9 of Hoffmann & Chiang, 2000 [2].

With the finite difference method, rectangular grids work best. But it is difficult to discretize the domain into rectangular grids when there are objects (like airfoils) within it.

In such a case, the domain can be discretized into a non-Cartesian structured mesh. Such a mesh may contain arbitrary quadrilateral grids, instead of rectangular grids. The approach is to map our real domain into another space that is discretized into rectangular grids and solve our problem in this space. The original domain is called *physical domain*, and the space in which we solve our problems is called *computational domain*. The following figure shows the concepts of physical and computational domains.

We define the coordinates in the computational domain to be $\xi$ and $\eta$ (called computational coordinates) and those in the physical doamin to be $x$ and $y$ (physical coordinates). Physical coordinates are functions of computational coordinates, and vice versa. This means that all calculations performed in the physical domain can also be performed in the computational domain through the relationship between $x$, $y$ and $\xi$, $\eta$, after their relationship has been determined.

If we fix the computational coordinates of grid points in a way that all grids in computational domain are rectangular grids with $\Delta\xi=\Delta\eta=1$, what we need to do next is to determine the physical coordinates of these computational grid points. Let $i$ and $j$ be the indices of grid points in the computational domain. What we want are the values $x_{i,j}, y_{i,j}$ representing the $x, y$ components in physical space of the grid points $(i,j)$ in the computational domain.

In this notebook, we use an O-type mesh on the physical domain, that is, the outer boundary of the physical domain is a circle, as shown in the figure below.

We divide the annular area between the arifoil and the outer boundary by a line $\overline{AC}$:

Then we "expand" the divided annular area to a rectangular domain, which will be our computational domain:

From the name of boundaries ($B_1$~$B_4$), we can understand the relationship between our physical domain and the computational domain.

```
In [8]:
```Nxi = 51
Neta = 21
eta, xi = np.meshgrid(np.linspace(0, Neta-1, Neta), np.linspace(0, Nxi-1, Nxi))

```
In [9]:
```def plotMesh(x, y):
for i in range(Nxi):
plt.plot(x[i, :], y[i, :], 'k.-', lw=2)
for i in range(Neta):
plt.plot(x[:, i], y[:, i], 'k.-', lw=2)

```
In [10]:
```plotMesh(xi, eta)

```
```

The idea of algebraic grid generation is to interpolate the physical coordinates of interior grid points between two known boundaries, i.e., the airfoil and the outer boundary in our case.

We know that the number of grid points on the airfoil and outer boundary should be the same with $N_{\xi}$, while on the dividing line, $\overline{AC}$, they should equal $N_{\eta}$. We also know the physical coordinates of grid points on $\eta=0$ and $\eta=N_{\eta}$, because they are located on the airfoil surface and the outer boundary and can be determined as follows:

$$ \begin{array}{ll} x_{i, j=0} = x_{airfoil}\text{,} && y_{i, j=0} = y_{airfoil}\\ x_{i, j=-1} = x_{outer\ boundary}\text{,} && y_{i, j=-1} = y_{outer\ boundary} \end{array} $$For convenience, we use Python indexing, where $j=-1$ represents the last grid points in the $\eta$ direction.

The $x$ and $y$ coordinates for interior points (i.e., $0\le i\le -1$ and $1\le j \le -2$, using Python indexing) can be obtained through interpolation. This is called algebraic grid-generation.

We now generate algebraic grids for a NACA 2412 airfoil. The ordering of grid points on the airfoil is clockwise and starts from the trailing edge. The center of the circular outer boundary is located at $0.5c$.

```
In [11]:
```rBC = 5.0 # the radius of outer boundary
m, p, t, c = 0.02, 0.4, 0.12, 1.0 # parameters of NACA 2412 airfoil
# Initialize x[i, j] and y[i, j]
x = np.empty((Nxi, Neta))
y = np.empty((Nxi, Neta))
# Generate grid points on airfoil surface
Nxc = (Nxi-1)/2 + 1
xc = np.linspace(0., 1., Nxc)
xU, yU = AF(xc, 0.12, 1, 0.02, 0.4, 1.0)
xL, yL = AF(xc, 0.12, -1, 0.02, 0.4, 1.0)
# Set x_{i, j=0} and y_{i, j=0}
x[:Nxc, 0] = xL[-1::-1].copy()
x[Nxc:, 0] = xU[1:].copy()
y[:Nxc, 0] = yL[-1::-1].copy()
y[Nxc:, 0] = yU[1:].copy()
# Generate grid points on circular outer boundary
# and set x_{i, j=-1}, y_{i, j=-1}
dr = 2. * np.pi / (Nxi -1)
th = - np.array([i * dr for i in range(Nxi)])
x[:, -1] = rBC * np.cos(th) + 0.5 * c
y[:, -1] = rBC * np.sin(th)

Now, we interpolate the coordinates of interior grid points.

```
In [12]:
```for i in range(Nxi):
x[i, 1:-1] = np.linspace(x[i, 0], x[i, -1], Neta)[1:-1]
y[i, 1:-1] = np.linspace(y[i, 0], y[i, -1], Neta)[1:-1]

Let's see what the mesh looks like.

```
In [13]:
```plt.figure(figsize=(8, 8), dpi=100)
plotMesh(x, y)
plt.axis('equal')
plt.xlim((-4.5, 5.5)); plt.ylim((-5, 5))

```
Out[13]:
```

Looking closer to the vicinity of the airfoil:

```
In [14]:
```plt.figure(figsize=(8, 8), dpi=100)
plotMesh(x, y)
plt.axis('equal')
plt.xlim((-0.5, 1.5)); plt.ylim((-1, 1))

```
Out[14]:
```

In the elliptic grid-generation method, we use elliptic PDEs to constrain the relationship between physical and computational coordinates. The elliptic PDEs used in this grid generation method are:

$$ \frac{\partial^2 \xi}{\partial x^2}+\frac{\partial^2 \xi}{\partial y^2}=0\\ $$$$ \frac{\partial^2 \eta}{\partial x^2}+\frac{\partial^2 \eta}{\partial y^2}=0\\ $$However, if we solve these PDEs (if we could), what we would obtain are the computational coordinates, i.e., $\xi=\xi(x, y)$ and $\eta=\eta(x, y)$, which are not what we want, because we already fixed the computational coordinates (they're simply on a rectangular grid with $\Delta\xi=\Delta\eta=1$).

Therefore, the actual PDEs that we are going to solve, which are derived from the above PDEs (see reference [2]), are

$$ a\frac{\partial^2 x}{\partial \xi^2}-2b\frac{\partial^2 x}{\partial \xi \partial \eta}+c\frac{\partial^2 x}{\partial \eta^2}=0 $$$$ a\frac{\partial^2 y}{\partial \xi^2}-2b\frac{\partial^2 y}{\partial \xi \partial \eta}+c\frac{\partial^2 y}{\partial \eta^2}=0 $$where

$$ a=\left(\frac{\partial x}{\partial \eta}\right)^2+\left(\frac{\partial y}{\partial \eta}\right)^2\\ b=\left(\frac{\partial x}{\partial \xi}\right)\left(\frac{\partial x}{\partial \eta}\right)+\left(\frac{\partial y}{\partial \xi}\right)\left(\frac{\partial y}{\partial \eta}\right)\\ c=\left(\frac{\partial x}{\partial \xi}\right)^2+\left(\frac{\partial y}{\partial \eta}\right)^2 $$We can see that, in order to solve the PDEs, we need two boundary conditions on each direction, This means that we need $x_{i,j=0}$, $x_{i,j=-1}$, $x_{i=0,j}$, $x_{i=-1,j}$, $y_{i,j=0}$, $y_{i,j=-1}$, $y_{i=0,j}$, and $y_{i=-1,j}$.

Luckly, $x_{i,j=0}$, $x_{i,j=-1}$, $y_{i,j=0}$, and $y_{i,j=-1}$ are known and are the physical coordinates of nodes on the airfoils and outer boundary, as described in the previous subsection.

From the illustration above, we also know that $x_{i=0,j}$ and $x_{i=-1,j}$ are actually the same, because they both represent the physical coordinates of nodes on the dividing line $\overline{AC}$. Similarly, $y_{i=0,j}$ and $y_{i=-1,j}$ are the same, too. This implies that ** periodic boundary conditions** can be applied on $x_{i=0,j}$, $x_{i=-1,j}$, $y_{i=0,j}$, and $y_{i=-1,j}$.

These PDEs are non-linear but we linearize the solution by letting the coefficients lag behind by one iteration. We'll use the basic Jacobi iterative method with the results obtained from the algebraic grid-generation method as initial guess.

We can use central differences on these PDEs, and with $\Delta\xi=\Delta\eta=1$, the discretized equations are very simple. Since the forms of these two PDEs are the same, their discretized PDEs are also the same. The discretized equation of $x$ are

$$ a(x_{i+1, j}-2x_{i,j}+x_{i-1, j}) - \frac{1}{2}b(x_{i+1, j+1}-x_{i+1, j-1}+x_{i-1,j-1}-x_{i-1,j+1})+c(x_{i, j+1}-2x_{i,j}+x_{i, j-1})=0 $$or

$$ x_{i,j}=\frac{1}{2}\left\{ a(x_{i+1, j}+x_{i-1, j}) - \frac{1}{2}b(x_{i+1, j+1}-x_{i+1, j-1}+x_{i-1,j-1}-x_{i-1,j+1})+c(x_{i, j+1}+x_{i, j-1}) \right\} / (a+c) $$where

$$ a=(\frac{x_{i,j+1}-x_{i,j-1}}{2})^2+(\frac{y_{i,j+1}-y_{i,j-1}}{2})^2\\ b=(\frac{x_{i+1,j}-x_{i-1,j}}{2})(\frac{x_{i,j+1}-x_{i,j-1}}{2})+(\frac{y_{i+1,j}-y_{i-1,j}}{2})(\frac{y_{i,j+1}-y_{i,j-1}}{2})\\ c=(\frac{x_{i+1,j}-x_{i-1,j}}{2})^2+(\frac{y_{i+1,j}-y_{i-1,j}}{2})^2 $$Replace $x$ with $y$ when we want to solve $y$.

We can now solve these equations. The criterion for stopping the iteration is that the maximum difference between a current result and the result of the previous iteration should be less than $10^{-6}$.

```
In [15]:
```def Solve_a_b_c(x, y):
'''
input:
x: the x coordinate of x_{i=i-1~i+1, j=j-1~j+1}, at least 3x3 array
y: the y coordinate of y_{i=i-1~i+1, j=j-1~j+1}, at least 3x3 array
output:
a, b, c: at least 1x1 float
'''
a = 0.25 * (((x[1:-1, 2:] - x[1:-1, :-2])**2) +
((y[1:-1, 2:] - y[1:-1, :-2])**2))
b = 0.25 * ((x[2:, 1:-1] - x[:-2, 1:-1]) *
(x[1:-1, 2:] - x[1:-1, :-2]) +
(y[2:, 1:-1] - y[:-2, 1:-1]) *
(y[1:-1, 2:] - y[1:-1, :-2]))
c = 0.25 * (((x[2:, 1:-1] - x[:-2, 1:-1])**2) +
((y[2:, 1:-1] - y[:-2, 1:-1])**2))
return a, b, c

```
In [16]:
```def SolveEq(a, b, c, U):
'''
input:
a, b, c: as described in the content
U: the result of the last iteration
output:
return the result of current iteration
'''
return 0.5 * (
a * (U[2:, 1:-1] + U[:-2, 1:-1]) +
c * (U[1:-1, 2:] + U[1:-1, :-2]) -
b * 0.5 * (U[2:, 2:] - U[2:, :-2] + U[:-2, :-2] - U[:-2, 2:])
) / (a + c)

```
In [17]:
```iters=0
while True:
# count the number of iterations
iters += 1
# backup the last result
xn = x.copy()
yn = y.copy()
# solve periodic BC first
tempx = np.append([x[-2, :].copy()], x[0:2, :].copy(), 0)
tempy = np.append([y[-2, :].copy()], y[0:2, :].copy(), 0)
a, b, c = Solve_a_b_c(tempx, tempy)
x[0, 1:-1] = SolveEq(a, b, c, tempx)
y[0, 1:-1] = SolveEq(a, b, c, tempy)
x[-1, 1:-1] = x[0, 1:-1].copy()
y[-1, 1:-1] = y[0, 1:-1].copy()
# solve interior
a, b, c = Solve_a_b_c(x, y)
x[1:-1, 1:-1] = SolveEq(a, b, c, x)
y[1:-1, 1:-1] = SolveEq(a, b, c, y)
# calculate difference between current and the last result
errx = np.abs(x - xn)
erry = np.abs(y - yn)
# adjudge whether the iteration should stop
if (errx.max() <= 1e-6) and (erry.max() <= 1e-6):
break

Let's see what the mesh looks like.

```
In [18]:
```plt.figure(figsize=(8, 8), dpi=100)
plotMesh(x, y)
plt.axis('equal')
plt.xlim((-4.5, 5.5)); plt.ylim((-5, 5))

```
Out[18]:
```

```
In [19]:
```plt.figure(figsize=(8, 8), dpi=100)
plotMesh(x, y)
plt.axis('equal')
plt.xlim((-0.5, 1.5)); plt.ylim((-1, 1))

```
Out[19]:
```

We know the governing equation of the stream function $\Psi$ in the physical domain is the Laplace equation, as described in the introduction section. To solve the stream function in the computational domain, we need a new governing equation. According to the derivation in reference [3], we obtain the governing equation for stream function in the computational domain as follows:

$$ a\frac{\partial^2 \Psi}{\partial \xi^2}-2b\frac{\partial^2 \Psi}{\partial \xi \partial \eta}+c\frac{\partial^2 \Psi}{\partial \eta^2}=0 $$where

$$ a=\left(\frac{\partial x}{\partial \eta}\right)^2+\left(\frac{\partial y}{\partial \eta}\right)^2\\ b=\left(\frac{\partial x}{\partial \xi}\right)\left(\frac{\partial x}{\partial \eta}\right)+\left(\frac{\partial y}{\partial \xi}\right)\left(\frac{\partial y}{\partial \eta}\right)\\ c=\left(\frac{\partial x}{\partial \xi}\right)^2+\left(\frac{\partial y}{\partial \eta}\right)^2 $$The form of the governing equation is the same as we solved in the elliptic grid-generation method! This means we can use the functions we wrote previously for that purpose.

Though the form of the governing equation is the same as the PDEs for the elliptic grid-generation method, the boundary conditions are not. We now examine the boundary conditions.

On the outer boundary, we expect the flow velocity to be equal to the free-stream velocity, that is:

$$ \left\{ \begin{array}{l} u=\frac{\partial \Psi}{\partial y}=V_\infty\cos(AOA)\\ v=-\frac{\partial \Psi}{\partial x}=V_\infty\sin(AOA) \end{array}\right.,\ on\ outer\ boundary $$where $AOA$ represents angle of attact, and $V_\infty$ is the free stream velocity.

Since $V_\infty$ and $AOA$ are constant, we can simply integrate the above boundary conditions for the stream function on the outer boundary

$$ \Psi=-xV_\infty\sin(AOA)+yV_\infty\cos(AOA)+C_1,\ on\ outer\ boundary $$In potential flow, the values of $\Psi$ are not important; what matters are the derivatives of the stream function. Since the constant $C_1$ in the above equation will disappear in the derivatives, we can simply make $C_1=0$ here. Now we can apply Dirichlet boundary condition on the outer boundary:

$$ \Psi_{i,j=-1}=-x_{i,j=-1}V_\infty\sin(AOA)+y_{i,j=-1}V_\infty\cos(AOA) $$Next, we examine the boundary conditions for nodes on the airfoil. For walls in potential flows, the no-slip boundary condition does not apply, given the lack of viscosity. The boundary condition for walls in potential flows is that the flow can not penetrate the wall (no-thru BC). This implies that the wall itself is a streamline.

$$ \Psi=C_2,\ on\ the\ airfoil $$or

$$ \Psi_{i,j=0}=C_2 $$Since we have determined the values of $\Psi$ on the outer boundary, we cannot determine the value of $\Psi$ on the airfoil, $C_2$ (because we don't know which streamline on the outer boundary corresponds to the streamline on the airfoil surface). Fortunately, we will resolve this problem later after we introduce the Kutta condition.

Finally, the boundary condition on the dividing line $\overline{AC}$ is again given from periodic boundary condition.

$$ \Psi_{i=0,j}=\Psi_{i=-1,j} $$In potential flow, due to the lack of viscosity, the flow under an airfoil will flow back to the upper surface after it has passed the trailing edge, even at a low angle of attact and low free stream velocity. This is not physical. The Kutta condition ensures this situation won't happen: it makes the flow smoothly pass by the trailing edge. More details can be found in any textbook covering inviscid flow or potential flow. Here we only introduce the implementation of the Kutta condition in our code.

According to reference [3], for the airfoils with finite-angle trailing edge, like the NACA 2412, the Kutta condition results in the value of $\Psi$ on the node next to the trailing edge being the same as that on the trailing edge, in order to provide a smooth streamline near the trailing edge. In other words,

$$ \Psi_{i=0, j=1}=\Psi_{i=0, j=0} $$Recall the boundary condition on the airfoil surface: $\Psi_{i,j=0}=C_2$. We can determine $C_2$ and $\Psi_{i,j=0}$ now. In each iteration when solving the governing equation, let

$$ \Psi_{i,j=0}=C_2=\Psi_{i=0,j=1} $$The problem of the boundary condition on the airfoil surface is then resolved. We can start to solve for the stream functions!

```
In [20]:
```# set angle of attact
AOA = 15. / 180. * np.pi
# set free stream velocity
Vinf = 70.
# initialize stram functions
stream = np.zeros((Nxi, Neta))
# set up the BCs on the outer boundary
stream[:, -1] = - x[:, -1] * Vinf * np.sin(AOA) + y[:, -1] * Vinf * np.cos(AOA)
# solve the PDE by iterative method
iters = 0
while True:
# count the number of current interation
iters += 1
# backup the last result
stream_n = stream.copy()
# apply periodic BC on dividing line
temp = np.append([stream[-2, :].copy()], stream[:2, :].copy(), 0)
tempx = np.append([x[-2, :].copy()], x[:2, :].copy(), 0)
tempy = np.append([y[-2, :].copy()], y[:2, :].copy(), 0)
a, b, c = Solve_a_b_c(tempx, tempy)
stream[0, 1:-1] = SolveEq(a, b, c, temp)
stream[-1, :] = stream[0, :].copy()
# apply Kutta condition
# and set the value of stream function on the airfoil surface
stream[:, 0] = stream[0, 1]
# solve interior
a, b, c = Solve_a_b_c(x, y)
stream[1:-1, 1:-1] = SolveEq(a, b, c, stream)
# calculate difference between current and the last result
err = np.abs(stream - stream_n)
# adjudge whether the iteration should stop
if (err.max() <= 1e-6):
break

```
In [21]:
```stream = stream - stream[0, 0]

Finally, the streamlines are shown in the following figure.

```
In [22]:
```# set the contour lines with negative values to be solid lines
import matplotlib as mpl
mpl.rcParams['contour.negative_linestyle'] = 'solid'
# contour
plt.figure(figsize=(10, 8), dpi=100)
cs = plt.contour(x, y, stream, 100, colors='k')
plt.clabel(cs)
plt.plot(x[:, 0], y[:, 0], 'k-', lw=1) # plot the airfoil
plt.xlim((-2., 3.))
plt.ylim((-1.5, 1.5))

```
Out[22]:
```

[1] *NACA airfoil*, Wikipedia, http://en.wikipedia.org/wiki/NACA_airfoil

[2] K.A. Hoffmann, and S.T. Chiang, *Computational fluid dynamics, Vol. 1*, Wichita, KS: Engineering Education System (2000).

[3] F. Mohebbi and M. Sellier. *On the Kutta Condition in Potential Flow over Airfoil*, Journal of Aerodynamics (2014).

```
In [23]:
```from IPython.core.display import HTML
css_file = '../../styles/numericalmoocstyle.css'
HTML(open(css_file, "r").read())

```
Out[23]:
```