PANEL METHOD CODE

List of Contents

1. Introduction
2. Mathematical Background
3. Code Structure
   1. Panel Method Module
   2. Streamlit App
   3. External Packages
4. Next Steps
5. References

---

Introduction

This code is a panel method solver for 2D airfoils and multi airfoil wings. It is based on the potential flow theory and uses a source and vortex distribution to model the flow around the airfoil. The code is composed the panel_method module, which contains the main functions for the panel method solver and a streamlit app.

Mathematical Background

The method assumes a potential flow. The flow is modeled using a combination of source and vortex distributions along the surface of the airfoil. The goal of the solver is to determine the strength of these distributions such that the following boundary conditions are satisfied:

1. The flow tangency condition: The flow velocity at each panel must be tangent to the surface of the airfoil.
2. The Kutta condition: The tangential velocity at the trailing edge's upper and lower panel must be zero.

To achieve this, a system of linear equations is formulated based on the influence of each panel on every other panel. To assemble this system, it is, first, necessary to calculate some geometrical parameters for the whole wing:

$$ n_{te} = \text{Number of trailing edges (flaps)} $$ $$ n_{pts} = \text{Total number of points} $$ $$ n_{pan} = \text{Total number of panels} $$ $$ n_{pts, list} = \text{Number of points in each airfoil} $$ $$ n_{pan, list} = \text{Number of panels in each airfoil} $$ $$n_{kutta} = \text{Indices of the panels at the trailing edges for each airfoil} $$

In the code the following index are used:

• $i$: Index of the control panel
• $j$: Index of the panel influenceing the control panel
• $k$: Index of x and y coordinates (0 for x and 1 for y)

The basic input is a list of points that define the airfoil shape. In the code, this information is represented by the following np.array object:

$$ X_{ik} = \begin{bmatrix} x_1 & y_1 \\ x_2 & y_2 \\ \vdots & \vdots \\ x_{n_{pts}} & y_{n_{pts}} \end{bmatrix} = \begin{bmatrix} X_{11} & X_{12} \\ X_{21} & X_{22} \\ \vdots & \vdots \\ X_{n_{pts}1} & X_{n_{pts}2} \end{bmatrix} $$

Then, the control point of each panel is calculated as the midpoint between the two endpoints:

$$ X^c_{ik} = \frac{1}{2} \cdot \left( X_{ik} + X_{(i+1)k} \right) $$

The length of each panel is calculated as:

$$ dS_i = \sqrt{(x_{i+1} - x_i)^2 + (y_{i+1} - y_i)^2} $$

Then, the tangent and normal vectors of each panel are calculated as:

$$ t_{ik} = \left( X_{(i+1)k} - X_{ik} \right) / dS_i $$ $$ n_{ik} = -\vec{t}_{i0} \cdot \vec{e_x} + \vec{t}_{i1} \cdot \vec{e_y} $$

Where $\vec{e}_x$ and $\vec{e}_y$ are the unit vectors in the x and y directions, respectively.

