Examples

These examples show how to use VortexLattice for various geometries, flow conditions, and analyses. Many of these examples also provide a verification for the implementation of the vortex lattice method in this package.

• Examples
• • Steady State Analysis of a Wing
  • Steady State Analysis of a Wing with Dihedral
  • Steady State Analysis of a Wing and Tail
  • Sudden Acceleration of a Rectangular Wing into a Constant-Speed Forward Flight
  • Heaving Oscillations of a Rectangular Wing
  • OpenVSP Geometry Import

Steady State Analysis of a Wing

This example shows how to calculate aerodynamic coefficients and stability derivatives for a symmetric planar wing.

using VortexLattice

# geometry (right half of the wing)
xle = [0.0, 0.4]
yle = [0.0, 7.5]
zle = [0.0, 0.0]
chord = [2.2, 1.8]
theta = [2.0*pi/180, 2.0*pi/180]
phi = [0.0, 0.0]
fc = fill((xc) -> 0, 2) # camberline function for each section

# discretization parameters
ns = 12
nc = 6
spacing_s = Uniform()
spacing_c = Uniform()

# reference parameters
Sref = 30.0
cref = 2.0
bref = 15.0
rref = [0.50, 0.0, 0.0]
Vinf = 1.0
ref = Reference(Sref, cref, bref, rref, Vinf)

# freestream parameters
alpha = 1.0*pi/180
beta = 0.0
Omega = [0.0; 0.0; 0.0]
fs = Freestream(Vinf, alpha, beta, Omega)

# construct surface
grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
    fc = fc, spacing_s=spacing_s, spacing_c=spacing_c)

# create vector containing all surfaces
surfaces = [surface]

# we can use symmetry since the geometry and flow conditions are symmetric about the X-Z axis
symmetric = true

# perform steady state analysis
system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

# retrieve near-field forces
CF, CM = body_forces(system; frame=Wind())

# perform far-field analysis
CDiff = far_field_drag(system)

CD, CY, CL = CF
Cl, Cm, Cn = CM

The aerodynamic coefficients predicted by VortexLattice are nearly identical to those predicted by AVL.

We can also generate files to visualize the results in Paraview using the function write_vtk.

properties = get_surface_properties(system)

write_vtk("symmetric-planar-wing", surfaces, properties; symmetric)

For asymmetric flow conditions and/or to obtain accurate asymmetric stability derivatives we can use the keyword argument mirror when constructing the geometry to reflect the geometry across the X-Z plane prior to the analysis. We also set the symmetric flag to false since we are no longer using symmetry in the analysis.

# construct geometry with mirror image
grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
    fc=fc, spacing_s=spacing_s, spacing_c=spacing_c, mirror=true)

# symmetry is not used in the analysis
symmetric = false

# create vector containing all surfaces
surfaces = [surface]

# perform steady state analysis
system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

# retrieve near-field forces
CF, CM = body_forces(system; frame=Wind())

# perform far-field analysis
CDiff = far_field_drag(system)

CD, CY, CL = CF
Cl, Cm, Cn = CM

Once again, the aerodynamic coefficients predicted by VortexLattice are nearly identical to those predicted by AVL.

The stability derivatives are also very close to those predicted by AVL.

dCF, dCM = stability_derivatives(system)

CDa, CYa, CLa = dCF.alpha
Cla, Cma, Cna = dCM.alpha
CDb, CYb, CLb = dCF.beta
Clb, Cmb, Cnb = dCM.beta
CDp, CYp, CLp = dCF.p
Clp, Cmp, Cnp = dCM.p
CDq, CYq, CLq = dCF.q
Clq, Cmq, Cnq = dCM.q
CDr, CYr, CLr = dCF.r
Clr, Cmr, Cnr = dCM.r

Visualizing the geometry now shows the circulation distribution across the entire wing.

properties = get_surface_properties(system)

write_vtk("mirrored-planar-wing", surfaces, properties; symmetric)

