Flutter Analysis of Goland’s Wing


This is simple rectangular wing with constant properties originally proposed by Martin Goland in 1945 as a flutter benchmark case. This example generates a beam/UVLM model with the original properties in SHARPy and computes flutter around an arbitrary equilibrium point using the following general process:

  • Calculate steady aerodynamic forces and deflections for a given angle of attack using a nonlinear solver

  • Linearise the dynamic aeroelastic equations about this reference condition

  • Reduce the size of the (linearized) structural subsystem using modal projection

  • Reduce the size of the (linearized) aerodyanic subsystem using a Krylov projection

  • Evaluate the stability of the linearised aeroelastic system at different velocities and plot the results.

Wing properties

Span: b=20 ft., chord: c=6 ft., weight: 4 psi., radius of gyration about CG: .25c, elastic axis: 0.33c, centre of gravity: 0.43c


Goland, M. (1945). The Flutter o a Uniorm Cantilever Wing. Journal of Applied Mechanics,12(4):197-208.

Maraniello, S., & Palacios, R. (2019). State-Space Realizations and Internal Balancing in Potential-Flow Aerodynamics with Arbitrary Kinematics. AIAA Journal, 57(6):1–14. https://doi.org/10.2514/1.J058153

Latest update, RPN 14.05.2023

Required Packages

import numpy as np
import matplotlib.pyplot as plt
import os
import sys
import sharpy.cases.templates.flying_wings as wings  # See this package for the Goland wing structural and aerodynamic definition
import sharpy.sharpy_main  # used to run SHARPy from Jupyter

Problem Set-up

Operatig conditions

The UVLM is assembled in normalised time at a velocity of \(1 m/s\). The only matrices that need updating then with free stream velocity are the structural matrices, which is significantly cheaper to do than to update the UVLM.

u_inf = 1.
alpha_deg = 2.  # Define angle of attack for static aeroelastic analysis.
rho = 1.02      # Air density.


Note: To achieve convergence of the flutter results with the ones found in the literature, a fine grid is needed. If you are running this notebook for the first time, set M = 4 initially to verify that your system can perform!

M = 16             # Number of chordwise panels
N = 32             # Number of spanwise panels
M_star_fact = 10   # Length of the wake in chords.
num_modes =  8     # Number of vibration modes retained in the structural model.


A moment-matching (Krylov subspace) model order reduction technique is employed. This ROM method offers the ability to interpolate the transfer functions at a desired point in the complex plane. See the ROM documentation pages for more info. Note that this ROM method matches the transfer function but does not guarantee stability. Therefore the resulting system may be unstable. These unstable modes may appear far in the right hand plane but will not affect the flutter speed calculations.

c_ref = 1.8288 # Goland wing reference chord. Used for reduced frequency normalisation
rom_settings = dict()
rom_settings['algorithm'] = 'mimo_rational_arnoldi'  # reduction algorithm
rom_settings['r'] = 6  # Krylov subspace order
frequency_continuous_k = np.array([0.])  # Interpolation point in the complex plane with reduced frequency units
frequency_continuous_w = 2 * u_inf * frequency_continuous_k / c_ref
rom_settings['frequency'] = frequency_continuous_w

Case Admin

case_name = 'goland_cs'
case_nlin_info = 'M%dN%dMs%d_nmodes%d' % (M, N, M_star_fact, num_modes)
case_rom_info = 'rom_MIMORA_r%d_sig%04d_%04dj' % (rom_settings['r'], frequency_continuous_k[-1].real * 100,
                                                  frequency_continuous_k[-1].imag * 100)

case_name += case_nlin_info + case_rom_info
route_test_dir = os.path.abspath('')

print('The case to run will be: %s' % case_name)
print('Case files will be saved in ./cases/%s' %case_name)
print('Output files will be saved in ./output/%s/' %case_name)
The case to run will be: goland_csM16N32Ms10_nmodes8rom_MIMORA_r6_sig0000_0000j
Case files will be saved in ./cases/goland_csM16N32Ms10_nmodes8rom_MIMORA_r6_sig0000_0000j
Output files will be saved in ./output/goland_csM16N32Ms10_nmodes8rom_MIMORA_r6_sig0000_0000j/

