# Flutter Analysis of a Goland Wing using the SHARPy Linear Solver¶

This is an example using SHARPy to find the flutter speed of a Goland wing by:

• Calculating aerodynamic forces and deflections using a nonlinear solver

• Creating a reduced order model of the linearised aerodynamics

• Evaluate the stability of the linearised aeroelastic system at different velocities

## References¶

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

### Required Packages¶

[1]:

import numpy as np
import matplotlib.pyplot as plt
import os
import sys
import 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¶

#### Velocity¶

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.

[2]:

u_inf = 1.
alpha_deg = 0.
rho = 1.02
num_modes = 4


#### Discretisation¶

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

[3]:

M = 16
N = 32
M_star_fact = 10


#### ROM¶

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: 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.

[4]:

c_ref = 1.8288 # Goland wing reference chord. Used for 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


[5]:

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_nmodes4rom_MIMORA_r6_sig0000_0000j
Case files will be saved in ./cases/goland_csM16N32Ms10_nmodes4rom_MIMORA_r6_sig0000_0000j
Output files will be saved in ./output/goland_csM16N32Ms10_nmodes4rom_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 cases.templates.flying_wings.GolandControlSurface for more info on the geometrical, structural and aerodynamic definition of the Goland wing here used.

[6]:

ws = wings.GolandControlSurface(M=M,
N=N,
Mstar_fact=M_star_fact,
u_inf=u_inf,
alpha=alpha_deg,
cs_deflection=[0, 0],
rho=rho,
sweep=0,
physical_time=2,
n_surfaces=2,
route=route_test_dir + '/cases',
case_name=case_name)

ws.clean_test_files()
ws.update_derived_params()
ws.set_default_config_dict()

ws.generate_aero_file()
ws.generate_fem_file()


#### 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.

[7]:

ws.config['SHARPy'] = {
'flow':
'StaticCoupled',
'AerogridPlot',
'BeamPlot',
'Modal',
'LinearAssembler',
'AsymptoticStability',
],
'case': ws.case_name, 'route': ws.route,
'write_screen': 'on', 'write_log': 'on',
'log_folder': route_test_dir + '/output/',
'log_file': ws.case_name + '.log'}


[8]:

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


[9]:

ws.config['AerogridLoader'] = {
'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¶

[10]:

ws.config['StaticCoupled'] = {
'print_info': 'on',
'max_iter': 200,
'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_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,
'delta_curved': 1e-1,
'min_delta': 1e-10,
'gravity_on': 'on',
'gravity': 9.81}}


## AerogridPlot Settings¶

[11]:

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


## BeamPlot Settings¶

[12]:

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


## Linear System Assembly Settings¶

[14]:

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¶

[15]:

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

[16]:

ws.config.write()


### Run SHARPy¶

[17]:

sharpy.sharpy_main.main(['', ws.route + ws.case_name + '.sharpy'])