Steady State Analysis of a Wing with Dihedral

This example shows how to calculate aerodynamic coefficients and stability derivatives for a simple wing with dihedral.

using VortexLattice

xle = [0.0, 0.4]
yle = [0.0, 7.5]
zle = [0.0, 3.0]
chord = [2.2, 1.8]
theta = [2.0*pi/180, 2.0*pi/180]
phi = [0.0, 0.0]
fc = fill((xc) -> 0, 2) #camberline function for each section

ns = 12
nc = 6
spacing_s = Uniform()
spacing_c = Uniform()
mirror = false
symmetric = true

Sref = 30.0
cref = 2.0
bref = 15.0
rref = [0.50, 0.0, 0.0]
Vinf = 1.0
ref = Reference(Sref, cref, bref, rref, Vinf)

alpha = 1.0*pi/180
beta = 0.0
Omega = [0.0; 0.0; 0.0]
fs = Freestream(Vinf, alpha, beta, Omega)

# declare symmetry
symmetric = true

# construct surface
grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
    fc = fc, spacing_s=spacing_s, spacing_c=spacing_c)

# create vector containing all surfaces
surfaces = [surface]

# perform steady state analysis
system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

# retrieve near-field forces
CF, CM = body_forces(system; frame=Wind())

# perform far-field analysis
CDiff = far_field_drag(system)

CD, CY, CL = CF
Cl, Cm, Cn = CM

The results predicted by VortexLattice are close to those predicted by AVL, with the difference primarily explained by the manner in which the normal vector is defined in VortexLattice and AVL, respectively.

If we set the normal vectors in VortexLattice equal to those used in AVL, the results are even closer, though not necessarily more accurate.

using LinearAlgebra

# function to construct a normal vector the way AVL does
#  - `ds` is a line representing the leading edge
#  - `theta` is the incidence angle, taken as a rotation (+ by RH rule) about
#        the surface's spanwise axis projected onto the Y-Z plane.
function avl_normal_vector(ds, theta)

    st, ct = sincos(theta)

    # bound vortex vector
    bhat = ds/norm(ds)

    # chordwise strip normal vector
    shat = [0, -ds[3], ds[2]]/sqrt(ds[2]^2+ds[3]^2)

    # camberline vector
    chat = [ct, -st*shat[2], -st*shat[3]]

    # normal vector perpindicular to camberline and bound vortex for entire chordwise strip
    ncp = cross(chat, ds)
    return ncp / norm(ncp) # normal vector used by AVL
end

# new normal vector
ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180)

# overwrite normal vector for each panel
for i = 1:length(surface)
    surface[i] = set_normal(surface[i], ncp)
end

# create vector containing all surfaces
surfaces = [surface]

# perform steady state analysis
system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

# retrieve near-field forces
CF, CM = body_forces(system; frame=Wind())

# perform far-field analysis
CDiff = far_field_drag(system)

CD, CY, CL = CF
Cl, Cm, Cn = CM

properties = get_surface_properties(system)

write_vtk("wing-with-dihedral", surfaces, properties; symmetric)

Steady State Analysis of a Wing and Tail

This example shows how to calculate aerodynamic coefficients and stability derivatives for multiple lifting surfaces.

using VortexLattice

# wing
xle = [0.0, 0.2]
yle = [0.0, 5.0]
zle = [0.0, 1.0]
chord = [1.0, 0.6]
theta = [2.0*pi/180, 2.0*pi/180]
phi = [0.0, 0.0]
fc = fill((xc) -> 0, 2) # camberline function for each section
ns = 12
nc = 6
spacing_s = Uniform()
spacing_c = Uniform()
mirror = false

