Edit on GitHub

Objectives

Here we will show you how to perform an aerodynamic analysis of an airfoil.

Recipe

  1. Compute the coordinates of a NACA 4-digit series airfoil.
  2. Plot its camber, thickness, upper and lower surface representations.
  3. Perform an aerodynamic analysis at a given angle of attack.
  4. Plot its aerodynamic properties.

For this, we will need to import some packages which will be convenient for plotting.

using AeroFuse      # Main package
using Plots         # Plotting library
gr(dpi = 300)       # Plotting backend
using LaTeXStrings  # For LaTeX printing in plots

Your First Airfoil

You can define a NACA 4-digit series airfoil using the following function.

airfoil = naca4((2,4,1,2), 60)
NACA 2412 Foil{Float64} with 119 points.

You can access the $x$- and $y$-coordinates as the following fields.

airfoil.x, airfoil.y
([1.0, 0.9992981882625027, 0.9971946084073808, 0.9936948246213843, 0.9888081041263779, 0.9825474070711719, 0.9749293714656669, 0.9659742923191958, 0.9557060940082984, 0.9441522948292751  …  0.9431997802272964, 0.9549286788466147, 0.9653618851331425, 0.9744682370027143, 0.9822204617433441, 0.9885952856914888, 0.9935735300999603, 0.9971401918027563, 0.9992845074142591, 1.0], [-1.6616460838224803e-17, 0.00014994327027663046, 0.0005984546661103868, 0.0013416010129608839, 0.002372899768779709, 0.0036834253181048772, 0.005261952261861909, 0.007095130492785934, 0.009167685786814105, 0.011462638875020038  …  -0.004305264850848705, -0.0034318413055497835, -0.002648506782021171, -0.001959587327218415, -0.0013690972131801405, -0.00088064261196314, -0.0004973272941477818, -0.00022166619369339412, -5.551211404810675e-5, 1.6616460838224803e-17])

You can obtain the coordinates as an array by calling the following function.

coordinates(airfoil)
119×2 Matrix{Float64}:
 1.0       -1.66165e-17
 0.999298   0.000149943
 0.997195   0.000598455
 0.993695   0.0013416
 0.988808   0.0023729
 0.982547   0.00368343
 0.974929   0.00526195
 0.965974   0.00709513
 0.955706   0.00916769
 0.944152   0.0114626
 ⋮         
 0.954929  -0.00343184
 0.965362  -0.00264851
 0.974468  -0.00195959
 0.98222   -0.0013691
 0.988595  -0.000880643
 0.993574  -0.000497327
 0.99714   -0.000221666
 0.999285  -5.55121e-5
 1.0        1.66165e-17

Geometric Representations

You can convert these coordinates into the camber-thickness representation.

xcamthick = camber_thickness(airfoil, 60)
61×3 Matrix{Float64}:
 0.0          0.0           0.0
 0.000685233  0.00165715    0.00867606
 0.00273905   0.0017672     0.0180042
 0.00615583   0.00203342    0.0269113
 0.0109262    0.00242718    0.0355243
 0.0170371    0.00293959    0.0438502
 0.0244717    0.00356255    0.0518745
 0.0332098    0.00428708    0.0595753
 0.0432273    0.00510318    0.066927
 0.0544967    0.0059998     0.0739022
 ⋮                         
 0.956773     0.0028295     0.0122458
 0.96679      0.00219383    0.00946855
 0.975528     0.00162951    0.0070167
 0.982963     0.00114212    0.00490866
 0.989074     0.000736515   0.00316062
 0.993844     0.000416739   0.0017863
 0.997261     0.000186      0.00079662
 0.999315     4.66184e-5    0.000199565
 1.0          0.0          -3.32329e-17

You can split the coordinates into their upper and lower surfaces.

upper, lower = split_surface(airfoil);

x_upper, y_upper = @views upper[:,1], upper[:,2]
x_lower, y_lower = @views lower[:,1], lower[:,2];

Plotting

You can plot the airfoil by calling plot from Plots.jl and passing the Foil object as the first argument. Optionally, you can enable flags to plot the camber and thickness distributions.

plot(
    airfoil,
    camber = true,
    thickness = true,
    aspect_ratio = 1,
    xlabel=L"(x/c)",
    ylabel = L"y"
)

Your First Doublet-Source Analysis

Now we have an airfoil, and we would like to analyze its aerodynamic characteristics. The potential flow panel method for inviscid analyses of airfoils, which you may have studied in your course on aerodynamics, provides decent estimations of the lift generated by the airfoil.

Our analysis also requires boundary conditions, which is the freestream flow defined by a magnitude $V_\infty$ and angle of attack $\alpha$. We provide these to the analysis by defining variables and feeding them to a Uniform2D type, corresponding to uniform flow in 2 dimensions.

V       = 1.0
alpha   = 4.0 # degrees
uniform = Uniform2D(V, alpha)
Uniform2D{Float64}(1.0, 0.06981317007977318)