--------------------------------------------------------------------------------
######  ##     ##    ###    ########  ########  ##    ##
##    ## ##     ##   ## ##   ##     ## ##     ##  ##  ##
##       ##     ##  ##   ##  ##     ## ##     ##   ####
######  ######### ##     ## ########  ########     ##
## ##     ## ######### ##   ##   ##           ##
##    ## ##     ## ##     ## ##    ##  ##           ##
######  ##     ## ##     ## ##     ## ##           ##
--------------------------------------------------------------------------------
Aeroelastics Lab, Aeronautics Department.
Running SHARPy from /home/ng213/2TB/pazy_code/pazy-sharpy/lib/sharpy/docs/source/content/example_notebooks
SHARPy being run is in /home/ng213/2TB/pazy_code/pazy-sharpy/lib/sharpy
The branch being run is dev_setting_error
The version and commit hash are: v1.2.1-339-g156c731-156c731
SHARPy output folder set
/home/ng213/2TB/pazy_code/pazy-sharpy/lib/sharpy/docs/source/content/example_notebooks/output//goland_csM16N32Ms10_nmodes4rom_MIMORA_r6_sig0000_0000j/
Variable for_pos has no assigned value in the settings file.
will default to the value: [0.0, 0, 0]
Variable control_surface_deflection has no assigned value in the settings file.
will default to the value: []
Variable control_surface_deflection_generator_settings has no assigned value in the settings file.
will default to the value: {}
Variable dx1 has no assigned value in the settings file.
will default to the value: -1.0
Variable ndx1 has no assigned value in the settings file.
will default to the value: 1
Variable r has no assigned value in the settings file.
will default to the value: 1.0
Variable dxmax has no assigned value in the settings file.
will default to the value: -1.0
The aerodynamic grid contains 2 surfaces
Surface 0, M=16, N=16
Wake 0, M=160, N=16
Surface 1, M=16, N=16
Wake 1, M=160, N=16
In total: 512 bound panels
In total: 5120 wake panels
Total number of panels = 5632
Generating an instance of StaticCoupled
Variable correct_forces_method has no assigned value in the settings file.
will default to the value:
Variable runtime_generators has no assigned value in the settings file.
will default to the value: {}
Generating an instance of NonLinearStatic
Variable newmark_damp has no assigned value in the settings file.
will default to the value: 0.0001
Variable gravity_dir has no assigned value in the settings file.
will default to the value: [0.0, 0.0, 1.0]
Variable relaxation_factor has no assigned value in the settings file.
will default to the value: 0.3
Variable dt has no assigned value in the settings file.
will default to the value: 0.01
Variable num_steps has no assigned value in the settings file.
will default to the value: 500
Variable initial_position has no assigned value in the settings file.
will default to the value: [0. 0. 0.]
Generating an instance of StaticUvlm
Variable iterative_solver has no assigned value in the settings file.
will default to the value: False
Variable iterative_tol has no assigned value in the settings file.
will default to the value: 0.0001
Variable iterative_precond has no assigned value in the settings file.
will default to the value: False
Variable cfl1 has no assigned value in the settings file.
will default to the value: True
Variable vortex_radius has no assigned value in the settings file.
will default to the value: 1e-06
Variable vortex_radius_wake_ind has no assigned value in the settings file.
will default to the value: 1e-06
Variable rbm_vel_g has no assigned value in the settings file.
will default to the value: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Variable centre_rot_g has no assigned value in the settings file.
will default to the value: [0.0, 0.0, 0.0]

|=====|=====|============|==========|==========|==========|==========|==========|==========|
|iter |step | log10(res) |    Fx    |    Fy    |    Fz    |    Mx    |    My    |    Mz    |
|=====|=====|============|==========|==========|==========|==========|==========|==========|
|  0  |  0  |  0.00000   |  0.0000  | -0.0000  |-4271.0417|  0.0000  | 781.0842 |  0.0000  |
|  1  |  0  | -11.89144  |  0.0000  | -0.0000  |-4271.0039|  0.0000  | 781.0906 |  0.0000  |
Generating an instance of AerogridPlot
Variable include_forward_motion has no assigned value in the settings file.
will default to the value: False
Variable include_unsteady_applied_forces has no assigned value in the settings file.
will default to the value: False
Variable name_prefix has no assigned value in the settings file.
will default to the value:
Variable u_inf has no assigned value in the settings file.
will default to the value: 0.0
Variable dt has no assigned value in the settings file.
will default to the value: 0.0
Variable include_velocities has no assigned value in the settings file.
will default to the value: False
Variable include_incidence_angle has no assigned value in the settings file.
will default to the value: False
Variable num_cores has no assigned value in the settings file.
will default to the value: 1
Variable vortex_radius has no assigned value in the settings file.
will default to the value: 1e-06
...Finished
Generating an instance of BeamPlot
Variable include_FoR has no assigned value in the settings file.
will default to the value: False
Variable include_applied_moments has no assigned value in the settings file.
will default to the value: True
Variable name_prefix has no assigned value in the settings file.
will default to the value:
Variable output_rbm has no assigned value in the settings file.
will default to the value: True
...Finished
Generating an instance of Modal
Variable print_info has no assigned value in the settings file.
will default to the value: True
Variable delta_curved has no assigned value in the settings file.
will default to the value: 0.01
Variable use_custom_timestep has no assigned value in the settings file.
will default to the value: -1
Structural eigenvalues