# horizontal stabilizer
xle_h = [0.0, 0.14]
yle_h = [0.0, 1.25]
zle_h = [0.0, 0.0]
chord_h = [0.7, 0.42]
theta_h = [0.0, 0.0]
phi_h = [0.0, 0.0]
fc_h = fill((xc) -> 0, 2) #camberline function for each section
ns_h = 6
nc_h = 3
spacing_s_h = Uniform()
spacing_c_h = Uniform()
mirror_h = false

# vertical stabilizer
xle_v = [0.0, 0.14]
yle_v = [0.0, 0.0]
zle_v = [0.0, 1.0]
chord_v = [0.7, 0.42]
theta_v = [0.0, 0.0]
phi_v = [0.0, 0.0]
fc_v = fill((xc) -> 0, 2) #camberline function for each section
ns_v = 5
nc_v = 3
spacing_s_v = Uniform()
spacing_c_v = Uniform()
mirror_v = false

Sref = 9.0
cref = 0.9
bref = 10.0
rref = [0.5, 0.0, 0.0]
Vinf = 1.0
ref = Reference(Sref, cref, bref, rref, Vinf)

alpha = 5.0*pi/180
beta = 0.0
Omega = [0.0; 0.0; 0.0]
fs = Freestream(Vinf, alpha, beta, Omega)

symmetric = [true, true, false]

# generate surface panels for wing
wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
    mirror=mirror, fc = fc, spacing_s=spacing_s, spacing_c=spacing_c)

# generate surface panels for horizontal tail
hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h;
    mirror=mirror_h, fc=fc_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h)
translate!(hgrid, [4.0, 0.0, 0.0])
translate!(htail, [4.0, 0.0, 0.0])

# generate surface panels for vertical tail
vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v;
    mirror=mirror_v, fc=fc_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v)
translate!(vgrid, [4.0, 0.0, 0.0])
translate!(vtail, [4.0, 0.0, 0.0])

grids = [wgrid, hgrid, vgrid]
surfaces = [wing, htail, vtail]
surface_id = [1, 2, 3]

system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, surface_id=surface_id)

CF, CM = body_forces(system; frame=Wind())

CDiff = far_field_drag(system)

CD, CY, CL = CF
Cl, Cm, Cn = CM

The results predicted by VortexLattice are close to those predicted by AVL (with the finite core model disabled in AVL), with the difference primarily explained by the manner in which the normal vector is defined in VortexLattice and AVL, respectively.

If we set the normal vectors in VortexLattice equal to those used in AVL, the results are closer, though not necessarily more accurate.

using LinearAlgebra

# function to construct a normal vector the way AVL does
#  - `ds` is a line representing the leading edge
#  - `theta` is the incidence angle, taken as a rotation (+ by RH rule) about
#        the surface's spanwise axis projected onto the Y-Z plane.
function avl_normal_vector(ds, theta)

    st, ct = sincos(theta)

    # bound vortex vector
    bhat = ds/norm(ds)

    # chordwise strip normal vector
    shat = [0, -ds[3], ds[2]]/sqrt(ds[2]^2+ds[3]^2)

    # camberline vector
    chat = [ct, -st*shat[2], -st*shat[3]]

    # normal vector perpindicular to camberline and bound vortex for entire chordwise strip
    ncp = cross(chat, ds)
    return ncp / norm(ncp) # normal vector used by AVL
end

# new normal vector for the wing
ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180)

# overwrite normal vector for each wing panel
for i = 1:length(wing)
    wing[i] = set_normal(wing[i], ncp)
end
surfaces[1] = wing

# perform steady state analysis
system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

# retrieve near-field forces
CF, CM = body_forces(system; frame=Wind())

# perform far-field analysis
CDiff = far_field_drag(system)

CD, CY, CL = CF
Cl, Cm, Cn = CM

To achieve a theoretically identical setup as AVL we can place all our panels in the X-Y plane and then set the normal vector manually to match the actual lifting geometry. In our case this involves removing the small amount of twist on the wing when creating the wing surface panels.

using VortexLattice