Simulation Set-Up

Goland Wing

ws is an instance of a Goland wing with a control surface. Reference the template file sharpy.cases.templates.flying_wings.GolandControlSurface for more info on the geometrical, structural and aerodynamic definition of the Goland wing here used.

ws = wings.GolandControlSurface(M=M,
                                cs_deflection=[0, 0],
                                route=route_test_dir + '/cases',



Simulation Settings

The settings for each of the solvers are now set. For a detailed description on them please reference their respective documentation pages

SHARPy Settings

The most important setting is the flow list. It tells SHARPy which solvers to run and in which order.

ws.config['SHARPy'] = {
        ['BeamLoader', 'AerogridLoader',
    'case': ws.case_name, 'route': ws.route,
    'write_screen': 'off', 'write_log': 'on',    # Change to 'on' as neded.
    'log_folder': route_test_dir + '/output/',
    'log_file': ws.case_name + '.log'}

Beam Loader Settings

ws.config['BeamLoader'] = {
    'unsteady': 'off',
    'orientation': ws.quat}

Aerogrid Loader Settings

ws.config['AerogridLoader'] = {
    'unsteady': 'off',
    'aligned_grid': 'on',
    'mstar': ws.Mstar_fact * ws.M,
    'freestream_dir': ws.u_inf_direction,
    'wake_shape_generator': 'StraightWake',
    'wake_shape_generator_input': {'u_inf': ws.u_inf,
                                    'u_inf_direction': ws.u_inf_direction,
                                    'dt': ws.dt}}

Static Coupled Solver

ws.config['StaticCoupled'] = {
    'print_info': 'on',
    'max_iter': 200,
    'n_load_steps': 1,
    'tolerance': 1e-10,
    'relaxation_factor': 0.,
    'aero_solver': 'StaticUvlm',
    'aero_solver_settings': {
        'rho': ws.rho,
        'print_info': 'off',
        'horseshoe': 'off',
        'num_cores': 4,
        'n_rollup': 0,
        'rollup_dt': ws.dt,
        'rollup_aic_refresh': 1,
        'rollup_tolerance': 1e-4,
        'velocity_field_generator': 'SteadyVelocityField',
        'velocity_field_input': {
            'u_inf': ws.u_inf,
            'u_inf_direction': ws.u_inf_direction}},
    'structural_solver': 'NonLinearStatic',
    'structural_solver_settings': {'print_info': 'off',
                                   'max_iterations': 150,
                                   'num_load_steps': 4,
                                   'delta_curved': 1e-1,
                                   'min_delta': 1e-10,
                                   'gravity_on': 'on',
                                   'gravity': 9.81}}

AerogridPlot Settings

ws.config['AerogridPlot'] = {'include_rbm': 'off',
                             'include_applied_forces': 'on',
                             'minus_m_star': 0}

BeamPlot Settings

ws.config['BeamPlot'] = {'include_rbm': 'off',
                         'include_applied_forces': 'on'}

Linear System Assembly Settings

ws.config['LinearAssembler'] = {'linear_system': 'LinearAeroelastic',
                                'linear_system_settings': {
                                    'beam_settings': {'modal_projection': 'on',
                                                      'inout_coords': 'modes',
                                                      'discrete_time': 'on',
                                                      'newmark_damp': 0.5e-4,
                                                      'discr_method': 'newmark',
                                                      'dt': ws.dt,
                                                      'proj_modes': 'undamped',
                                                      'use_euler': 'off',
                                                      'num_modes': num_modes,
                                                      'print_info': 'on',
                                                      'gravity': 'on',
                                                      'remove_sym_modes': 'on',
                                                      'remove_dofs': []},
                                    'aero_settings': {'dt': ws.dt,
                                                      'ScalingDict': {'length': 0.5 * ws.c_ref,
                                                                      'speed': u_inf,
                                                                      'density': rho},
                                                      'integr_order': 2,
                                                      'density': ws.rho,
                                                      'remove_predictor': 'on',
                                                      'use_sparse': 'on',
                                                      'remove_inputs': ['u_gust'],
                                                      'rom_method': ['Krylov'],
                                                      'rom_method_settings': {'Krylov': rom_settings}},

Asymptotic Stability Analysis Settings

ws.config['AsymptoticStability'] = {'print_info': True,
                                    'velocity_analysis': [100, 180, 81],
                                   'modes_to_plot': []}

Write solver settings config file



sharpy.sharpy_main.main(['', ws.route + ws.case_name + '.sharpy'])
fatal: not a git repository (or any of the parent directories): .git
/home/rpalacio/anaconda3/envs/sharpy/lib/python3.10/site-packages/scipy/sparse/_index.py:146: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient.
  self._set_arrayXarray(i, j, x)
/home/rpalacio/anaconda3/envs/sharpy/lib/python3.10/site-packages/sharpy/linear/src/lingebm.py:313: UserWarning: Euler parametrisation not implemented - Either rigid body modes are not being used or this method has already been called.
  warnings.warn('Euler parametrisation not implemented - Either rigid body modes are not being used or this '
/home/rpalacio/anaconda3/envs/sharpy/lib/python3.10/site-packages/sharpy/rom/krylov.py:242: UserWarning: Reduced Order Model Unstable
  warn.warn('Reduced Order Model Unstable')
<sharpy.presharpy.presharpy.PreSharpy at 0x7f90aedb16c0>


Nonlinear equilibrium

The nonlinear equilibrium condition can be visualised and analysed by opening, with Paraview, the files in the /output/<case_name>/aero and /output/<case_name>/beam folders to see the deflection and aerodynamic forces acting on the wing.


The stability of the Goland wing is now analysed under changing free stream velocity. The flutter modes involves the two lowest frequency modes near the imaginary axis (1st bending and 1st torsion if aerodynamics is removed). The two modes are seen quite separated at 100 m/s. As speed is increased, the damping of the torsion mode decreases until it crosses the imaginary axis onto the right hand plane and flutter begins.

file_name = './output/%s/stability/velocity_analysis_min1000_max1800_nvel0081.dat' % case_name

velocity_analysis = np.loadtxt(file_name)
u_inf = velocity_analysis[:, 0]
eigs_r = velocity_analysis[:, 1]
eigs_i = velocity_analysis[:, 2]
fig = plt.figure()
plt.scatter(eigs_r, eigs_i, c=u_inf, cmap='Blues')
cbar = plt.colorbar()
cbar.set_label('Free Stream Velocity, $u_\infty$ [m/s]')

plt.xlim(-50, 5)
plt.ylim(0, 300)
plt.xlabel('Real Part, $\lambda$ [rad/s]')
plt.ylabel('Imag Part, $\lambda$ [rad/s]');
fig = plt.figure()
natural_frequency = np.sqrt(eigs_r ** 2 + eigs_i ** 2)
damping_ratio = eigs_r / natural_frequency
cond = (eigs_r>-25) * (eigs_r<10) * (natural_frequency<100) # filter unwanted eigenvalues for this plot (mostly aero modes)

plt.scatter(u_inf[cond], damping_ratio[cond], color='k', marker='s', s=9)

plt.ylim(-0.25, 0.25)
plt.xlabel('Free Stream Velocity, $u_\infty$ [m/s]')
plt.ylabel('Damping Ratio, $\zeta$ [-]');
fig = plt.figure()
cond = (eigs_r>-25) * (eigs_r<10) # filter unwanted eigenvalues for this plot (mostly aero modes)
plt.scatter(u_inf[cond], natural_frequency[cond], color='k', marker='s', s=9)

plt.ylim(40, 100)
plt.xlabel('Free Stream Velocity, $u_\infty$ [m/s]')
plt.ylabel('Natural Frequency, $\omega_n$ [rad/s]');