Aeroelastic Stability Analysis of a Highly Flexible Wing

In this example, we demonstrate how to perform a three-dimensional aeroelastic analysis of a highly flexible cantilever wing.

Tip

This example is also available as a Jupyter notebook: cantilever.ipynb.

The wing we are considering in this example was created by modifying Daedalus aircraft data and is therefore representative of a high-altitude long-endurance wing. It has a 16 meter span (from root to tip) and a 1 meter chord. To model the wing's aerodynamics, we use a lifting line model. To model the wing's structure, we use a geometrically exact beam theory model.

using Aeroelasticity, GXBeam, DifferentialEquations, LinearAlgebra

# --- Initial Setup --- #

# discretization
N = 8 # number of elements

# geometric properties
span = 16 # m
chord = 1 # m
xref = 0.5 # normalized reference location (from leading edge)
xcg = 0.5 # center of gravity (from leading edge)

# stiffness properties
GJ = 1e4 # N*m^2 (torsional rigidity)
EIyy = 2e4 # N*m^2 (flat bending rigidity)
EIzz = 4e6 # N*m^2 (chord bending rigidity)

# inertial properties
mu = 0.75 # kg/m (mass per unit length)
i11 = 0.1 # kg*m (moment of inertia about elastic axis)
i22 = 0.0375 # moment of inertia about beam y-axis
i33 = 0.0625 # moment of inertia about beam z-axis

# freestream properties
Vinf = 10.0 # m/s (velocity)
rho = 0.0889 # kg/m^3 (air density at 20 km)
alpha = 2*pi/180 # angle of attack
c = 343.0 # m/s (air speed of sound)
beta = sqrt(1 - Vinf^2/c^2) # Prandtl-Glauert compressibility correction factor

# aerodynamic section properties
a = xref - 0.5 # normalized reference location (relative to semi-chord)
b = chord / 2 # m (semi-chord)
a0 = 2*pi # lift slope (for each section)
alpha0 = 0 # zero lift angle of attack (for each section)
cd0 = 0
cm0 = 0

# define geometry (assume NED coordinate frame)
xpt = range(0, 0, length=N+1) # point x-coordinates (in body frame)
ypt = range(0, span, length=N+1) # point y-coordinates (in body frame)
zpt = range(0, 0, length=N+1) # point z-coordinates (in body frame)
points = [[xpt[i],ypt[i],zpt[i]] for i = 1:N+1]
start = 1:N # starting point of each beam element
stop = 2:N+1 # ending point of each beam element
e1 = [0, 1,  0] # beam x-axis
e2 = [1, 0,  0] # beam y-axis
e3 = [0, 0, -1] # beam z-axis
frames = fill([e1 e2 e3], N) # local to body frame transformation
compliance = fill(Diagonal([0, 0, 0, 1/GJ, 1/EIyy, 1/EIzz]), N) # compliance matrix
mass = fill(Diagonal([mu, mu, mu, i11, i22, i33]), N) # mass matrix
damping = fill(fill(1e-4, 6), N) # stiffness proportional structural damping coefficients
assembly = GXBeam.Assembly(points, start, stop; frames, compliance, mass, damping)

prescribed_conditions = Dict(
    # fixed left edge
    1 => GXBeam.PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
)

# define GXBeam system
system = DynamicSystem(assembly)

# --- Define Submodels --- #

# construct section models (we use Peters' finite state model in this case)
section_models = fill(Peters{6}(), N)

# construct aerodynamic model
aerodynamic_model = LiftingLine(section_models)

# construct structural model
structural_model = GXBeamAssembly(system; structural_damping=true)

# define submodels
submodels = (aerodynamic_model, structural_model)

# --- Define Initial Parameters --- #

V = [-Vinf*cos(alpha), 0.0, -Vinf*sin(alpha)] # m/s (freestream velocity)

# define parameters for each lifting line section
section_parameters = fill([a, b, a0, alpha0, cd0, cm0], N)

# define parameters for the lifting line model
liftingline_parameters = LiftingLineParameters(section_parameters)

# define parameters for the geometrically exact beam theory model
gxbeam_parameters = GXBeamParameters(assembly)

# define parameters for the coupling
coupling_parameters = LiftingLineGXBeamParameters(V, rho, beta;
    prescribed_conditions = prescribed_conditions,
    gravity = [-9.81*sin(alpha), 0, 9.81*cos(alpha)])

# combine parameters
parameters = (liftingline_parameters, gxbeam_parameters, coupling_parameters)

# --- Define Coupled Model --- #

model = CoupledModel(submodels, parameters; symbolic=false)

# --- Perform Analysis --- #

# loop through freestream velocities
Vinf = vcat(0.1, 0.25:0.25:50)

# eigenvalue/eigenvector storage
nev = 12*N
λ = zeros(ComplexF64, nev, length(Vinf))
Uλ = zeros(ComplexF64, nev, number_of_states(model), length(Vinf))
Vλ = zeros(ComplexF64, number_of_states(model), nev, length(Vinf))

# initial guess for state variables
x0 = zeros(number_of_states(model))