|==============|==============|==============|==============|==============|==============|==============|
|     mode     |  eval_real   |  eval_imag   | freq_n (Hz)  | freq_d (Hz)  |   damping    |  period (s)  |
|==============|==============|==============|==============|==============|==============|==============|
|      0       |   0.000000   |  48.067396   |   7.650164   |   7.650164   |  -0.000000   |   0.130716   |
|      1       |   0.000000   |  48.067398   |   7.650164   |   7.650164   |  -0.000000   |   0.130716   |
|      2       |   0.000000   |  95.685736   |  15.228858   |  15.228858   |  -0.000000   |   0.065665   |
|      3       |   0.000000   |  95.685754   |  15.228861   |  15.228861   |  -0.000000   |   0.065665   |
|      4       |   0.000000   |  243.144471  |  38.697644   |  38.697644   |  -0.000000   |   0.025841   |
|      5       |   0.000000   |  243.144477  |  38.697645   |  38.697645   |  -0.000000   |   0.025841   |
|      6       |   0.000000   |  343.801136  |  54.717650   |  54.717650   |  -0.000000   |   0.018276   |
|      7       |   0.000000   |  343.801137  |  54.717650   |  54.717650   |  -0.000000   |   0.018276   |
|      8       |   0.000000   |  443.324608  |  70.557303   |  70.557303   |  -0.000000   |   0.014173   |
|      9       |   0.000000   |  443.324619  |  70.557304   |  70.557304   |  -0.000000   |   0.014173   |
|      10      |   0.000000   |  461.992869  |  73.528449   |  73.528449   |  -0.000000   |   0.013600   |
|      11      |   0.000000   |  461.992869  |  73.528449   |  73.528449   |  -0.000000   |   0.013600   |
|      12      |   0.000000   |  601.126871  |  95.672313   |  95.672313   |  -0.000000   |   0.010452   |
|      13      |   0.000000   |  601.126873  |  95.672313   |  95.672313   |  -0.000000   |   0.010452   |
|      14      |   0.000000   |  782.997645  |  124.617946  |  124.617946  |  -0.000000   |   0.008025   |
|      15      |   0.000000   |  782.997649  |  124.617946  |  124.617946  |  -0.000000   |   0.008025   |
|      16      |   0.000000   |  917.191257  |  145.975522  |  145.975522  |  -0.000000   |   0.006850   |
|      17      |   0.000000   |  917.191259  |  145.975523  |  145.975523  |  -0.000000   |   0.006850   |
|      18      |   0.000000   |  975.005694  |  155.176976  |  155.176976  |  -0.000000   |   0.006444   |
|      19      |   0.000000   |  975.005699  |  155.176977  |  155.176977  |  -0.000000   |   0.006444   |
Generating an instance of LinearAssembler
Variable linearisation_tstep has no assigned value in the settings file.
will default to the value: -1
Variable modal_tstep has no assigned value in the settings file.
will default to the value: -1
Variable inout_coordinates has no assigned value in the settings file.
will default to the value:
Variable retain_inputs has no assigned value in the settings file.
will default to the value: []
Variable retain_outputs has no assigned value in the settings file.
will default to the value: []
Generating an instance of LinearAeroelastic
Variable uvlm_filename has no assigned value in the settings file.
will default to the value:
Variable track_body has no assigned value in the settings file.
will default to the value: True
Variable use_euler has no assigned value in the settings file.
will default to the value: False
Generating an instance of LinearUVLM
Variable gust_assembler has no assigned value in the settings file.
will default to the value:
Variable vortex_radius has no assigned value in the settings file.
will default to the value: 1e-06
Variable cfl1 has no assigned value in the settings file.
will default to the value: True
Variable velocity_field_generator has no assigned value in the settings file.
will default to the value: SteadyVelocityField
Variable velocity_field_input has no assigned value in the settings file.
will default to the value: {}
Variable physical_model has no assigned value in the settings file.
will default to the value: True
Variable track_body has no assigned value in the settings file.
will default to the value: False
Variable track_body_number has no assigned value in the settings file.
will default to the value: -1
Initialising Static linear UVLM solver class...
...done in 1.35 sec
Generating an instance of Krylov
Variable print_info has no assigned value in the settings file.
will default to the value: True
Variable single_side has no assigned value in the settings file.
will default to the value:
Variable tangent_input_file has no assigned value in the settings file.
will default to the value:
Variable restart_arnoldi has no assigned value in the settings file.
will default to the value: False
Initialising Krylov Model Order Reduction
State-space realisation of UVLM equations started...
Computing wake propagation matrix with CFL1=True