You can analyze this airfoil for the given freestream flow as shown.

system  = solve_case(
                     airfoil, uniform;
                     num_panels = 80
                    );

This will run the analysis and return a system which can be used to obtain the aerodynamic quantities of interest and post-processing. The panels used for the analysis can be accessed as follows.

panels = system.surface_panels
78-element view(::Vector{Panel2D{Float64}}, 78:-1:1) with eltype Panel2D{Float64}:
 Panel2D{Float64}([1.0, 1.6616460838224803e-17], [0.9983786540671049, -0.00012570291449866332])
 Panel2D{Float64}([0.9983786540671049, -0.00012570291449866332], [0.9935251313189564, -0.0005010539079096189])
 Panel2D{Float64}([0.9935251313189564, -0.0005010539079096189], [0.985470908713026, -0.0011200400170492484])
 Panel2D{Float64}([0.985470908713026, -0.0011200400170492484], [0.9742682209735727, -0.001974719065561243])
 Panel2D{Float64}([0.9742682209735727, -0.001974719065561243], [0.9599897218294121, -0.0030518536407216934])
 Panel2D{Float64}([0.9599897218294121, -0.0030518536407216934], [0.9427280128266049, -0.004340093269408703])
 Panel2D{Float64}([0.9427280128266049, -0.004340093269408703], [0.9225950427718974, -0.005821397857625334])
 Panel2D{Float64}([0.9225950427718974, -0.005821397857625334], [0.8997213817017505, -0.0074835618802499525])
 Panel2D{Float64}([0.8997213817017505, -0.0074835618802499525], [0.8742553740855505, -0.00930660513155278])
 Panel2D{Float64}([0.8742553740855505, -0.00930660513155278], [0.8463621767547997, -0.011278505974986367])
 ⋮
 Panel2D{Float64}([0.8742553740855505, 0.024528991290121795], [0.8997213817017505, 0.01992712863731834])
 Panel2D{Float64}([0.8997213817017505, 0.01992712863731834], [0.9225950427718974, 0.015634911050710687])
 Panel2D{Float64}([0.9225950427718974, 0.015634911050710687], [0.9427280128266049, 0.01174051529843651])
 Panel2D{Float64}([0.9427280128266049, 0.01174051529843651], [0.9599897218294121, 0.008303069108260726])
 Panel2D{Float64}([0.9599897218294121, 0.008303069108260726], [0.9742682209735727, 0.005397295196262384])
 Panel2D{Float64}([0.9742682209735727, 0.005397295196262384], [0.985470908713026, 0.0030714609375638983])
 Panel2D{Float64}([0.985470908713026, 0.0030714609375638983], [0.9935251313189564, 0.0013774132702494887])
 Panel2D{Float64}([0.9935251313189564, 0.0013774132702494887], [0.9983786540671049, 0.00034600027495394753])
 Panel2D{Float64}([0.9983786540671049, 0.00034600027495394753], [1.0, -1.6616460838224803e-17])

You can compute the lift coefficient from the system.

cl = lift_coefficient(system)
0.7127871936239761

You can also compute the sectional lift, moment and pressure coefficients.

cls, cms, cps = surface_coefficients(system)
([0.0012701252075160581, 0.0019206412409908823, 0.0023203321585153563, 0.002634584429986258, 0.002901018757611492, 0.0031109619395088147, 0.0032792715891740393, 0.003398731138494381, 0.003490577088474002, 0.0035446744988400346  …  0.004071006912730074, 0.0025531183912331588, 0.0012263974049625287, 5.943071294792467e-5, -0.0008734792485322162, -0.0015927539699747367, -0.0020340330255393, -0.0021885688225591654, -0.0019770495579706157, -0.0012850621154423412], [-0.0012649769503887032, -0.0019035570488670395, -0.0022810497497611853, -0.002560539309580555, -0.002778164196321185, -0.002925646835817775, -0.003018069892994145, -0.0030504621582679655, -0.0030442220981581358, -0.0029927703844662636  …  -0.003437153094695431, -0.0022266402456688848, -0.0011007281018684829, -5.469691685783374e-5, 0.0008214474652249358, 0.00152530280658214, 0.0019768664308494424, 0.00215151712081944, 0.001959463611274343, 0.0012798533138566627], [0.3923231560582138, 0.2975910175377454, 0.24098267358203973, 0.20678079024304608, 0.18395185193928165, 0.16637933041425834, 0.15249603856882032, 0.1406135833794867, 0.13082924224853076, 0.12215717981664109  …  -0.14039643463913087, -0.09575579955394198, -0.0507741737797236, -0.0027653672646179572, 0.04674405527490677, 0.10104527417090836, 0.159722559243883, 0.22737428967272533, 0.30641978613201315, 0.3969927129606229])

Note the difference between the lift coefficient computed and the sum of the sectional lift coefficients; this is due to numerical errors in the solution procedure and modeling.

