AERO40003 Computing & Numerical Methods 1

Part 1:  MATLAB

The aim of this coursework is for you to write a few more complex programs than what you typically have encountered so far in tutorials.   Unlike most tutorial questions, you will have to come up with the structure of the program yourself.  This is what programming will involve for the rest of your career after all.

The following brief presents the background information and instructions necessary for you to complete this part of the coursework.  You must submit all the files (the  .m files for these scripts/functions and figures) through this submission box in a  .zip file.

The coursework will be marked based on completeness, readability, ele- gance, consistency, correctness and obviously on whether the programs run!

Panel Codes

1    Theoretical Background

The flow field around an arbitrary lifting body in potential flow (incom- pressible, inviscid and irrotational) can be predicted computationally using a Panel Code. Panel codes work by breaking up the surface of an arbitrary body into a number of panels, as seen in figure 1.

In order to apply a Panel Code to an airfoil, first we need to discretise its surface into N panels.  The NACA 4-series of airfoils is a well known family of sections that can be generated using the parameters defined by the digits of their numerical code.

The first digit represents the maximum camber (m), in percent of chord; the second digit indicates the location of maximum camber (p), in tenths of chord; and the final two digits correspond to the maximum aerofoil thickness (t), which for this family occurs at 30% of the chord from the leading edge.

Figure 1: Dicretization of a NACA 2415 airfoil using 10 panels on its surface.

The numbers shown correspond to the order in which the panel endpoints are generated and stored in the arrays x and z of the function panelgen().

The NACA 2415 example would therefore have a maximum of 2% camber at 0.4 of the chord, and a maximum thickness of 15% of the chord. Symmetrical airfoils will start with 00, as they have no camber.

Aerofoil sections can therefore be constructed by adding a thickness func- tion over a mean camber line.  The mean camber line for this family of aero- foils can then be obtained from the following expression (note, x is the non- dimensional coordinate, 0 at the leading edge and 1 at the trailing edge):

y卜 =  ,．  × r) ),                                0 5 r 5 s,

．．(1 × s))((1 × 2s) + 2sr × r) ),    s 5 r 5 1

Similarly, the thickness function for this family is given by:

y≠ = 5u(0.2969,r × 0.126r × 0.3516r) + 0.2843r4 × 0.1015r= )

The thickness needs to be apply perpendicular to the camber line, therefore the coordinates of any point on the airfoil surface become:

r扌 = r × y≠sin 9,

rf = r + y≠sin 9,

where 9 is given by:

dy卜

9 = arctan

with

The panel end point r can be calculated using a cosine distribution, which gives the first value at the trailing edge:

r = 1 × 0.5 ╱ 1 × cos 2m、

An additional panel will be placed starting at the trailing edge and extend- ing to infinity (some large number) in the direction of the flow to represent the wake of the airfoil.

Each 2-dimensional panel can be made up of sources, sinks, doublets or vortices, hence affecting the flowfield around it. In this case we will be using doublet elements, each of which will affect the flowfield around it by imparting a horizontal velocity v and a vertical velocity z, proportional to the element’s strength u. Panel codes work by finding what strength each panel around the body should have so that a certain condition is met. The condition we will be using, is that the velocity normal to the wall and therefore also at the centre of each panel must be zero (non-permeability).

Therefore if we can find the total velocities U;  and v;  imparted by the freestream flow and all the panels at the centre of the the i≠A  panel, we can ensure that the total velocity normal to it is zero.  Mathematically the total velocity at the centre of the i≠A  panel can be expressed as:

衣＋1

U; =U○ cos a +       u〉v;〉

〉二1

衣＋1

v; =U○ sin a +       u〉z;〉

〉二1

where U○  is the speed of the freestream airflow, a is the airfoil’s angle of attack, u〉  is the strength of the j≠A panel and N is the total number of panels.

The velocities imparted at the centre of the i≠A panel by the j≠A panel, v;〉 and z;〉, can be found using the function cdoublet().

We can also find that the velocity component normal to the i≠A  panel at its midpoint is:

(U≠); = × U;sin 8; + v;cos 8;                                               (9)

where 8;  is the panel angle given by

8; = tan ← 1  ╱ 、 = tan ← 1  ╱ 、

where [r;, z;] and [r;＋1, z;＋1] are the coordinates of the i≠A  panel’s endpoints. Note that 8;  should be estimated in such a way that 0 5 8;  5 2m, depending on the sign of dr and dz.(you may wish to use the atan2() function such that 0 5 8;  5 2m).

