{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# A first test case with SHARPy\n", "This notebook explains how to create a prismatic cantilever and flexible wing in SHARPy. It aims to provide a big picture about the simulations available in SHARPy describing the very fundamental concepts. \n", "\n", "This notebook requires around two hours to be completed." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generation of the aeroelastic model\n", "This section explains how to generate the input files for SHARPy including the structure, the aerodynamics and the simulation details." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Loading of the used packages\n", "import numpy as np # basic mathematical and array functions\n", "import os # Functions related to the operating system\n", "import matplotlib.pyplot as plt # Plotting library\n", "\n", "import sharpy.sharpy_main # Run SHARPy inside jupyter notebooks\n", "import sharpy.utils.plotutils as pu # Plotting utilities\n", "from sharpy.utils.constants import deg2rad # Constant to conver degrees to radians\n", "\n", "import sharpy.utils.generate_cases as gc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The [generate cases](https://ic-sharpy.readthedocs.io/en/master/includes/utils/generate_cases/index.html) module of SHARPy includes a number of templates for the creation of simple cases using common parameters as inputs. It also removes clutter from this notebook.\n", "\n", "The following cell configures the plotting library to show the next graphics, it is not part of SHARPy." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "! pip install plotly" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define the problem parameters\n", "In this section, we define the basic parameters to generate a model for the wing shown in the image below. We use SI units for all variables." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ " " ], "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%config InlineBackend.figure_format = 'svg'\n", "from IPython.display import Image\n", "url = ('https://raw.githubusercontent.com/ImperialCollegeLondon/sharpy/dev_example/docs/' + \n", " 'source/content/example_notebooks/images/simple_wing_scheme.png')\n", "Image(url=url, width=800)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Geometry\n", "chord = 1. # Chord of the wing\n", "aspect_ratio = 16. # Ratio between lenght and chord: aspect_ratio = length/chord\n", "wake_length = 50 # Length of the wake in chord lengths\n", "\n", "# Discretization\n", "num_node = 21 # Number of nodes in the structural discretisation\n", " # The number of nodes will also define the aerodynamic panels in the \n", " # spanwise direction\n", "num_chord_panels = 4 # Number of aerodynamic panels in the chordwise direction\n", "num_points_camber = 200 # The camber line of the wing will be defined by a series of (x,y) \n", " # coordintes. Here, we define the size of the (x,y) vectors\n", "\n", "# Structural properties of the beam cross section\n", "mass_per_unit_length = 0.75 # Mass per unit length\n", "mass_iner_x = 0.1 # Mass inertia around the local x axis\n", "mass_iner_y = 0.05 # Mass inertia around the local y axis\n", "mass_iner_z = 0.05 # Mass inertia around the local z axis\n", "pos_cg_B = np.zeros((3)) # position of the centre of mass with respect to the elastic axis\n", "EA = 1e7 # Axial stiffness\n", "GAy = 1e6 # Shear stiffness in the local y axis\n", "GAz = 1e6 # Shear stiffness in the local z axis\n", "GJ = 1e4 # Torsional stiffness\n", "EIy = 2e4 # Bending stiffness around the flapwise direction\n", "EIz = 5e6 # Bending stiffness around the edgewise direction\n", "\n", "# Operation\n", "WSP = 2. # Wind speed \n", "air_density = 0.1 # Air density \n", "\n", "# Time discretization\n", "end_time = 5.0 # End time of the simulation\n", "dt = chord/num_chord_panels/WSP # Always keep one timestep per panel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we are going to compute the static equilibrium of the wing at an angle of attack (aoa_ini_deg). Next, we will change the angle of attack to aoa_end_deg and we will compute the dynamic response of the wing." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "aoa_ini_deg = 2. # Angle of attack at the beginning of the simulation \n", "aoa_end_deg = 1. # Angle of attack at the end of the simulation " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For now, we keep the size of the wake panels equal to the distance covered by the flow in one time step. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Structural model\n", "\n", "This section creates a class called AeroelasticInformation. This class is not directly used by SHARPy, it has been thought as an intermediate step between common engineering inputs and the input information that SHARPy needs. \n", "\n", "It has many functionalities (rotate objects, assembly simple strucures together ...) which can be looked up in the [documentation](https://ic-sharpy.readthedocs.io/en/master/includes/utils/generate_cases/index.html).\n", "\n", "Let's initialise an Aeroelastic system that will include a StructuralInformation class and an AerodynamicInformation class:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "wing = gc.AeroelasticInformation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The attibutes that have to be defined (not all of them are compulsory) can be displayed with the code below. They constitute the input parameters to SHARPy and their names are intuitive, however, a complete description of the input variables to SHARPy can be found in the documentation: \n", "[structural inputs](https://ic-sharpy.readthedocs.io/en/master/content/casefiles.html#fem-file) and [aerodynamic inputs](https://ic-sharpy.readthedocs.io/en/master/content/casefiles.html#aerodynamics-file)\n", "\n", "For example, the structural properties are:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "dict_keys(['num_node_elem', 'num_node', 'num_elem', 'coordinates', 'connectivities', 'elem_stiffness', 'stiffness_db', 'elem_mass', 'mass_db', 'frame_of_reference_delta', 'structural_twist', 'boundary_conditions', 'beam_number', 'body_number', 'app_forces', 'lumped_mass_nodes', 'lumped_mass', 'lumped_mass_inertia', 'lumped_mass_position', 'lumped_mass_mat', 'lumped_mass_mat_nodes'])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "wing.StructuralInformation.__dict__.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example, the connectivities between the nodes required by the finite element solver are empty so far:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "None\n" ] } ], "source": [ "print(wing.StructuralInformation.connectivities)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For a list of methods that allow us to modify the structure, Check the [documentation ](https://ic-sharpy.readthedocs.io/en/master/includes/utils/generate_cases/index.html).\n", "\n", "First, we need to define the basic characteristics of the equations as the number of nodes, the number of nodes per element, the number of elements and the location of the nodes in the space:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# Define the number of nodes and the number of nodes per element\n", "wing.StructuralInformation.num_node = num_node\n", "wing.StructuralInformation.num_node_elem = 3\n", "# Compute the number of elements assuming basic connections\n", "wing.StructuralInformation.compute_basic_num_elem()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0. 0. 0. ]\n", " [ 0. 0.8 0. ]\n", " [ 0. 1.6 0. ]\n", " [ 0. 2.4 0. ]\n", " [ 0. 3.2 0. ]\n", " [ 0. 4. 0. ]\n", " [ 0. 4.8 0. ]\n", " [ 0. 5.6 0. ]\n", " [ 0. 6.4 0. ]\n", " [ 0. 7.2 0. ]\n", " [ 0. 8. 0. ]\n", " [ 0. 8.8 0. ]\n", " [ 0. 9.6 0. ]\n", " [ 0. 10.4 0. ]\n", " [ 0. 11.2 0. ]\n", " [ 0. 12. 0. ]\n", " [ 0. 12.8 0. ]\n", " [ 0. 13.6 0. ]\n", " [ 0. 14.4 0. ]\n", " [ 0. 15.2 0. ]\n", " [ 0. 16. 0. ]]\n" ] } ], "source": [ "# Generate an array with the location of the nodes\n", "node_r = np.zeros((num_node, 3))\n", "node_r[:,1] = np.linspace(0, chord*aspect_ratio, num_node)\n", "print(node_r)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following function creates a uniform beam from the previous parameters. On top of assigning the previous parameters, it defines other variables such as the connectivities between nodes" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "wing.StructuralInformation.generate_uniform_beam(node_r,\n", " mass_per_unit_length,\n", " mass_iner_x,\n", " mass_iner_y,\n", " mass_iner_z,\n", " pos_cg_B,\n", " EA,\n", " GAy,\n", " GAz,\n", " GJ,\n", " EIy,\n", " EIz,\n", " num_node_elem = wing.StructuralInformation.num_node_elem,\n", " y_BFoR = 'x_AFoR',\n", " num_lumped_mass=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we now show the connectivities between nodes, the function has created them for us. [Further information](https://ic-sharpy.readthedocs.io/en/master/content/casefiles.html?highlight=connectivities#fem-file)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0 2 1]\n", " [ 2 4 3]\n", " [ 4 6 5]\n", " [ 6 8 7]\n", " [ 8 10 9]\n", " [10 12 11]\n", " [12 14 13]\n", " [14 16 15]\n", " [16 18 17]\n", " [18 20 19]]\n" ] } ], "source": [ "print(wing.StructuralInformation.connectivities)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's define the boundary conditions as clamped for node 0 and free for the last node (-1) [Further information](https://ic-sharpy.readthedocs.io/en/master/content/casefiles.html?highlight=boundary#fem-file)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "wing.StructuralInformation.boundary_conditions = 1\n", "wing.StructuralInformation.boundary_conditions[-1] = -1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Aerodynamics\n", "We need to define the number of panels in the wake (wake_panels) and the camber line of the wing which, in this case, is flat. [Further information](https://ic-sharpy.readthedocs.io/en/master/content/casefiles.html#aerodynamics-file)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "# Compute the number of panels in the wake (streamwise direction) based on the previous paramete\n", "wake_panels = int(wake_length*chord/dt)\n", "\n", "# Define the coordinates of the camber line of the wing\n", "wing_camber = np.zeros((1, num_points_camber, 2))\n", "wing_camber[0, :, 0] = np.linspace(0, 1, num_points_camber)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following function creates an aerodynamic surface uniform on top of the beam that we have already created" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# Generate blade aerodynamics\n", "wing.AerodynamicInformation.create_one_uniform_aerodynamics(wing.StructuralInformation,\n", " chord = chord,\n", " twist = 0.,\n", " sweep = 0.,\n", " num_chord_panels = num_chord_panels,\n", " m_distribution = 'uniform',\n", " elastic_axis = 0.5,\n", " num_points_camber = num_points_camber,\n", " airfoil = wing_camber)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Summary\n", "Now, we have all the inputs that we need for SHARPy. In this section, a first simulation with SHARPy is run. However, it will perform no computation, it will just load the data so we can plot the system we have just created.\n", "\n", "SHARPy runs a series of [solvers](https://ic-sharpy.readthedocs.io/en/master/content/solvers.html) and [postprocessors](https://ic-sharpy.readthedocs.io/en/master/content/postproc.html) in the order indicated by the flow variable in the SHARPy dictionary.\n", "\n", "The generate cases module also allows us to show all the avaiable solvers and all the parameters they accept as inputs. To do so, we need to run the following commands:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "_BaseStructural\n", "AerogridLoader\n", "BeamLoader\n", "DynamicCoupled\n", "DynamicUVLM\n", "LinDynamicSim\n", "LinearAssembler\n", "Modal\n", "NoAero\n", "NonLinearDynamic\n", "NonLinearDynamicCoupledStep\n", "NonLinearDynamicMultibody\n", "NonLinearDynamicPrescribedStep\n", "NonLinearStatic\n", "PrescribedUvlm\n", "RigidDynamicCoupledStep\n", "RigidDynamicPrescribedStep\n", "StaticCoupled\n", "StaticCoupledRBM\n", "StaticTrim\n", "StaticUvlm\n", "StepLinearUVLM\n", "StepUvlm\n", "_BaseTimeIntegrator\n", "NewmarkBeta\n", "GeneralisedAlpha\n", "Trim\n", "LiftDistribution\n", "Cleanup\n", "AerogridPlot\n", "BeamPlot\n", "UDPout\n", "PlotFlowField\n", "PickleData\n", "SaveParametricCase\n", "FrequencyResponse\n", "StallCheck\n", "WriteVariablesTime\n", "CreateSnapshot\n", "SHARPy\n", "SaveData\n", "AeroForcesCalculator\n", "AsymptoticStability\n", "BeamLoads\n", "FloatingForces\n", "ShearVelocityField\n", "DynamicControlSurface\n", "GustVelocityField\n", "TrajectoryGenerator\n", "TurbVelocityField\n", "SteadyVelocityField\n", "ModifyStructure\n", "StraightWake\n", "HelicoidalWake\n", "GridBox\n", "BumpVelocityField\n", "TurbVelocityFieldBts\n" ] } ], "source": [ "# Gather data about available solvers\n", "SimInfo = gc.SimulationInformation() # Initialises the SimulationInformation class\n", "SimInfo.set_default_values() # Assigns the default values to all the solvers\n", "\n", "# Print the available solvers and postprocessors\n", "for key in SimInfo.solvers.keys():\n", " print(key)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's output as an example, the input parameters of the BeamLoader solver. This solver is in charge of loading the structural infromation" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'unsteady': True, 'orientation': [1.0, 0, 0, 0], 'for_pos': [0.0, 0, 0]}" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "SimInfo.solvers['BeamLoader']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following dictionary defines the basic inputs of SHARPy including the solvers to be run (flow variable), the case name and the route in the file system. In this case we are turning off the screen output because it is useless because we are not running any computation. " ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "SimInfo.solvers['SHARPy']['flow'] = ['BeamLoader',\n", " 'AerogridLoader']\n", "\n", "SimInfo.solvers['SHARPy']['case'] = 'plot'\n", "SimInfo.solvers['SHARPy']['route'] = './'\n", "SimInfo.solvers['SHARPy']['write_screen'] = 'off'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We do not modify any of the default input paramteters in BeamLoader but we need to define the initial wake shape that we want and its size (wake_panels)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "SimInfo.solvers['AerogridLoader']['unsteady'] = 'on'\n", "SimInfo.solvers['AerogridLoader']['mstar'] = wake_panels\n", "SimInfo.solvers['AerogridLoader']['freestream_dir'] = np.array([1.,0.,0.])\n", "SimInfo.solvers['AerogridLoader']['wake_shape_generator'] = 'StraightWake'\n", "SimInfo.solvers['AerogridLoader']['wake_shape_generator_input'] = {'u_inf': WSP,\n", " 'u_inf_direction' : np.array(\n", " [np.cos(aoa_ini_deg*deg2rad),\n", " 0.,\n", " np.sin(aoa_ini_deg*deg2rad)]),\n", " 'dt': dt}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following functions write the input files needed by SHARPy" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "gc.clean_test_files(SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])\n", "wing.generate_h5_files(SimInfo.solvers['SHARPy']['route'], SimInfo.solvers['SHARPy']['case'])\n", "SimInfo.generate_solver_file()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following line of code runs SHARPy inside the jupyter notebook. It is equivalent to run in a terminal: sharpy plot.sharpy being \"plot\" the case name defined above" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "sharpy_output = sharpy.sharpy_main.main(['',\n", " SimInfo.solvers['SHARPy']['route'] +\n", " SimInfo.solvers['SHARPy']['case'] +\n", " '.sharpy'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot the case we have created. The function below provides a quick way of generating an interactive plot but it is not suitable for large cases. Luckily SHARPy has other more complex plotting capabilities.\n", "\n", "If you download the original [jupyter notebook](https://ic-sharpy.readthedocs.io/en/master/content/examples.html), change thee geometric parameters above and rerun the jupyter notebook you will see its impact on the geometry.\n", "\n", "Only the first 6 wake panels are plotted for efficiency." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/html": [ "