# loop through each velocity
for i = 1:length(Vinf)

    # --- Update Parameters --- #

    V = [-Vinf[i]*cos(alpha), 0.0, -Vinf[i]*sin(alpha)] # m/s (freestream velocity)

    # define parameters for each lifting line section
    section_parameters = fill([a, b, a0, alpha0, cd0, cm0], N)

    # define parameters for the lifting line model
    liftingline_parameters = LiftingLineParameters(section_parameters)

    # define parameters for the geometrically exact beam theory model
    gxbeam_parameters = GXBeamParameters(assembly)

    # define parameters for the coupling
    coupling_parameters = LiftingLineGXBeamParameters(V, rho, beta;
        prescribed_conditions = prescribed_conditions,
        gravity = [-9.81*sin(alpha), 0, 9.81*cos(alpha)])

    # combine parameters
    parameters = (liftingline_parameters, gxbeam_parameters, coupling_parameters)

    # --- Perform Analysis --- #

    # define an ODEFunction for this model
    f = ODEFunction(model, parameters)

    # find equilibrium point
    sol = solve(NonlinearProblem(SteadyStateProblem(f, x0, parameters)))

    # use state variables from steady state operating conditions
    x = sol.u

    # linearize about steady state operating conditions
    K, M = linearize(model, x, parameters)

    # perform linear stability analysis
    λi, Uλi, Vλi = sparse_eigen(K, M; nev=nev)

    # --- Correlate Eigenvalues --- #

    if i > 1
        # previous left eigenvector matrix
        Uλpi = Uλ[:,:,i-1]

        # use correlation matrix to correlate eigenmodes
        perm, corruption = Aeroelasticity.correlate_eigenmodes(Uλpi, M, Vλi)

        # re-arrange eigenmodes
        λi = λi[perm]
        Uλi = Uλi[perm,:]
        Vλi = Vλi[:,perm]
    end

    # save eigenvalues/eigenvectors
    λ[:,i] = λi
    Uλ[:,:,i] = Uλi
    Vλ[:,:,i] = Vλi

    # update initial guess for the state variables
    x0 .= x
end

┌ Warning: Assignment to `V` in soft scope is ambiguous because a global variable by the same name exists: `V` will be treated as a new local. Disambiguate by using `local V` to suppress this warning or `global V` to assign to the existing global variable.
└ @ cantilever-stability.md:144
┌ Warning: Assignment to `section_parameters` in soft scope is ambiguous because a global variable by the same name exists: `section_parameters` will be treated as a new local. Disambiguate by using `local section_parameters` to suppress this warning or `global section_parameters` to assign to the existing global variable.
└ @ cantilever-stability.md:147
┌ Warning: Assignment to `liftingline_parameters` in soft scope is ambiguous because a global variable by the same name exists: `liftingline_parameters` will be treated as a new local. Disambiguate by using `local liftingline_parameters` to suppress this warning or `global liftingline_parameters` to assign to the existing global variable.
└ @ cantilever-stability.md:150
┌ Warning: Assignment to `gxbeam_parameters` in soft scope is ambiguous because a global variable by the same name exists: `gxbeam_parameters` will be treated as a new local. Disambiguate by using `local gxbeam_parameters` to suppress this warning or `global gxbeam_parameters` to assign to the existing global variable.
└ @ cantilever-stability.md:153
┌ Warning: Assignment to `coupling_parameters` in soft scope is ambiguous because a global variable by the same name exists: `coupling_parameters` will be treated as a new local. Disambiguate by using `local coupling_parameters` to suppress this warning or `global coupling_parameters` to assign to the existing global variable.
└ @ cantilever-stability.md:156
┌ Warning: Assignment to `parameters` in soft scope is ambiguous because a global variable by the same name exists: `parameters` will be treated as a new local. Disambiguate by using `local parameters` to suppress this warning or `global parameters` to assign to the existing global variable.
└ @ cantilever-stability.md:161

To identify the flutter speed and frequency, we can plot the results.

using Plots
pyplot()

sp1 = plot(
    xlim = (0, 50),
    xtick = 0:5:50,
    xlabel = "Velocity (m/s)",
    ylim = (0, 350),
    ytick = 0:50:350,
    ylabel = "Frequency (rad/s)",
    framestyle = :zerolines,
    grid = false,
    titlefontsize = 14,
    guidefontsize = 14,
    legendfontsize = 11,
    tickfontsize = 11,
    legend = :topright,
    foreground_color_legend = nothing,
    background_color_legend = nothing,
    minorgrid=false)

sp2 = plot(
    xlim = (0, 50),
    xtick = 0:5:50,
    xlabel = "Velocity (m/s)",
    ylim = (-4, 2),
    ytick = -4:2:2,
    ylabel = "Damping (1/s)",
    framestyle = :zerolines,
    grid = false,
    titlefontsize = 14,
    guidefontsize = 14,
    legendfontsize = 11,
    tickfontsize = 11,
    legend = :topleft,
    foreground_color_legend = nothing,
    background_color_legend = nothing,
    minorgrid=false)

for i = 1:size(λ, 1)

    Vi = Vinf[:]
    λi = λ[i,:]

    if any(abs.(λi) .<= 1e4)
        plot!(sp1, Vi, imag.(λi),
            label = "",
            color = i,
            markersize = 3,
            markerstrokewidth = 0,
            )
    end

end

for i = 1:size(λ, 1)

    Vi = Vinf[:]
    λi = λ[i,:]

    if any(abs.(λi) .<= 1e4)
        plot!(sp2, Vi,
            real.(λi),#./sqrt.(real.(λi).^2 + abs.(λi).^2)*100,
            label = "",
            color = i,
            markersize = 3,
            markerstrokewidth = 0,
            )
    end
end

p1 = plot(sp1, sp2, layout = (1, 2), size = (800, 300))

nothing

These results are similar to those presented by Patil, Hodges, and Cesnik in "Nonlinear Aeroelasticity and Flight Dynamics of High Altitude Long-Endurance Aircraft" and therefore serve as a verification case for this coupled model.

This page was generated using Literate.jl.