/home/ng213/anaconda3/envs/sharpy_env/lib/python3.7/site-packages/scipy/sparse/_index.py:126: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient.
self._set_arrayXarray(i, j, x)

     state-space model produced in form:
h_{n+1} = A h_{n} + B u_{n}
with:
x_n = h_n + Bp u_n
...done in 17.25 sec
Scaling UVLM system with reference time 0.914400s
Non-dimensional time step set (0.125000)
System scaled in 31.698308s
Generating an instance of LinearBeam
Warning, projecting system with damping onto undamped modes

Linearising gravity terms...
M = 7.26 kg
X_CG A -> 0.00 0.00 -0.00
Node  1      -> B -0.000 -0.116 0.000
-> A 0.116 0.381 -0.000
-> G 0.116 0.381 -0.000
Node mass:
Matrix: 14.5125
Node  2      -> B -0.000 -0.116 0.000
-> A 0.116 0.762 -0.000
-> G 0.116 0.762 -0.000
Node mass:
Matrix: 7.2563
Node  3      -> B -0.000 -0.116 0.000
-> A 0.116 1.143 -0.000
-> G 0.116 1.143 -0.000
Node mass:
Matrix: 14.5125
Node  4      -> B -0.000 -0.116 0.000
-> A 0.116 1.524 -0.001
-> G 0.116 1.524 -0.001
Node mass:
Matrix: 7.2563
Node  5      -> B -0.000 -0.116 0.000
-> A 0.116 1.905 -0.001
-> G 0.116 1.905 -0.001
Node mass:
Matrix: 14.5125
Node  6      -> B -0.000 -0.116 0.000
-> A 0.116 2.286 -0.001
-> G 0.116 2.286 -0.001
Node mass:
Matrix: 7.2563
Node  7      -> B -0.000 -0.116 0.000
-> A 0.116 2.667 -0.002
-> G 0.116 2.667 -0.002
Node mass:
Matrix: 14.5125
Node  8      -> B -0.000 -0.116 0.000
-> A 0.116 3.048 -0.002
-> G 0.116 3.048 -0.002
Node mass:
Matrix: 7.2563
Node  9      -> B -0.000 -0.116 0.000
-> A 0.116 3.429 -0.003
-> G 0.116 3.429 -0.003
Node mass:
Matrix: 14.5125
Node 10      -> B -0.000 -0.116 0.000
-> A 0.116 3.810 -0.003
-> G 0.116 3.810 -0.003
Node mass:
Matrix: 7.2563
Node 11      -> B -0.000 -0.116 0.000
-> A 0.116 4.191 -0.004
-> G 0.116 4.191 -0.004
Node mass:
Matrix: 14.5125
Node 12      -> B -0.000 -0.116 0.000
-> A 0.116 4.572 -0.004
-> G 0.116 4.572 -0.004
Node mass:
Matrix: 7.2563
Node 13      -> B -0.000 -0.116 0.000
-> A 0.116 4.953 -0.005
-> G 0.116 4.953 -0.005
Node mass:
Matrix: 14.5125
Node 14      -> B -0.000 -0.116 0.000
-> A 0.116 5.334 -0.005
-> G 0.116 5.334 -0.005
Node mass:
Matrix: 7.2563
Node 15      -> B -0.000 -0.116 0.000
-> A 0.116 5.715 -0.006
-> G 0.116 5.715 -0.006
Node mass:
Matrix: 14.5125
Node 16      -> B -0.000 -0.116 0.000
-> A 0.116 6.096 -0.006
-> G 0.116 6.096 -0.006
Node mass:
Matrix: 3.6281
Node 17      -> B -0.000 -0.116 -0.000
-> A 0.116 -6.096 -0.006
-> G 0.116 -6.096 -0.006
Node mass:
Matrix: 3.6281
Node 18      -> B -0.000 -0.116 -0.000
-> A 0.116 -5.715 -0.006
-> G 0.116 -5.715 -0.006
Node mass:
Matrix: 14.5125
Node 19      -> B -0.000 -0.116 -0.000
-> A 0.116 -5.334 -0.005
-> G 0.116 -5.334 -0.005
Node mass:
Matrix: 7.2563
Node 20      -> B -0.000 -0.116 -0.000
-> A 0.116 -4.953 -0.005
-> G 0.116 -4.953 -0.005
Node mass:
Matrix: 14.5125
Node 21      -> B -0.000 -0.116 -0.000
-> A 0.116 -4.572 -0.004
-> G 0.116 -4.572 -0.004
Node mass:
Matrix: 7.2563
Node 22      -> B -0.000 -0.116 -0.000
-> A 0.116 -4.191 -0.004
-> G 0.116 -4.191 -0.004
Node mass:
Matrix: 14.5125
Node 23      -> B -0.000 -0.116 -0.000
-> A 0.116 -3.810 -0.003
-> G 0.116 -3.810 -0.003
Node mass:
Matrix: 7.2563
Node 24      -> B -0.000 -0.116 -0.000
-> A 0.116 -3.429 -0.003
-> G 0.116 -3.429 -0.003
Node mass:
Matrix: 14.5125
Node 25      -> B -0.000 -0.116 -0.000
-> A 0.116 -3.048 -0.002
-> G 0.116 -3.048 -0.002
Node mass:
Matrix: 7.2563
Node 26      -> B -0.000 -0.116 -0.000
-> A 0.116 -2.667 -0.002
-> G 0.116 -2.667 -0.002
Node mass:
Matrix: 14.5125
Node 27      -> B -0.000 -0.116 -0.000
-> A 0.116 -2.286 -0.002
-> G 0.116 -2.286 -0.002
Node mass:
Matrix: 7.2563
Node 28      -> B -0.000 -0.116 -0.000
-> A 0.116 -1.905 -0.001
-> G 0.116 -1.905 -0.001
Node mass:
Matrix: 14.5125
Node 29      -> B -0.000 -0.116 -0.000
-> A 0.116 -1.524 -0.001
-> G 0.116 -1.524 -0.001
Node mass:
Matrix: 7.2563
Node 30      -> B -0.000 -0.116 -0.000
-> A 0.116 -1.143 -0.000
-> G 0.116 -1.143 -0.000
Node mass:
Matrix: 14.5125
Node 31      -> B -0.000 -0.116 -0.000
-> A 0.116 -0.762 -0.000
-> G 0.116 -0.762 -0.000
Node mass:
Matrix: 7.2563
Node 32      -> B -0.000 -0.116 -0.000
-> A 0.116 -0.381 -0.000
-> G 0.116 -0.381 -0.000
Node mass:
Matrix: 14.5125
Updated the beam C, modal C and K matrices with the terms from the
gravity linearisation