# wing
xle = [0.0, 0.2]
yle = [0.0, 5.0]
zle = [0.0, 1.0]
chord = [1.0, 0.6]
theta = [0.0, 0.0]
phi = [0.0, 0.0]
fc = fill((xc) -> 0, 2) # camberline function for each section
ns = 12
nc = 6
spacing_s = Uniform()
spacing_c = Uniform()
mirror = false

# horizontal stabilizer
xle_h = [0.0, 0.14]
yle_h = [0.0, 1.25]
zle_h = [0.0, 0.0]
chord_h = [0.7, 0.42]
theta_h = [0.0, 0.0]
phi_h = [0.0, 0.0]
fc_h = fill((xc) -> 0, 2) # camberline function for each section
ns_h = 6
nc_h = 3
spacing_s_h = Uniform()
spacing_c_h = Uniform()
mirror_h = false

# vertical stabilizer
xle_v = [0.0, 0.14]
yle_v = [0.0, 0.0]
zle_v = [0.0, 1.0]
chord_v = [0.7, 0.42]
theta_v = [0.0, 0.0]
phi_v = [0.0, 0.0]
fc_v = fill((xc) -> 0, 2) # camberline function for each section
ns_v = 5
nc_v = 3
spacing_s_v = Uniform()
spacing_c_v = Uniform()
mirror_v = false

Sref = 9.0
cref = 0.9
bref = 10.0
rref = [0.5, 0.0, 0.0]
Vinf = 1.0
ref = Reference(Sref, cref, bref, rref, Vinf)

alpha = 5.0*pi/180
beta = 0.0
Omega = [0.0; 0.0; 0.0]
fs = Freestream(Vinf, alpha, beta, Omega)

symmetric = [true, true, false]

# generate surface panels for wing
wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
    mirror=mirror, fc=fc, spacing_s=spacing_s, spacing_c=spacing_c)

# generate surface panels for horizontal tail
hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h;
    mirror=mirror_h, fc=fc_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h)
translate!(hgrid, [4.0, 0.0, 0.0])
translate!(htail, [4.0, 0.0, 0.0])

# generate surface panels for vertical tail
vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v;
    mirror=mirror_v, fc=fc_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v)
translate!(vgrid, [4.0, 0.0, 0.0])
translate!(vtail, [4.0, 0.0, 0.0])

# now set normal vectors manually
ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180)

# overwrite normal vector for each wing panel
for i = 1:length(wing)
    wing[i] = set_normal(wing[i], ncp)
end

grids = [wgrid, hgrid, vgrid]
surfaces = [wing, htail, vtail]
surface_id = [1, 2, 3]

system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, surface_id=surface_id)

CF, CM = body_forces(system; frame=Stability())

CDiff = far_field_drag(system)

CD, CY, CL = CF
Cl, Cm, Cn = CM

The resulting aerodynamic coefficients now match very closely with AVL.

By comparing these results with previous results we can see exactly how much restricting surface panels in the X-Y plane changes the results from the vortex lattice method.

properties = get_surface_properties(system)

write_vtk("wing-tail", surfaces, properties; symmetric)

Sudden Acceleration of a Rectangular Wing into a Constant-Speed Forward Flight

This example shows how to predict the transient forces and moments on a rectangular wing when suddenly accelerated into forward flight at a five degree angle.

# Katz and Plotkin: Figures 13.34 and 13.35
# AR = [4, 8, 12, 20, ∞]
# Vinf*Δt/c = 1/16
# α = 5°

using VortexLattice

AR = [4, 8, 12, 20, 1e3] # last aspect ratio is essentially infinite

system = Vector{Any}(undef, length(AR))
surface_history = Vector{Any}(undef, length(AR))
property_history = Vector{Any}(undef, length(AR))
wake_history = Vector{Any}(undef, length(AR))
CF = Vector{Vector{Vector{Float64}}}(undef, length(AR))
CM = Vector{Vector{Vector{Float64}}}(undef, length(AR))