The previous calculations are made with respect to the global coordinate system. However, to calculate the influence coefficients, it is easier to transform the coordinates to a local system aligned with each panel. This can be done by using the normal and/or tangent vectors to define the transformation (the $'$ symbol represents the local coordinate system). To do this, first a tensor $\bar{\bar{\bar{R}}}$ where:

$$ R_{ijk} = X^c_{ik} - X_{jk} $$

Tis tensor stores the vector from the j-th panel to the i-th control point. Then, the local coordinates are calculated as:

$$ X'^c_{ij0} = x'^c_{ij} = R_{ijk} \cdot t_{jk} $$ $$ X'^c_{ij1} = y'^c_{ij} = R_{ijk} \cdot n_{jk} $$

With all this information, it is possible to calculate the influence coefficients of each panel on every other panel. Two different influences are considered: the source influence and the vortex influence. Both of them impose a $\Delta \vec{V} = \Delta u \cdot \vec{e}_x + \Delta v \cdot \vec{e}_y$ on the control point of the i-th panel due to the j-th panel. As said before, the influence is first computed with respect to the local coordinate system. Therefore, the source influence is calculated as:

$$ \Delta u'_s = \lambda_j \cdot \frac{1}{4\pi} \cdot \ln{\left[ \frac{(x'^c_{ij} + \frac{1}{2}dS_j)^2 + (y'^c_{ij})^2}{(x'^c_{ij} - \frac{1}{2}dS_j)^2 + (y'^c_{ij})^2} \right]} $$

$$ \Delta v'_s = \lambda_j \cdot \frac{1}{2\pi} \cdot \left[ \arctan{\left( \frac{x'^c_{ij} + \frac{1}{2}dS_j}{y'^c_{ij}}\right)} - \arctan{\left( \frac{x'^c_{ij} - \frac{1}{2}dS_j}{y'^c_{ij}}\right)} \right] $$

Analogously, the vortex influence can be calculated as:

$$ \Delta u'_v = \gamma_j \cdot \frac{1}{2\pi} \cdot \left[ \arctan{\left( \frac{x'^c_{ij} + \frac{1}{2}dS_j}{y'^c_{ij}}\right)} - \arctan{\left( \frac{x'^c_{ij} - \frac{1}{2}dS_j}{y'^c_{ij}}\right)} \right] $$

$$ \Delta v'_v = \gamma_j \cdot \frac{-1}{4\pi} \cdot \ln{\left[ \frac{(x'^c_{ij} + \frac{1}{2}dS_j)^2 + (y'^c_{ij})^2}{(x'^c_{ij} - \frac{1}{2}dS_j)^2 + (y'^c_{ij})^2} \right]} $$

Now, it is necessary to transform the influence back to the global coordinate system. This can be done using the following transformations:

$$ \Delta V_{ijk} = \Delta u'_{ij} \cdot t_{jk} + \Delta v'_{ij} \cdot n_{jk} $$

Now, that all the influence coefficients are calculated, it is possible to assemble the system of linear equations:

$$ [A] \cdot \begin{bmatrix} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_{n_{pan}} \\ \gamma_1 \\ \vdots \\ \gamma_{n_{te}} \end{bmatrix} = \{b\} $$

That means that we have to solve for $n_{pan} + n_{te}$ unknowns (the source strengths and the vortex strengths at the trailing edges). The first $n_{pan}$ equations are obtained by applying the flow tangency condition at each control point. Meanwhile, the last $n_{te}$ equations are obtained by applying the Kutta condition at each trailing edge.

The matrix $[A]$ and vector ${b}$ can be build as:

$$ [A] = \begin{bmatrix} \begin{array}{c|c} [A_{11}]_{n_{pan} \times n_{pan}} & [A_{12}]_{n_{pan} \times n_{te}} \\\ \\ \hline \\\ [A_{21}]_{n_{te} \times n_{pan}} & [A_{22}]_{n_{te} \times n_{te}} \end{array} \end{bmatrix} \qquad and \qquad \{b\} = \begin{Bmatrix} \{b_1\}_{n_{pan} \times 1} \\\ \\ \hline \\\ \{b_2\}_{n_{te} \times 1} \end{Bmatrix} $$

Therefore, each of the sub-matrices and sub-vectors can be build as:

A11 = np.einsum('ik,ijk->ij', N, dVs)
A12 = np.add.reduceat(np.einsum('ik,ijk->ij', N, dVv), ids_reduce, axis=1)
A21 = np.einsum('ik,ijk->ij', T[n_kutta[:,0], :], dVs[n_kutta[:,0], :, :]) + np.einsum('ik,ijk->ij', T[n_kutta[:,1], :], dVs[n_kutta[:,1], :, :])
A22 = np.add.reduceat(np.einsum('ik,ijk->ij', T[n_kutta[:,0], :], dVv[n_kutta[:,0], :, :]) + np.einsum('ik,ijk->ij', T[n_kutta[:,1], :], dVv[n_kutta[:,1], :, :]), ids_reduce, axis=1)

And:

b1 = -Vinf.dot(N.T).T
b2 = -(Vinf.dot(T[n_kutta[:,0], :].T) + Vinf.dot(T[n_kutta[:,1], :].T)).T

Code Structure

Panel Method module

Setup

This code is responsible to provide the class used as input to the solver which allow us to simulate.

Geometry

This code is responsable to generate the wing. It works with 2 base base class: Airfoil e Wing. The first is responsible to adjust the flap so that the solver can interpret it without problems. That means that it has to do a few operations besides the collection, translation and rotation of points. Those operations can be seen in the code itself. The Wing class by the other side is responsible to collect all the airfoils and create a list with all of them and a few other things...

Solver

This is the run code. In it you can see the solver for the metho and the class responsable for holding the results provided by tghe functions

Post Processing

This code is resposable for the analysis of the results previously obtained. It is runly resposable for the contour and grid generation so that we can have a better visualization of the flow.

Report Generation

This code is responsible to make the results visualization better in order to understand the results in a more efficient way.

Streamlit App

The streamlit app is composed by the main file Home.py and a few pages in the pages folder. This main file is responsible to provide the general layout of the app and show this README, and the pages are responsible to provide the functionalities for the simulation.

run streamlit Home.py

External packages

All the packages used in this code are listed in the requirements.txt file. To install them just run the following command in your terminal:

pip install requirements.txt

Next steps

Althought this code is functional, it still have a long way to go in terms of performance and funcionalities. In this section I list a few od the next few steps I believe can be done in order to make this code better.

1. Add more pre biuild funcionalities. Add more functions and improve the app. Including the Profile Editor page to edit the points file.

2. Validate the results comparing then to ones obtained on xFoil. I left this for last because it should be the first one I remove.

3. Add turbulent viscous flow model to the code. This will allow the results to be much closer to real ones and use higher angle of attacks with more fidelity on the results at a larger range.

References

Below are some useful references used during the development of this code:

1. Katz, Plotkin: Low Speed Aerodynamics
2. Oregon State - Intermediate Fluid Mechanics (specially chapter 4, 5 and 6)
3. Science Direct - Panel Method Overview