<jupyter_start><jupyter_text>Flutter Analysis of a Goland Wing using the SHARPy Linear SolverThis is an example using SHARPy to find the flutter speed of a Goland wing by:* Calculating aerodynamic forces and deflections using a nonlinear solver* Linearising about this reference condition* Creating a reduced order model of the linearised aerodynamics* Evaluate the stability of the linearised aeroelastic system at different velocities  References<NAME>., & <NAME>. (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<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Problem Set-up  VelocityThe 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.<jupyter_code>u_inf = 1.
alpha_deg = 0.
rho = 1.02
num_modes = 4<jupyter_output><empty_output><jupyter_text>DiscretisationNote: 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!<jupyter_code>M = 16
N = 32
M_star_fact = 10<jupyter_output><empty_output><jupyter_text>ROMA 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.<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Case Admin<jupyter_code>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)<jupyter_output>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/<jupyter_text>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.<jupyter_code>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()<jupyter_output>Surface0
Surface1<jupyter_text>Simulation SettingsThe settings for each of the solvers are now set. For a detailed description on them please reference their respective documentation pages  SHARPy SettingsThe most important setting is the `flow` list. It tells SHARPy which solvers to run and in which order.<jupyter_code>ws.config['SHARPy'] = {
    'flow':
        ['BeamLoader', 'AerogridLoader',
         'StaticCoupled',
         'AerogridPlot',
         'BeamPlot',
         'Modal',
         'LinearAssembler',
         'FrequencyResponse',
         'AsymptoticStability',
         ],
    'case': ws.case_name, 'route': ws.route,
    'write_screen': 'on', 'write_log': 'on',
    'log_folder': route_test_dir + '/output/' + ws.case_name + '/',
    'log_file': ws.case_name + '.log'}<jupyter_output><empty_output><jupyter_text>Beam Loader Settings<jupyter_code>ws.config['BeamLoader'] = {
    'unsteady': 'off',
    'orientation': ws.quat}<jupyter_output><empty_output><jupyter_text>Aerogrid Loader Settings<jupyter_code>ws.config['AerogridLoader'] = {
    'unsteady': 'off',
    'aligned_grid': 'on',
    'mstar': ws.Mstar_fact * ws.M,
    'freestream_dir': ws.u_inf_direction
}<jupyter_output><empty_output><jupyter_text>Static Coupled Solver<jupyter_code>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}}<jupyter_output><empty_output><jupyter_text>AerogridPlot Settings<jupyter_code>ws.config['AerogridPlot'] = {'folder': route_test_dir + '/output/',
                             'include_rbm': 'off',
                             'include_applied_forces': 'on',
                             'minus_m_star': 0}<jupyter_output><empty_output><jupyter_text>BeamPlot Settings<jupyter_code>ws.config['BeamPlot'] = {'folder': route_test_dir + '/output/',
                         'include_rbm': 'off',
                         'include_applied_forces': 'on'}<jupyter_output><empty_output><jupyter_text>Modal Solver Settings<jupyter_code>ws.config['Modal'] = {'folder': route_test_dir + '/output/',
                      'NumLambda': 20,
                      'rigid_body_modes': 'off',
                      'print_matrices': 'on',
                      'keep_linear_matrices': 'on',
                      'write_dat': 'off',
                      'rigid_modes_cg': 'off',
                      'continuous_eigenvalues': 'off',
                      'dt': 0,
                      'plot_eigenvalues': False,
                      'max_rotation_deg': 15.,
                      'max_displacement': 0.15,
                      'write_modes_vtk': True,
                      'use_undamped_modes': True}<jupyter_output><empty_output><jupyter_text>Linear System Assembly Settings<jupyter_code>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',
                                                      'rigid_body_motion': 'off',
                                                      'use_euler': 'off',
                                                      'remove_inputs': ['u_gust'],
                                                      'rom_method': ['Krylov'],
                                                      'rom_method_settings': {'Krylov': rom_settings}},
                                    'rigid_body_motion': False}}<jupyter_output><empty_output><jupyter_text>Asymptotic Stability Analysis Settings<jupyter_code>ws.config['AsymptoticStability'] = {'print_info': True,
                                    'folder': route_test_dir + '/output/',
                                    'velocity_analysis': [100, 180, 81],
                                   'modes_to_plot': []}<jupyter_output><empty_output><jupyter_text>Write solver settings config file<jupyter_code>ws.config.write()<jupyter_output><empty_output><jupyter_text>Run SHARPy<jupyter_code>sharpy.sharpy_main.main(['', ws.route + ws.case_name + '.solver.txt'])<jupyter_output>--------------------------------------------------------------------------------[0m
            ######  ##     ##    ###    ########  ########  ##    ##[0m
           ##    ## ##     ##   ## ##   ##     ## ##     ##  ##  ##[0m
           ##       ##     ##  ##   ##  ##     ## ##     ##   ####[0m
            ######  ######### ##     ## ########  ########     ##[0m
                 ## ##     ## ######### ##   ##   ##           ##[0m
           ##    ## ##     ## ##     ## ##    ##  ##           ##[0m
            ######  ##     ## ##     ## ##     ## ##           ##[0m
--------------------------------------------------------------------------------[0m
Aeroelastics Lab, Aeronautics Department.[0m
    Copyright (c), Imperial College London.[0m
    All rights reserved. [0m
    License available at https://github.com/imperialcollegelondon/sharpy[0m
[36mRunning SHARPy from /home/ng213/code/sharpy/docs/source/content/example_notebooks[0m
[36mSHARPy being run is in /home/ng213/co[...]<jupyter_text>Analysis  Nonlinear equilibriumThe nonlinear equilibrium condition can be visualised and analysed by opening, with Paraview, the files in the `/output//aero` and `/output//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.<jupyter_code>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.grid()
plt.xlim(-10, 10)
plt.ylim(-150, 150)
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.grid()
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.grid()
plt.ylim(40, 100)
plt.xlabel('Free Stream Velocity, $u_\infty$ [m/s]')
plt.ylabel('Natural Frequency, $\omega_n$ [rad/s]');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12267997750281215, 'fraction_numerical': 0.016943194600674914, 'mean_word_length': 2.0703647744442044, 'pattern_counts': {'":': 0, '<': 75, '<?xml version=': 0, '>': 74, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 69, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '197480', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4368, 'n_tokens_neox': 3996, 'n_words': 967, 'path': 'docs/source/content/example_notebooks/linear_goland_flutter.ipynb', 'repo_name': 'LMNS3d/sharpy', 'size': '14236'}