# non-dimensional time (t*Vinf/c)
t = range(0.0, 10.0, step=1/16)

# chord length
c = 1

# time step
dt = [t[i+1]-t[i] for i = 1:length(t)-1]

for i = 1:length(AR)

    # span length
    b = AR[i]*c

    # planform area
    S = b*c

    # geometry
    xle = [0.0, 0.0]
    yle = [-b/2, b/2]
    zle = [0.0, 0.0]
    chord = [c, c]
    theta = [0.0, 0.0]
    phi = [0.0, 0.0]
    fc = fill((xc) -> 0, 2) # camberline function for each section
    ns = 13
    nc = 4
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = false
    symmetric = false

    # reference parameters
    cref = c
    bref = b
    Sref = S
    rref = [0.0, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    # freestream parameters
    alpha = 5.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    # create vortex rings
    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, fc=fc, spacing_s=spacing_s, spacing_c=spacing_c)

    # create vector containing surfaces
    surfaces = [surface]

    # run analysis
    system[i], surface_history[i], property_history[i], wake_history[i] =
        unsteady_analysis(surfaces, ref, fs, dt; symmetric, wake_finite_core = false)

    # extract forces at each time step
    CF[i], CM[i] = body_forces_history(system[i], surface_history[i],
        property_history[i]; frame=Wind())

end

We can visualize the solution using the write_vtk function.

write_vtk("acceleration-AR4", surface_history[1], property_history[1],
    wake_history[1], dt; symmetric=false)

The transient lift and drag coefficients are similar to those shown in Figures 13.34 and 13.35 of Low-Speed Aerodynamics by Katz and Plotkin.

using Plots
pyplot()

# lift coefficient plot
plot(
    xlim = (0.0, 10.0),
    xticks = 0.0:1.0:10.0,
    xlabel = "\$ \\frac{U_\\infty t}{c} \$",
    ylim = (0.0, 0.55),
    yticks = 0.0:0.1:0.5,
    ylabel = "\$ C_{L} \$",
    grid = false,
    overwrite_figure=false
    )

for i = 1:length(AR)
    CL = [CF[i][j][3] for j = 1:length(CF[i])]
    plot!(t[2:end], CL, label="AR = $(AR[i])")
end

plot!(show=true)

# drag coefficient plot
plot(
    xlim = (0.0, 10.0),
    xticks = 0.0:1.0:10.0,
    xlabel = "\$ \\frac{U_\\infty t}{c} \$",
    ylim = (0.0, 0.030),
    yticks = 0.0:0.005:0.03,
    ylabel = "\$ C_{D} \$",
    grid = false,
    overwrite_figure=false
    )

for i = 1:length(AR)
    CD = [CF[i][j][1] for j = 1:length(CF[i])]
    plot!(t[2:end], CD, label="AR = $(AR[i])")
end

plot!(show=true)

We modeled the problem in the body-fixed reference frame (which for this problem is more straightforward), but we could have also modeled the problem in the global reference frame.

# Katz and Plotkin: Figures 13.34 and 13.35
# AR = [4, 8, 12, 20, ∞]
# Vinf*Δt/c = 1/16
# α = 5°

using VortexLattice

AR = [4, 8, 12, 20, 1e3] # last aspect ratio is essentially infinite

system_t = Vector{Any}(undef, length(AR))
surface_history_t = Vector{Any}(undef, length(AR))
property_history_t = Vector{Any}(undef, length(AR))
wake_history_t = Vector{Any}(undef, length(AR))
CF_t = Vector{Vector{Vector{Float64}}}(undef, length(AR))
CM_t = Vector{Vector{Vector{Float64}}}(undef, length(AR))

# non-dimensional time (t*Vinf/c)
t = range(0.0, 10.0, step=1/16)

# chord length
c = 1

# time step
dt = [t[i+1]-t[i] for i = 1:length(t)-1]