Scaling beam according to reduced time...
Setting the beam time step to (0.1250)
Updating C and K matrices and natural frequencies with new normalised time...
Model Order Reduction in progress...
Moment Matching Krylov Model Reduction
Construction Algorithm:
mimo_rational_arnoldi
Interpolation points:
sigma = 0.000000 + 0.000000j [rad/s]

Krylov order:
r = 6
Constructing controllability space
Constructing observability space
Deflating column 23
Deflating column 25
Deflating column 28
Deflating column 33
ROM is stable
DT Eigenvalues:
mu = 0.992076 + -0.000000j
mu = 0.971409 + -0.000000j
mu = 0.957534 + -0.038927j
mu = 0.957534 + 0.038927j
mu = 0.954926 + -0.072820j
mu = 0.954926 + 0.072820j
mu = 0.954282 + -0.000000j
mu = 0.935272 + 0.000000j
mu = 0.927865 + 0.092012j
mu = 0.927865 + -0.092012j
mu = 0.927075 + -0.066636j
mu = 0.927075 + 0.066636j
mu = 0.925588 + 0.038813j
mu = 0.925588 + -0.038813j
mu = 0.922836 + -0.004482j
mu = 0.922836 + 0.004482j
mu = 0.915491 + -0.025196j
mu = 0.915491 + 0.025196j
mu = 0.908662 + 0.053066j
mu = 0.908662 + -0.053066j
mu = 0.881127 + -0.056743j
mu = 0.881127 + 0.056743j
mu = 0.882642 + -0.000000j
mu = 0.867726 + -0.000000j
mu = 0.717587 + 0.263893j
mu = 0.717587 + -0.263893j
mu = 0.697184 + 0.000000j
mu = 0.360807 + 0.000000j
mu = -0.000018 + -0.002663j
mu = -0.000018 + 0.002663j
mu = 0.000155 + -0.000000j
mu = -0.000154 + 0.000000j