Combining equations (8) and (9), a system of N + 1 linear equations, with N + 1 unknowns is created. The first N relations are given by:

衣＋1

(U≠); = U○ sin(a × 8;) +       u〉(z;〉cos 8; × v;〉sin 8;) = 0         (11)

〉二1

and therefore

衣＋1

u〉(z;〉cos 8; × v;〉sin 8;) = × U○ sin(a × 8;)                 (12) 〉二1

The final relation comes from the Kutta condition and simply states that circulation at the trailing edge must be zero and therefore:

u衣＋1 + u1 × u衣 = 0                                        (13)

The system of equations can be solved computationally by posing the prob- lem in the form:

Ax = B                                                  (14)

where A is an (N +1) · (N +1) square matrix and x and B are column vectors of length (N + 1).

Once the system has been solved and the unknown panel strengths (u1 to u衣＋1) have been found, equation (8) can be used to estimate what the velocity field around the airfoil is. The lift coefficient of the airfoil can also be crudely estimated using the strength of the wake panel doublet:

Cf = ×

2    Implementation

First write a function panelgen() that takes the airfoil code as a string (i.e.   ’1410’), the number of panels to be used N  and the angle of attack ApA.  The function should use the method described above to calculate the coordinates of all the panel end points.  The function should return two one dimensional arrays r and z .

Write a script that asks the user which airfoil to use, the freestream velocity U○ , the angle of attack a and the number of panels to be used N . The code should:

● Use the function panelgen() to generate and discretize the airfoil.

●  Generate the matrices A and B and solve the system of equations to find the unknown panel strengths.  You may find that the matrix A is close to singular, but that is common for this kind of problem.

●  Output the airfoil’s lift coefficient.

●  Plot the velocity field vectors (arrows) and streamlines around the airfoil in the domain where ×0.2 5 r 5 1.2 and ×0.7 5 z 5 0.7.  To achieve that you will have to evaluate the velocities U and v not on the airfoil body but at various points in the domain surrounding the airfoil. Note that you will need a much finer grid to accurately draw the streamlines than to draw vectors. The inpolygon() function may also prove useful for identifying which points are inside your airfoil and should thus not be plotted.

●  Save the plots programmatically.

However, if the airfoil is NACA  2412, the script should:

●  Calculate the lift coefficient for a range of angles of attack between 0= and 10=  and for different values of N : 50, 100 and 200.

●  Read the XFOIL data file for this airfoil provided.

●  Plot the Cf  vs.  a for the different cases of N against the XFOIL data for the range of angles calculated.

●  Plot the velocity field vectors and streamlines for the last case.

●  Save your plots programmatically.

3    Sample case

Try your program for different NACA 4-series airfoils, at various angles of attack.  Check your results against the theoretical lift curve slope of 2m for those airfoils and the flowfields seen in figures 2 and 3.

0.6

0.4

0.2

0

−0.2

−0.4

−0.6

−0.2          0          0.2        0.4        0.6        0.8          1          1.2

Figure 2: Streamlines around a NACA 2405 at a = 5=  and using N = 250.

4    What to submit

You should submit the following items in a folder called CW_Pt1_username, replacing username with your own, and zip the folder before uploading:

● Your panelgen() function.

0.6

0.4

0.2

0

−0.2

−0.4

−0.6

−0.2          0          0.2        0.4        0.6        0.8          1          1.2

Figure 3: Streamlines around a NACA 2117 at a = 10=  and using N = 250.

● Your script running the panel code you have written.

● Arrow and streamline plots for two different NACA 4-series airfoils of your choice, each at two different angles of attack and at the speed and discretization of your choice (but label accordingly).

●  Cf vs. a plot for NACA 2412 and arrow and streamline plots at 15m/t, a = 10=  and N = 200.

5    Formatting guidelines

As you have more time for this project, it is expected that you will put some care into presenting both your code and plots neatly and professionally. Please refer to the mark scheme to make sure your code complies with the ex- pectations we will be looking when marking. Correctness is indeed important. But other criteria such as readability or elegance are also important.

For this exercise, it is recommended that you work on an .m file.  Save the file with some meaningful name.   Parts of your code can be provided in functions, which should also be named meaningfully so the markers can identify them.  Functions might make your code more elegant and readable, but only one is required.

Additionally, the following criteria are expected for the formatting of the graphs:

● Label your axes and add titles to your plots.

● For line plots with multiple lines, add a legend.