Example Usage

This guide demonstrate how to use the basic capabilities of this package. Further information about using this package may be found in the Public API.

Loading Airfoil Geometry

Before any analysis can be performed, we need to load our airfoil geometry into XFOIL. This can be accomplished using the set_coordinates function.

using Xfoil, Plots, Printf
pyplot()

# read airfoil coordinates from a file
x, y = open("naca2412.dat", "r") do f
    x = Float64[]
    y = Float64[]
    for line in eachline(f)
        entries = split(chomp(line))
        push!(x, parse(Float64, entries[1]))
        push!(y, parse(Float64, entries[2]))
    end
    x, y
end

# load airfoil coordinates into XFOIL
Xfoil.set_coordinates(x,y)

# plot the airfoil geometry
scatter(x, y, label="", framestyle=:none, aspect_ratio=1.0, show=true)

Refining the Airfoil Geometry

It is often a good idea to refine airfoil discretizations prior to performing analyses using XFOIL. This may be accomplished using the pane function.

# repanel using XFOIL's `PANE` command
xr, yr = Xfoil.pane()

# plot the refined airfoil geometry
scatter(xr, yr, label="", framestyle=:none, aspect_ratio=1.0, show=true)

Defining Operating Conditions

We now need to define some operating conditions for our analysis. For inviscid analyses, only the angle of attack is required. For viscous analyses, the Reynolds number must also be specified.

# set operating conditions
alpha = -9:1:14 # range of angle of attacks, in degrees
re = 1e5 # Reynolds number

The solve_alpha function may now be used to perform an analysis to obtain the airfoil coefficients $c_l$, $c_d$, $c_{d_p}$, and $c_m$. Note that $c_{d_p}$ is profile drag. Skin friction drag may be obtained by subtracting the profile drag coefficient from the total drag coefficient i.e., $c_{d_f} = c_d - c_{d_p}$.

# initialize outputs
n_a = length(alpha)
c_l = zeros(n_a)
c_d = zeros(n_a)
c_dp = zeros(n_a)
c_m = zeros(n_a)
converged = zeros(Bool, n_a)

# determine airfoil coefficients across a range of angle of attacks
for i = 1:n_a
    c_l[i], c_d[i], c_dp[i], c_m[i], converged[i] = Xfoil.solve_alpha(alpha[i], re; iter=100, reinit=true)
end

# print results
println("Angle\t\tCl\t\tCd\t\tCm\t\tConverged")
for i = 1:n_a
    @printf("%8f\t%8f\t%8f\t%8f\t%d\n",alpha[i],c_l[i],c_d[i],c_m[i],converged[i])
end

# plot results
plot(alpha, c_l, label="", xlabel="Angle of Attack (degrees)", ylabel="Lift Coefficient", show=true)
plot(alpha, c_d, label="", xlabel="Angle of Attack (degrees)", ylabel="Drag Coefficient",
    overwrite_figure=false, show=true)
plot(alpha, c_m, label="", xlabel="Angle of Attack (degrees)", ylabel="Moment Coefficient",
    overwrite_figure=false, show=true)

Angle		Cl		Cd		Cm		Converged
-9.000000	-0.650846	0.064745	-0.049415	1
-8.000000	-0.662429	0.046357	-0.043245	1
-7.000000	-0.635637	0.034118	-0.035047	1
-6.000000	-0.565933	0.026936	-0.029727	1
-5.000000	-0.479794	0.022554	-0.025883	1
-4.000000	-0.390118	0.019767	-0.023503	1
-3.000000	-0.301726	0.018318	-0.021741	1
-2.000000	-0.150388	0.017026	-0.030251	1
-1.000000	0.071766	0.017390	-0.049186	1
0.000000	0.263166	0.016828	-0.065553	1
1.000000	0.404422	0.015968	-0.069629	1
2.000000	0.512138	0.015591	-0.066432	1
3.000000	0.611810	0.015749	-0.061841	1
4.000000	0.708466	0.016397	-0.057255	1
5.000000	0.804453	0.017169	-0.052536	1
6.000000	0.893460	0.018133	-0.046931	1
7.000000	0.969537	0.019458	-0.039801	1
8.000000	1.020497	0.023512	-0.030233	1
9.000000	1.087529	0.029120	-0.023964	1
10.000000	1.173720	0.035266	-0.020845	1
11.000000	1.242721	0.043276	-0.016251	1
12.000000	1.256900	0.054220	-0.007123	1
13.000000	1.197794	0.066990	0.003651	1
14.000000	0.849596	0.148350	-0.038913	1

Note that the order in which viscous analyses are performed matters since XFOIL uses boundary layer parameters corresponding to the last previously converged solution as an initial guess when solving for the current boundary layer parameters. This behavior can be disabled by passing the keyword argument pair reinit=true to solve_alpha

Sensitivity Analysis

Suppose we want to find the derivative of $c_l$, $c_d$, and $c_m$ with respect to the angle of attack. One approach to calculate these derivatives would be to use the finite difference method.

using Xfoil, Printf