for i = 1:length(AR)

    # span length
    b = AR[i]*c

    # planform area
    S = b*c

    # geometry
    xle = [0.0, 0.0]
    yle = [-b/2, b/2]
    zle = [0.0, 0.0]
    chord = [c, c]
    theta = [0.0, 0.0]
    phi = [0.0, 0.0]
    fc = fill((xc) -> 0, 2) # camberline function for each section
    ns = 13
    nc = 4
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = false
    symmetric = false

    # reference parameters
    cref = c
    bref = b
    Sref = S
    rref = [0.0, 0.0, 0.0]
    Vinf = 1.0 # reference velocity is 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    # freestream parameters
    Vinf = 0.0 # freestream velocity is 0.0
    alpha = 5.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    # create vortex rings
    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, fc=fc, spacing_s=spacing_s, spacing_c=spacing_c)

    # create vector containing surfaces at each time step
    surfaces = [[VortexLattice.translate(surface,
        -t[it]*[cos(alpha), 0, sin(alpha)])] for it = 1:length(t)]

    # run analysis
    system_t[i], surface_history_t[i], property_history_t[i], wake_history_t[i] =
        unsteady_analysis(surfaces, ref, fs, dt; symmetric, wake_finite_core = false)

    # extract forces at each time step
    CF_t[i], CM_t[i] = body_forces_history(system_t[i], surface_history_t[i],
        property_history_t[i]; frame=Wind())

end

As can be seen, the transient lift and drag coefficients for the two setups are identical.

using Plots
pyplot()

# lift coefficient plot
plot(
    xlim = (0.0, 10.0),
    xticks = 0.0:1.0:10.0,
    xlabel = "\$ \\frac{U_\\infty t}{c} \$",
    ylim = (0.0, 0.55),
    yticks = 0.0:0.1:0.5,
    ylabel = "\$ C_{L} \$",
    grid = false,
    overwrite_figure=false
    )

for i = 1:length(AR)
    CL = [CF_t[i][j][3] for j = 1:length(CF_t[i])]
    plot!(t[2:end], CL, label="AR = $(AR[i])")
end

plot!(show=true)

# drag coefficient plot
plot(
    xlim = (0.0, 10.0),
    xticks = 0.0:1.0:10.0,
    xlabel = "\$ \\frac{U_\\infty t}{c} \$",
    ylim = (0.0, 0.030),
    yticks = 0.0:0.005:0.03,
    ylabel = "\$ C_{D} \$",
    grid = false,
    overwrite_figure=false
    )

for i = 1:length(AR)
    CD = [CF_t[i][j][1] for j = 1:length(CF_t[i])]
    plot!(t[2:end], CD, label="AR = $(AR[i])")
end

plot!(show=true)

Visualizing the solution shows the movement of the body in the global reference frame.

write_vtk("acceleration-AR4-moving", surface_history_t[1], property_history_t[1],
    wake_history_t[1], dt; symmetric=false)

For infinite aspect ratios, the problem degenerates into the analysis of the sudden acceleration of a 2D flat plate, for which we have an analytical solution through the work of Herbert Wagner.

# See Katz and Plotkin: Figure 13.37
# AR = ∞
# Vinf*Δt/c = 1/16
# α = 5°

# essentially infinite aspect ratio
AR = 1e3

# chord length
c = 1

# span length
b = AR*c

# planform area
S = b*c

# geometry
xle = [0.0, 0.0]
yle = [-b/2, b/2]
zle = [0.0, 0.0]
chord = [c, c]
theta = [0.0, 0.0]*pi/180
phi = [0.0, 0.0]
fc = fill((xc) -> 0, 2) # camberline function for each section
ns = 1
nc = 4
spacing_s = Uniform()
spacing_c = Uniform()
mirror = false
symmetric = false

# reference parameters
cref = c
bref = b
Sref = S
rref = [0.0, 0.0, 0.0]
Vinf = 1.0
ref = Reference(Sref, cref, bref, rref, Vinf)