cl, sum(cls)
(0.7127871936239761, 0.7121901349530607)
Info

Support for drag prediction with boundary layer calculations will be added soon. For now, try out the amazing Webfoil developed by the MDOLab at University of Michigan – Ann Arbor!

Visualization

Let's see what the pressure and lift distribution curves look like over the airfoil. AeroFuse provides more helper functions for post-processing data. For example, you can make your fancy plots by segregating the values depending on the locations of the panels by defining the following function.

cp_upper, cp_lower = get_surface_values(panels, cps)
([(0.0008106729664475176, 0.1533261632456412), (0.004048107306969306, -0.8925692272375141), (0.010501979984008786, -1.3168612651425402), (0.020130435156700627, -1.427331062216862), (0.03287102859850757, -1.3621585058923351), (0.048641132671991494, -1.3173211755764096), (0.0673384722007489, -1.2429762772031356), (0.08884178776317603, -1.1958803576457742), (0.11301162210634938, -1.138138697501467), (0.13969122457982486, -1.095654237409902)  …  (0.8603087754201751, -0.14039643463913087), (0.8869883778936505, -0.09575579955394198), (0.911158212236824, -0.0507741737797236), (0.9326615277992512, -0.0027653672646179572), (0.9513588673280085, 0.04674405527490677), (0.9671289714014923, 0.10104527417090836), (0.9798695648432993, 0.159722559243883), (0.9894980200159912, 0.22737428967272533), (0.9959518926930306, 0.30641978613201315), (0.9991893270335525, 0.3969927129606229)], [(0.9959518926930306, 0.3923231560582138), (0.9894980200159912, 0.2975910175377454), (0.9798695648432993, 0.24098267358203973), (0.9671289714014923, 0.20678079024304608), (0.9513588673280085, 0.18395185193928165), (0.9326615277992512, 0.16637933041425834), (0.911158212236824, 0.15249603856882032), (0.8869883778936505, 0.1406135833794867), (0.8603087754201751, 0.13082924224853076), (0.8312924322762442, 0.12215717981664109)  …  (0.13969122457982486, 0.09672889133303653), (0.11301162210634938, 0.1237986002531195), (0.08884178776317603, 0.1557196293510733), (0.0673384722007489, 0.21058345890170038), (0.048641132671991494, 0.2829531459154514), (0.03287102859850757, 0.39260938760434805), (0.020130435156700627, 0.5461527496971887), (0.010501979984008786, 0.753638677029852), (0.004048107306969306, 0.966543918281474), (0.0008106729664475176, 0.914072404054249)])

Now let's plot the results!

# Pressure coefficients
plot(yflip = true, xlabel = L"(x/c)", ylabel = L"C_p", lw = 2)
plot!(cp_upper, label = L"$C_p$ Upper",
      ls = :dash, lw = 2, c = :cornflowerblue)
plot!(cp_lower, label = L"$C_p$ Lower",
      ls = :dash, lw = 2, c = :orange)
plot!(x_upper, -y_upper, label = "$(airfoil.name) Upper",
      ls = :solid, lw = 2, c = :cornflowerblue)
plot!(x_lower, -y_lower, label = "$(airfoil.name) Lower",
      ls = :solid, lw = 2, c = :orange)
# Lift coefficients
cl_upper, cl_lower = get_surface_values(panels, cls)

cl_plot = plot(xlabel = L"(x/c)", ylabel = L"C_l", lw = 2)
plot!(cl_upper, label = L"$C_l$ Upper",
      ls = :dash, lw = 2, c = :cornflowerblue)
plot!(cl_lower, label = L"$C_l$ Lower",
      ls = :dash, lw = 2, c = :orange)
plot!(x_lower, y_lower, label = "$(airfoil.name) Lower",
      ls = :solid, lw = 2, c = :orange)
plot!(x_upper, y_upper, label = "$(airfoil.name) Upper",
      ls = :solid, lw = 2, c = :cornflowerblue)
# Moment coefficients
cm_upper, cm_lower = get_surface_values(panels, cms)

cm_plot = plot(xlabel = L"(x/c)", ylabel = L"C_m")
plot!(cm_upper, label = L"$C_m$ Upper",
      ls = :dash, lw = 2, c = :cornflowerblue)
plot!(cm_lower, label = L"$C_m$ Lower",
      ls = :dash, lw = 2, c = :orange)
plot!(x_upper, y_upper, label = L"$C_m$ Upper",
      ls = :solid, lw = 2, c = :cornflowerblue)
plot!(x_lower, y_lower, label = L"$C_m$ Lower",
      ls = :solid, lw = 2, c = :orange)

Great! We've created our first airfoil and run an aerodynamic analysis in 2 dimensions. Now we can move on to analyzing an aircraft configuration in the Aircraft Aerodynamic Analysis tutorial.

This page was generated using Literate.jl.