# read airfoil into XFOIL
open("naca2412.dat", "r") do f
    x = Float64[]
    y = Float64[]
    for line in eachline(f)
        entries = split(chomp(line))
        push!(x, parse(Float64, entries[1]))
        push!(y, parse(Float64, entries[2]))
    end
    Xfoil.set_coordinates(x,y)
end

# repanel using XFOIL's `PANE` command
Xfoil.pane()

# set operating conditions
alpha = -9:1:14
re = 1e5
mach = 0.0

# set step size
h = 1e-6

# initialize outputs
n_a = length(alpha)
c_l_a = zeros(n_a)
c_d_a = zeros(n_a)
c_dp_a = zeros(n_a)
c_m_a = zeros(n_a)
converged = zeros(Bool, n_a)

for i = 1:n_a
    c_l1, c_d1, c_dp1, c_m1, converged[i] = Xfoil.solve_alpha(alpha[i], re; mach, iter=100, reinit=true)
    c_l2, c_d2, c_dp2, c_m2, converged[i] = Xfoil.solve_alpha(alpha[i]+h, re; mach, iter=100, reinit=true)
    c_l_a[i] = (c_l2 - c_l1)/h * 180/pi
    c_d_a[i] = (c_d2 - c_d1)/h * 180/pi
    c_m_a[i] = (c_m2 - c_m1)/h * 180/pi
end

# print results
println("Angle\t\tdClda\t\tdCdda\t\tdCmda\t\tConverged")
for i = 1:n_a
  @printf("%8f\t%8f\t%8f\t%8f\t%d\n",alpha[i],c_l_a[i],c_d_a[i],c_m_a[i],converged[i])
end

Angle		dClda		dCdda		dCmda		Converged
-9.000000	4.368002	-1.420669	0.050753	1
-8.000000	-0.595924	-0.987825	0.492841	1
-7.000000	3.092874	-0.537248	0.414648	1
-6.000000	4.582598	-0.250074	0.264979	1
-5.000000	5.114411	-0.079249	0.248431	1
-4.000000	5.137603	-0.072158	0.141209	1
-3.000000	5.279715	-0.026029	0.058558	1
-2.000000	8.234179	-0.042702	-0.136227	1
-1.000000	10.039118	-0.015135	-0.885506	1
0.000000	12.644209	-0.082045	-1.124300	1
1.000000	9.212672	-0.093061	-0.324636	1
2.000000	3.641180	0.067090	0.521405	1
3.000000	6.896780	-0.016704	0.126283	1
4.000000	5.982546	0.027492	0.212131	1
5.000000	6.180839	0.018997	0.228296	1
6.000000	5.383315	0.055856	0.295899	1
7.000000	3.761211	0.103208	0.470835	1
8.000000	2.940949	0.295101	0.506869	1
9.000000	4.686078	0.293471	0.276193	1
10.000000	5.637746	0.418566	0.028201	1
11.000000	3.889931	0.448186	0.215136	1
12.000000	-1.889114	0.640222	0.760204	1
13.000000	-5.154260	0.958681	0.299039	1
14.000000	1.558546	1.051198	-0.058802	1

A better approach might be to use the complex step method. To use this approach, however, we must use the complex-step enabled version of XFOIL provided by this package.

The complex-step version of each function is denoted by appending _cs to each function name. Note that there is no interaction between the two versions of XFOIL wrapped by this package, so if you wish to use the complex step version of the code you must append _cs to all function names.

For the complex step method to work, we also need to ensure that the imaginary portion of variables from previous iterations does not affect the solution for the current iteration. This may be achieved by setting reinit=true when calling solve_alpha.

using Xfoil, Printf

# read airfoil into XFOIL
open("naca2412.dat", "r") do f
    x = Float64[]
    y = Float64[]
    for line in eachline(f)
        entries = split(chomp(line))
        push!(x, parse(Float64, entries[1]))
        push!(y, parse(Float64, entries[2]))
    end
    Xfoil.set_coordinates_cs(x,y)
end

# repanel using XFOIL's `PANE` command
Xfoil.pane_cs()

# set operating conditions
alpha = -9:1:14
re = 1e5
mach = 0.0

# set step size
h = 1e-12im

# initialize outputs
n_a = length(alpha)
c_l_a = zeros(n_a)
c_d_a = zeros(n_a)
c_dp_a = zeros(n_a)
c_m_a = zeros(n_a)
converged = zeros(Bool, n_a)

for i = 1:n_a
    c_l, c_d, c_dp, c_m, converged[i] = Xfoil.solve_alpha_cs(alpha[i]+h, re; mach, iter=100, reinit=true)
    c_l_a[i] = imag(c_l)/imag(h) * 180/pi
    c_d_a[i] = imag(c_d)/imag(h) * 180/pi
    c_m_a[i] = imag(c_m)/imag(h) * 180/pi
end

# print results
println("Angle\t\tdClda\t\tdCdda\t\tdCmda\t\tConverged")
for i = 1:n_a
  @printf("%8f\t%8f\t%8f\t%8f\t%d\n",alpha[i],c_l_a[i],c_d_a[i],c_m_a[i],converged[i])