# freestream parameters
alpha = 5.0*pi/180
beta = 0.0
Omega = [0.0; 0.0; 0.0]
fs = Freestream(Vinf, alpha, beta, Omega)

# non-dimensional time (t*Vinf/c)
t = range(0.0, 7.0, step=1/8)

# time step
dt = [(t[i+1]-t[i]) for i = 1:length(t)-1]

# create vortex rings
grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
    mirror=mirror, fc=fc, spacing_s=spacing_s, spacing_c=spacing_c)

# create vector containing all surfaces
surfaces = [surface]

# run steady analysis
system = steady_analysis(surfaces, ref, fs; symmetric)

# extract steady forces
CFs, CMs = body_forces(system; frame=Wind())

# run transient analysis
system, surface_history, property_history, wake_history = unsteady_analysis(
    surfaces, ref, fs, dt; symmetric=symmetric)

# extract transient forces
CF, CM = body_forces_history(system, surface_history, property_history; frame=Wind())

The results from VortexLattice compare very well with the analytical solution provided by Wagner. As discussed in Low Speed Aerodynamics by Katz and Plotkin, the difference between the curves can be attributed to the finite acceleration rate during the first time step, which increases the lift sharply during the acceleration and then increases it moderately later.

# lift coefficient plot
plot(
    xlim = (0.0, 7.0),
    xticks = 0.0:1.0:7.0,
    xlabel = "\$ \\frac{U_\\infty t}{c} \$",
    ylim = (0.0, 1.0),
    yticks = 0.0:0.1:1.0,
    ylabel = "\$ C_{L} \$",
    grid = false,
    overwrite_figure=false
    )

# Computational Results
CL = getindex.(CF, 3)
CLs = getindex(CFs, 3)
plot!(t[2:end], CL./CLs, label="VortexLattice")

# Wagner's Function (using approximation of R. T. Jones)
Φ(t) = 1 - 0.165*exp(-0.045*t) - 0.335*exp(-0.3*t)

plot!(t, Φ.(2*t), label = "Wagner's Function")

plot!(show=true)

Heaving Oscillations of a Rectangular Wing

This example shows how to predict the transient forces and moments for a heaving rectangular wing.

# Katz and Plotkin: Figures 13.38a
# AR = 4
# k = ω*c/(2*Vinf) = [0.5, 0.3, 0.1]
# c = [1.0, 0.6, 0.2]
# α = -5°

using VortexLattice

# forward velocity
Vinf = 1

# angle of attack
alpha = -5*pi/180

# aspect ratio
AR = 4

# chord lengths
c = [1.0, 0.6, 0.2]

# reduced frequency
k = [0.5, 0.3, 0.1]

t = Vector{Vector{Float64}}(undef, length(k))
CF = Vector{Vector{Vector{Float64}}}(undef, length(k))
CM = Vector{Vector{Vector{Float64}}}(undef, length(k))

for i = 1:length(k)

    # span length
    b = AR*c[i]

    # geometry
    xle = [0.0, 0.0]
    yle = [0.0, b/2]
    zle = [0.0, 0.0]
    chord = [c[i], c[i]]
    theta = [0.0, 0.0]
    phi = [0.0, 0.0]
    fc = fill((xc) -> 0, 2) # camberline function for each section
    ns = 13
    nc = 4
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = false
    symmetric = true

    # reference parameters
    cref = c[i]
    bref = b
    Sref = b*c[i]
    rref = [0.0, 0.0, 0.0]
    ref = Reference(Sref, cref, bref, rref, Vinf)

    # angular frequency
    ω = 2*Vinf*k[i]/c[i]

    # time
    t[i] = range(0.0, 9*pi/ω, length = 100)
    dt = t[i][2:end] - t[i][1:end-1]
    dt = Vinf*dt

    # heaving amplitude
    h = 0.1*c[i]

    # use forward and vertical velocity at beginning of each time step
    Xdot = Vinf*cos(alpha)
    Zdot = Vinf*sin(alpha) .- h*cos.(ω*t[i][1:end-1])

    # freestream parameters for each time step
    fs = trajectory_to_freestream(dt; Xdot, Zdot)

    # surface panels
    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, fc=fc, spacing_s=spacing_s, spacing_c=spacing_c)

    # create vector containing all surfaces
    surfaces = [surface]

    # run analysis
    system, surface_history, property_history, wake_history = unsteady_analysis(
        surfaces, ref, fs, dt; symmetric=symmetric, nwake = 50)

    # extract forces at each time step (uses instantaneous velocity as reference)
    CF[i], CM[i] = body_forces_history(system, surface_history, property_history; frame=Wind())