System reduced from order 6656 to
n = 32 states
...Completed Model Order Reduction in 5.16 s
Aeroelastic system assembled:
Aerodynamic states: 32
Structural states: 4
Total states: 36
Inputs: 8
Outputs: 6
Final system is:
State-space system
States: 36
Inputs: 8
Outputs: 6

Generating an instance of AsymptoticStability
Variable reference_velocity has no assigned value in the settings file.
will default to the value: 1.0
Variable frequency_cutoff has no assigned value in the settings file.
will default to the value: 0.0
Variable export_eigenvalues has no assigned value in the settings file.
will default to the value: False
Variable display_root_locus has no assigned value in the settings file.
will default to the value: False
Variable iterative_eigvals has no assigned value in the settings file.
will default to the value: False
Variable num_evals has no assigned value in the settings file.
will default to the value: 200
Variable postprocessors has no assigned value in the settings file.
will default to the value: []
Variable postprocessors_settings has no assigned value in the settings file.
will default to the value: {}
Dynamical System Eigenvalues
Calculating eigenvalues using direct method

|==============|==============|==============|==============|==============|==============|==============|
|     mode     |  eval_real   |  eval_imag   | freq_n (Hz)  | freq_d (Hz)  |   damping    |  period (s)  |
|==============|==============|==============|==============|==============|==============|==============|
|      0       |  -0.021637   |  24.315428   |   3.869922   |   3.869921   |   0.000890   |   0.258403   |
|      1       |  -0.021637   |  -24.315428  |   3.869922   |   3.869921   |   0.000890   |   0.258403   |
|      2       |  -0.069601   |  -0.000000   |   0.011077   |   0.000000   |   1.000000   |     inf      |
|      3       |  -0.105223   |  -21.320535  |   3.393310   |   3.393269   |   0.004935   |   0.294701   |
|      4       |  -0.105223   |  21.320535   |   3.393310   |   3.393269   |   0.004935   |   0.294701   |
|      5       |  -0.253787   |  -0.000000   |   0.040391   |   0.000000   |   1.000000   |     inf      |
|      6       |  -0.372425   |   0.355473   |   0.081940   |   0.056575   |   0.723379   |  17.675568   |
|      7       |  -0.372425   |  -0.355473   |   0.081940   |   0.056575   |   0.723379   |  17.675568   |
|      8       |  -0.378143   |   0.665882   |   0.121875   |   0.105978   |   0.493812   |   9.435879   |
|      9       |  -0.378143   |  -0.665882   |   0.121875   |   0.105978   |   0.493812   |   9.435879   |
|      10      |  -0.409413   |  -0.000000   |   0.065160   |   0.000000   |   1.000000   |     inf      |
|      11      |  -0.585461   |   0.000000   |   0.093179   |   0.000000   |   1.000000   |     inf      |
|      12      |  -0.612214   |  -0.864763   |   0.168631   |   0.137631   |   0.577812   |   7.265789   |
|      13      |  -0.612214   |   0.864763   |   0.168631   |   0.137631   |   0.577812   |   7.265789   |
|      14      |  -0.639935   |   0.627769   |   0.142673   |   0.099913   |   0.713860   |  10.008752   |
|      15      |  -0.639935   |  -0.627769   |   0.142673   |   0.099913   |   0.713860   |  10.008752   |
|      16      |  -0.668833   |   0.366654   |   0.121394   |   0.058355   |   0.876881   |  17.136529   |
|      17      |  -0.668833   |  -0.366654   |   0.121394   |   0.058355   |   0.876881   |  17.136529   |
|      18      |  -0.702465   |   0.042487   |   0.112005   |   0.006762   |   0.998176   |  147.883816  |
|      19      |  -0.702465   |  -0.042487   |   0.112005   |   0.006762   |   0.998176   |  147.883816  |
|      20      |  -0.769174   |   0.240726   |   0.128273   |   0.038313   |   0.954353   |  26.100974   |
|      21      |  -0.769174   |  -0.240726   |   0.128273   |   0.038313   |   0.954353   |  26.100974   |
|      22      |  -0.823096   |   0.510355   |   0.154138   |   0.081226   |   0.849886   |  12.311395   |
|      23      |  -0.823096   |  -0.510355   |   0.154138   |   0.081226   |   0.849886   |  12.311395   |
|      24      |  -1.089098   |   0.562639   |   0.195099   |   0.089547   |   0.888447   |  11.167353   |
|      25      |  -1.089098   |  -0.562639   |   0.195099   |   0.089547   |   0.888447   |  11.167353   |
|      26      |  -1.092180   |  -0.000000   |   0.173826   |   0.000000   |   1.000000   |     inf      |
|      27      |  -1.241288   |   0.000000   |   0.197557   |   0.000000   |   1.000000   |     inf      |
|      28      |  -2.348569   |   3.083089   |   0.616840   |   0.490689   |   0.605970   |   2.037951   |
|      29      |  -2.348569   |  -3.083089   |   0.616840   |   0.490689   |   0.605970   |   2.037951   |
|      30      |  -3.155780   |  -0.000000   |   0.502258   |   0.000000   |   1.000000   |     inf      |
|      31      |  -8.918115   |  -0.000000   |   1.419362   |   0.000000   |   1.000000   |     inf      |
|      32      |  -26.700406  |  -27.485500  |   6.098697   |   4.374453   |   0.696788   |   0.228600   |
|      33      |  -28.824123  |   0.000000   |   4.587502   |   0.000000   |   1.000000   |     inf      |
|      34      |  -35.766803  |  -27.485500  |   7.179135   |   4.374453   |   0.792918   |   0.228600   |
|      35      |  -36.676675  |  -0.000000   |   5.837274   |   0.000000   |   1.000000   |     inf      |
Velocity Asymptotic Stability Analysis
Initial velocity: 100.00 m/s
Final velocity: 180.00 m/s
Number of evaluations: 81
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 100.00 m/2   max. CT eig. real: -3.925687
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 101.00 m/2   max. CT eig. real: -3.951244
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 102.00 m/2   max. CT eig. real: -3.975951
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 103.00 m/2   max. CT eig. real: -3.999782
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 104.00 m/2   max. CT eig. real: -4.022709
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 105.00 m/2   max. CT eig. real: -4.044705
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 106.00 m/2   max. CT eig. real: -4.065744
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 107.00 m/2   max. CT eig. real: -4.085798
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 108.00 m/2   max. CT eig. real: -4.104841
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 109.00 m/2   max. CT eig. real: -4.122845
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 110.00 m/2   max. CT eig. real: -4.139781
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 111.00 m/2   max. CT eig. real: -4.155617
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 112.00 m/2   max. CT eig. real: -4.170322
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 113.00 m/2   max. CT eig. real: -4.183860
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 114.00 m/2   max. CT eig. real: -4.196192
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 115.00 m/2   max. CT eig. real: -4.207276
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 116.00 m/2   max. CT eig. real: -4.217065
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 117.00 m/2   max. CT eig. real: -4.225507
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 118.00 m/2   max. CT eig. real: -4.232544
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 119.00 m/2   max. CT eig. real: -4.238111
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 120.00 m/2   max. CT eig. real: -4.242138
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 121.00 m/2   max. CT eig. real: -4.244545
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 122.00 m/2   max. CT eig. real: -4.245246
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 123.00 m/2   max. CT eig. real: -4.244143
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 124.00 m/2   max. CT eig. real: -4.241130
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 125.00 m/2   max. CT eig. real: -4.236092
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 126.00 m/2   max. CT eig. real: -4.228899
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 127.00 m/2   max. CT eig. real: -4.219413
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 128.00 m/2   max. CT eig. real: -4.207482
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 129.00 m/2   max. CT eig. real: -4.192940
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 130.00 m/2   max. CT eig. real: -4.175607
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 131.00 m/2   max. CT eig. real: -4.155291
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 132.00 m/2   max. CT eig. real: -4.131780
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 133.00 m/2   max. CT eig. real: -4.104848
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 134.00 m/2   max. CT eig. real: -4.074252
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 135.00 m/2   max. CT eig. real: -4.039730
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 136.00 m/2   max. CT eig. real: -4.001000
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 137.00 m/2   max. CT eig. real: -3.957763
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 138.00 m/2   max. CT eig. real: -3.909697
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 139.00 m/2   max. CT eig. real: -3.856462
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 140.00 m/2   max. CT eig. real: -3.797697
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 141.00 m/2   max. CT eig. real: -3.733020
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 142.00 m/2   max. CT eig. real: -3.662031
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 143.00 m/2   max. CT eig. real: -3.584314
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 144.00 m/2   max. CT eig. real: -3.499436
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 145.00 m/2   max. CT eig. real: -3.406959
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 146.00 m/2   max. CT eig. real: -3.306437
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 147.00 m/2   max. CT eig. real: -3.197433
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 148.00 m/2   max. CT eig. real: -3.079522
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 149.00 m/2   max. CT eig. real: -2.952307
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 150.00 m/2   max. CT eig. real: -2.815436
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 151.00 m/2   max. CT eig. real: -2.668615
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 152.00 m/2   max. CT eig. real: -2.511634
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 153.00 m/2   max. CT eig. real: -2.344382
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 154.00 m/2   max. CT eig. real: -2.166866
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 155.00 m/2   max. CT eig. real: -1.979231
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 156.00 m/2   max. CT eig. real: -1.781770
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 157.00 m/2   max. CT eig. real: -1.574929
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 158.00 m/2   max. CT eig. real: -1.359301
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 159.00 m/2   max. CT eig. real: -1.135613
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 160.00 m/2   max. CT eig. real: -0.904707
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 161.00 m/2   max. CT eig. real: -0.667505
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 162.00 m/2   max. CT eig. real: -0.424982
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 163.00 m/2   max. CT eig. real: -0.178127
N unstab.: 000
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 164.00 m/2   max. CT eig. real: 0.072086
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       70.64   70.64
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 165.00 m/2   max. CT eig. real: 0.324729
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       70.41   70.41
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 166.00 m/2   max. CT eig. real: 0.578936
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       70.20   70.20
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 167.00 m/2   max. CT eig. real: 0.833925
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       69.99   69.99
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 168.00 m/2   max. CT eig. real: 1.088998
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       69.79   69.79
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 169.00 m/2   max. CT eig. real: 1.343552
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       69.61   69.61
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 170.00 m/2   max. CT eig. real: 1.597068
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       69.43   69.43
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 171.00 m/2   max. CT eig. real: 1.849117
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       69.26   69.26
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 172.00 m/2   max. CT eig. real: 2.099343
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       69.10   69.10
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 173.00 m/2   max. CT eig. real: 2.347461
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       68.94   68.94
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 174.00 m/2   max. CT eig. real: 2.593246
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       68.79   68.79
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 175.00 m/2   max. CT eig. real: 2.836529
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       68.65   68.65
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 176.00 m/2   max. CT eig. real: 3.077180
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       68.51   68.51
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 177.00 m/2   max. CT eig. real: 3.315113
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       68.38   68.38
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 178.00 m/2   max. CT eig. real: 3.550268
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       68.25   68.25
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 179.00 m/2   max. CT eig. real: 3.782616
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       68.13   68.13
Updating C and K matrices and natural frequencies with new normalised time...
LTI     u: 180.00 m/2   max. CT eig. real: 4.012145
N unstab.: 002
Unstable aeroelastic natural frequency CT(rad/s):       68.01   68.01
Saving velocity analysis results...
Successful
FINISHED - Elapsed time = 59.1569176 seconds
FINISHED - CPU process time = 94.3774846 seconds

[17]:

<sharpy.presharpy.presharpy.PreSharpy at 0x7f85b5231850>


### Analysis¶

#### 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

#### Stability¶

The stability of the Goland wing is now analysed under changing free stream velocity. The aeroelastic system is projected onto 2 structural modes (1st bending and 1st torsion). 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. This flutter mode is a bending-torsion mode, as seen from the natural frequency plot where the frequencies of each coalesce into this mode.

[18]:

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

u_inf = velocity_analysis[:, 0]
eigs_r = velocity_analysis[:, 1]
eigs_i = velocity_analysis[:, 2]

[19]:

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.grid()
plt.xlim(-10, 10)
plt.ylim(-150, 150)
plt.xlabel('Real Part, $\lambda$ [rad/s]')
plt.ylabel('Imag Part, $\lambda$ [rad/s]');

[20]:

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.grid()
plt.ylim(-0.25, 0.25)
plt.xlabel('Free Stream Velocity, $u_\infty$ [m/s]')
plt.ylabel('Damping Ratio, $\zeta$ [-]');

[21]:

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.grid()
plt.ylim(40, 100)
plt.xlabel('Free Stream Velocity, $u_\infty$ [m/s]')
plt.ylabel('Natural Frequency, $\omega_n$ [rad/s]');