end

Angle		dClda		dCdda		dCmda		Converged
-9.000000	4.368078	-1.420675	0.050750	1
-8.000000	-0.595926	-0.987825	0.492841	1
-7.000000	3.094020	-0.537099	0.414477	1
-6.000000	4.583204	-0.250262	0.264980	1
-5.000000	5.114413	-0.079250	0.248431	1
-4.000000	5.137602	-0.072158	0.141209	1
-3.000000	5.279711	-0.026030	0.058559	1
-2.000000	8.234162	-0.042703	-0.136221	1
-1.000000	10.039099	-0.015135	-0.885503	1
0.000000	12.644207	-0.082045	-1.124300	1
1.000000	9.212668	-0.093061	-0.324636	1
2.000000	3.641152	0.067091	0.521409	1
3.000000	6.896776	-0.016704	0.126284	1
4.000000	5.982540	0.027492	0.212131	1
5.000000	6.180833	0.018997	0.228296	1
6.000000	5.383310	0.055856	0.295900	1
7.000000	3.761214	0.103208	0.470834	1
8.000000	2.940954	0.295101	0.506868	1
9.000000	4.686077	0.293470	0.276194	1
10.000000	5.637742	0.418565	0.028201	1
11.000000	3.889928	0.448186	0.215137	1
12.000000	-1.889107	0.640223	0.760202	1
13.000000	-5.154275	0.958683	0.299040	1
14.000000	1.555321	1.049194	-0.058519	1

Note that since XFOIL was not originally designed for sensitivity analysis, there is a distinct possibility that sensitivities may be non-physical. We therefore always recommend checking that computed sensitivities are realistic. Sometimes adjusting the number of panels or the step size will fix the computed sensitivies, but if all else fails a coarse step size (e.g. 1.0 degree) may be used with central finite differencing to artificially smooth the computed sensitivities.

Automated Angle of Attack Sweep

For performing angle of attack sweeps, the function alpha_sweep may also be used.

using Xfoil, Printf

# extract geometry
x = Float64[]
y = Float64[]

f = open("naca2412.dat", "r")

for line in eachline(f)
    entries = split(chomp(line))
    push!(x, parse(Float64, entries[1]))
    push!(y, parse(Float64, entries[2]))
end

close(f)

# set operating conditions
alpha = -9:1:14
re = 1e5

c_l, c_d, c_dp, c_m, converged = Xfoil.alpha_sweep(x, y, alpha, re, iter=100, zeroinit=false, printdata=true, reinit=true)


Angle		Cl		Cd		Cm		Converged
-9.000000	-0.650846	0.064745	-0.049415	1
-8.000000	-0.662429	0.046357	-0.043245	1
-7.000000	-0.635637	0.034118	-0.035047	1
-6.000000	-0.565933	0.026936	-0.029727	1
-5.000000	-0.479794	0.022554	-0.025883	1
-4.000000	-0.390118	0.019767	-0.023503	1
-3.000000	-0.301726	0.018318	-0.021741	1
-2.000000	-0.150388	0.017026	-0.030251	1
-1.000000	0.071766	0.017390	-0.049186	1
0.000000	0.263166	0.016828	-0.065553	1
1.000000	0.404422	0.015968	-0.069629	1
2.000000	0.512138	0.015591	-0.066432	1
3.000000	0.611810	0.015749	-0.061841	1
4.000000	0.708466	0.016397	-0.057255	1
5.000000	0.804453	0.017169	-0.052536	1
6.000000	0.893460	0.018133	-0.046931	1
7.000000	0.969537	0.019458	-0.039801	1
8.000000	1.020497	0.023512	-0.030233	1
9.000000	1.087529	0.029120	-0.023964	1
10.000000	1.173720	0.035266	-0.020845	1
11.000000	1.242721	0.043276	-0.016251	1
12.000000	1.256900	0.054220	-0.007123	1
13.000000	1.197794	0.066990	0.003651	1
14.000000	0.849596	0.148350	-0.038913	1

A version of alpha_sweep has also been implemented for use with the complex step version of XFOIL.

using Xfoil, Printf

# extract geometry
x = Float64[]
y = Float64[]

f = open("naca2412.dat", "r")

for line in eachline(f)
    entries = split(chomp(line))
    push!(x, parse(Float64, entries[1]))
    push!(y, parse(Float64, entries[2]))
end

close(f)

# set operating conditions
alpha = -9:1:14
re = 1e5
mach = 0.0

# set step size
h = 1e-20im

c_l, c_d, c_dp, c_m, converged = Xfoil.alpha_sweep_cs(x, y, alpha .+ h,
    re, mach=mach, iter=100, zeroinit=false, printdata=true, reinit=true)

println("Angle\t\tdClda\t\tdCdda\t\tdCmda\t\tConverged")
for i = 1:length(alpha)
    @printf("%8f\t%8f\t%8f\t%8f\t%d\n", alpha[i], imag(c_l[i])/imag(h)*180/pi, imag(c_d[i])/imag(h)*180/pi, imag(c_m[i])/imag(h)*180/pi, converged[i])
end