end

Plotting the results reveals that the results are similar to the results in Figure 13.34 of Low-Speed Aerodynamic by Katz and Plotkin, which verifies the unsteady vortex lattice method implementation in VortexLattice.

using Plots
pyplot()

# lift coefficient plot
plot(
    xlim = (6*pi, 8*pi),
    xticks = ([6*pi, 13*pi/2, 7*pi, 15*pi/2, 8*pi], ["\$ 0 \$",
        "\$ \\frac{\\pi}{2} \$", "\$ \\pi \$", "\$ \\frac{3\\pi}{2} \$",
        "\$ 2\\pi \$"]),
    xlabel = "\$ ω \\cdot t \$",
    ylim = (-1.0, 0.1),
    yticks = -1.0:0.2:0.0,
    yflip = true,
    ylabel = "\$ C_{L} \$",
    grid = false,
    )

for i = 1:length(k)
    # extract ω
    ω = 2*Vinf*k[i]/c[i]

    # extract ω*t (use time at the beginning of the time step)
    ωt = ω*t[i][1:end-1]

    # extract CL
    CL = [CF[i][it][3] for it = 1:length(t[i])-1]

    plot!(ωt, CL, label="\$ k = \\frac{\\omega c}{2 U_\\infty} = $(k[i]) \$")
end

plot!(show=true)

Visualizing the k=0.5 case in ParaView yields the following animation.

OpenVSP Geometry Import

This example shows how to import a wing geometry created using OpenVSP into VortexLattice for analysis. We'll make use of the default swept wing inside OpenVSP with a few minor changes.

Start up OpenVSP and create the default wing. Change the airfoil to NACA 2412 sections so that our wing has a camber to it. VortexLattice will make use of the cambersurface computed by OpenVSP when simulating it. The number of spanwise panels was increased to 20 per semispan in this example.

Once the geometry has been created, write out a DegenGeom file by selecting DegenGeom in the Analysis tab in OpenVSP. We only require a DegenGeom file in the csv format. The example DegenGeom file named samplewing.csv provided in docs/src was created in this manner.

The DegenGeom file can be imported into VortexLattice by using the functions read_degengeom and import_vsp. The read_degengeom function reads the DegenGeom file into an array of components suitable for use inside Julia. The import_vsp function imports required components from the array as specified by the user.

In the following example code, a steady state analysis is performed on the sample wing imported from OpenVSP and results are visualized in Paraview.

using VortexLattice

Sref = 45.0
cref = 2.5
bref = 18.0
rref = [0.625, 0.0, 0.0]
Vinf = 1.0
ref = Reference(Sref, cref, bref, rref, Vinf)

alpha = 1.0*pi/180
beta = 0.0
Omega = [0.0; 0.0; 0.0]
fs = Freestream(Vinf, alpha, beta, Omega)

# Import components inside Degengeom file into Julia
comp = read_degengeom("samplewing.csv")

# Use the first (and only) imported component to create the lifting surface
grid, surface = import_vsp(comp[1]; mirror=true)

symmetric = false
surfaces = [surface]

system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)
properties = get_surface_properties(system)
write_vtk("samplewing", surfaces, properties; symmetric)