Python_validation_codes/nastran_panair_complete.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nastran_Panair_Complete tutorial\n",
    "\n",
    "This function solves a static aeroelastic problem using a partitioned approach Nastran95 and Panair for the NASA CRM wing. \n",
    "\n",
    "The method is based on the OpenMDAO package *aerostructures* and allows to change some wing parameters interpolating the nodes in the deformed geometry.\n",
    "\n",
    "The input files are the nastran .inp template *nastran_static_template.inp* and the list of nodes in the structural mesh *nastran_input_geometry.inp*.\n",
    "The other files required are related to the airfoil shape: *crm.eta65.unswept31.5deg.sharp.te.txt*, *aero_jig.msh* and *aero_jig.geo*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Main\n",
    "\n",
    "Import the *aerostructures* functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from openmdao.api import Problem, Group, IndepVarComp, ScipyGMRES, SqliteRecorder, ExecComp\n",
    "\n",
    "from aerostructures import NastranStatic, DisplacementTransfer, Panair, LoadTransfer, Interpolation, PanairMesher, StructureMesher, PlanformGeometry, StaticStructureProblemDimensions, AeroProblemDimensions, StaticStructureProblemParams, AeroProblemParams, NLGaussSeidel, GeometricalProperties\n",
    "\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interpolation setting, airfoil shape and cruise conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "if __name__ == \"__main__\":\n",
    "\n",
    "    #Interpolation function type and setup\n",
    "    function_type = 'linear'\n",
    "    bias_morph = (1.,1.,1.)\n",
    "    bias_inter = (1.,100.,1.)\n",
    "\n",
    "    #Symmetry plane index\n",
    "    sym_plane_index = 2\n",
    "    case_name = 'alpha_c'\n",
    "    #Number of wing sections\n",
    "    n_sec = 8\n",
    "\n",
    "    #Position (index) of the wing break\n",
    "    b_sec = 4\n",
    "\n",
    "    #Airfoil file\n",
    "    ref_airfoil_file = 'crm.eta65.unswept31.5deg.sharp.te.txt'\n",
    "\n",
    "    #Problem parameters\n",
    "    #Speed of sound\n",
    "    a = 297.4\n",
    "\n",
    "    Sw = 383.689555\n",
    "    V = 252.16168\n",
    "    Mach = V/a\n",
    "    rho_a = 0.38058496\n",
    "    # h = 10500.\n",
    "    alpha_0 = 1.340"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More parameters, b_0 and sweep_0 have been changed for example during the fidelity correlation analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "    b_0 = 58.7629\n",
    "    sweep_0 = 37.16\n",
    "    \n",
    "    b_baseline = 58.7629\n",
    "    c = 7.00532\n",
    "    E = 6.89e10\n",
    "    nu = 0.31\n",
    "    rho_s = 2795.67\n",
    "    #Reference aircraft weight (mass units)\n",
    "    W_ref = 226796.185\n",
    "    #Wing weight (full span) of the reference aircraft (mass units)\n",
    "    W_ref_wing = 26400.\n",
    "    #Airframe weight (complete aircraft excluding wing structure, mass units)\n",
    "    W_airframe = W_ref - W_ref_wing\n",
    "    #Yield stress (can also be used as ultimate stress if FS = 1.5)\n",
    "    sigma_y = 450.e6\n",
    "    #Factor of safety\n",
    "    FS = 1.\n",
    "    #Cruise load factor\n",
    "    n = 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wing planform parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "    #Sectional properties (that are not design variables)\n",
    "    y_le_baseline = np.array([0., 2.938145, 7.3453752, 10.8711746, 16.1598356, 20.5670658, 24.974296, 29.3815262])\n",
    "    z_le = np.array([4.424397971, 4.44511389, 4.476187859, 4.501047142, 4.538335797, 4.569409766, 4.600483735, 4.631557704])\n",
    "    c_0 = np.array([13.6189974, 11.9001794, 9.3216984,\n",
    "                    7.2588628, 5.9643264, 4.8855376, 3.8067488, 2.72796])\n",
    "    tc_0 = np.array([0.1542, 0.138, 0.1137, 0.1052, 0.0988, 0.0962, 0.0953, 0.095])\n",
    "    th = tc_0*c_0\n",
    "    camc = np.array([0.0003, 0.0012, 0.0037, 0.0095, 0.0146, 0.0158, 0.0161, 0.0009])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NASA CRM wing model values used to redefine the elements properties (stringers and spars in particular) for low fidelity meshes. GeometricalProperties takes this values as input and reassign properties in order to maintain the overall stiffness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "    #Parameters to define the rod section\n",
    "    t_0 = np.array([.00635, .005334, .004572, .003302, .00254,\n",
    "                  .001651, .01905, .01524, .0127, .009525, .00508, .00254]) \n",
    "    l_0 = np.array([.358097605120250,   .326274684415706,   .271148349471907,\n",
    "                  .221233096812120,   .179653543701401,   .158489131881446])\n",
    "\n",
    "    #Parameters to define the bar properties\n",
    "    s_0_tot = np.array([.0035992, .0074866, .0003097, .0005161, .024981, \n",
    "                    .0001064, .0134186, .0027664 ])\n",
    "    Ix_0_tot = np.array([1.68542e-07, 7.293e-07, 2.745e-08, 5.695e-07, \n",
    "                         8.1202e-06, 2.421e-08, 1.4807e-05, 6.2946e-07 ])\n",
    "    Iy_0_tot = np.array([1.25906e-06, 5.4472e-06, 2.051e-07, 7.624e-08, \n",
    "                         6.0654e-05, 3.24e-09, 1.98224e-06, 8.424e-08])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reading the model properties from the template file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "    structure_problem_dimensions = StaticStructureProblemDimensions()\n",
    "    aero_problem_dimensions = AeroProblemDimensions()\n",
    "\n",
    "    ns = structure_problem_dimensions.ns\n",
    "    ns_all = structure_problem_dimensions.ns_all\n",
    "    node_id = structure_problem_dimensions.node_id\n",
    "    node_id_all = structure_problem_dimensions.node_id_all\n",
    "    n_stress = structure_problem_dimensions.n_stress\n",
    "    \n",
    "    tn = structure_problem_dimensions.tn\n",
    "    mn = structure_problem_dimensions.mn\n",
    "    sn = structure_problem_dimensions.sn\n",
    "    an = structure_problem_dimensions.an\n",
    "    Ln = structure_problem_dimensions.Ln\n",
    "    Vn = structure_problem_dimensions.Vn\n",
    "\n",
    "    \n",
    "    structure_problem_params = StaticStructureProblemParams(node_id, node_id_all)\n",
    "    aero_problem_params = AeroProblemParams()\n",
    "\n",
    "    aero_problem_dimensions = AeroProblemDimensions()\n",
    "\n",
    "    na = aero_problem_dimensions.na\n",
    "    na_unique = aero_problem_dimensions.na_unique\n",
    "    network_info = aero_problem_dimensions.network_info\n",
    "\n",
    "    aero_problem_params = AeroProblemParams()\n",
    "\n",
    "    geom_propts = GeometricalProperties(Ln,Vn,tn,an,sn,t_0,l_0,s_0_tot,Ix_0_tot,Iy_0_tot)\n",
    "    \n",
    "    t = geom_propts.t\n",
    "    a = geom_propts.a\n",
    "    s = geom_propts.s\n",
    "    Ix = geom_propts.Ix\n",
    "    Iy = geom_propts.Iy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last wing parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "    theta = np.array([6.691738003, 4.545042708, 2.793550837, 1.673916686,\n",
    "                      0.754303126, 0.91369482, 1.136056807, 0.272576679])\n",
    "\n",
    "    cr_0 = 13.6189974\n",
    "\n",
    "    cb_0 = 7.2588628\n",
    "\n",
    "    ct_0 = 2.72796\n",
    "\n",
    "    \n",
    "    \n",
    "    #X-position of the leading edge at the root\n",
    "    xr = 22.9690676\n",
    "    \n",
    "    #Coordinates of aerodynamic and structure matching meshes\n",
    "    xa_b = aero_problem_params.apoints_coord_unique\n",
    "    xs_b = structure_problem_params.node_coord_all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Independent variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<openmdao.components.indep_var_comp.IndepVarComp at 0x9273588>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "    top = Problem()\n",
    "    top.root = root = Group()\n",
    "\n",
    "    #Add independent variables (parameters)\n",
    "    root.add('wing_area', IndepVarComp('Sw', Sw), promotes=['*'])\n",
    "    root.add('Airspeed', IndepVarComp('V', V), promotes=['*'])\n",
    "    root.add('air_density', IndepVarComp('rho_a', rho_a), promotes=['*'])\n",
    "    root.add('Mach_number', IndepVarComp('Mach', Mach), promotes=['*'])\n",
    "    root.add('baseline_wing_span', IndepVarComp('b_baseline', b_baseline), promotes=['*'])\n",
    "    root.add('wing_chord', IndepVarComp('c', c), promotes=['*'])\n",
    "    root.add('Youngs_modulus', IndepVarComp('E', E), promotes=['*'])\n",
    "    root.add('Poissons_ratio', IndepVarComp('nu', nu), promotes=['*'])\n",
    "    root.add('material_density', IndepVarComp('rho_s', rho_s), promotes=['*'])\n",
    "    root.add('airframe_mass', IndepVarComp('W_airframe', W_airframe), promotes=['*'])\n",
    "    root.add('Tensile_Yield_Strength', IndepVarComp('sigma_y', sigma_y), promotes=['*'])\n",
    "    root.add('factor_safety', IndepVarComp('FS', FS), promotes=['*'])\n",
    "    root.add('y_leading_edge_baseline', IndepVarComp('y_le_baseline', y_le_baseline), promotes=['*'])\n",
    "    root.add('z_leading_edge', IndepVarComp('z_le', z_le), promotes=['*'])\n",
    "    root.add('airfoil_thickness', IndepVarComp('th', th), promotes=['*'])\n",
    "    root.add('camber_chord_ratio', IndepVarComp('camc', camc), promotes=['*'])\n",
    "    root.add('base_aerodynamic_mesh', IndepVarComp('xa_b', xa_b), promotes=['*'])\n",
    "    root.add('base_structure_mesh', IndepVarComp('xs_b', xs_b), promotes=['*'])\n",
    "    root.add('cruise_load_factor', IndepVarComp('n', n), promotes=['*'])\n",
    "    root.add('root_leading_edge_x', IndepVarComp('xr', xr), promotes=['*'])\n",
    "    root.add('wing_twist', IndepVarComp('theta', theta), promotes=['*'])\n",
    "    root.add('root_chord', IndepVarComp('cr', cr_0), promotes=['*'])\n",
    "    root.add('break_chord', IndepVarComp('cb', cb_0), promotes=['*'])\n",
    "    root.add('tip_chord', IndepVarComp('ct', ct_0), promotes=['*'])\n",
    "    root.add('sweep_angle', IndepVarComp('sweep', sweep_0), promotes=['*'])\n",
    "    root.add('wing_span', IndepVarComp('b', b_0), promotes=['*'])\n",
    "    root.add('angle_of_attack', IndepVarComp('alpha', alpha_0), promotes=['*'])\n",
    "\n",
    "    \n",
    "    root.add('thicknesses', IndepVarComp('t', t), promotes=['*'])\n",
    "    if an > 0:\n",
    "        root.add('rod_sections', IndepVarComp('a', a), promotes=['*'])   \n",
    "    root.add('bar_sections', IndepVarComp('s', s), promotes=['*'])\n",
    "    root.add('first_moment_inertia', IndepVarComp('Ix', Ix), promotes=['*'])\n",
    "    root.add('second_moment_inertia', IndepVarComp('Iy', Iy), promotes=['*'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interpolation and remeshing components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<openmdao.components.exec_comp.ExecComp at 0x8e267f0>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "    #Interpolation Components\n",
    "    root.add('interp_struct_morph', Interpolation(ns_all, na_unique, function = function_type, bias = bias_morph))\n",
    "    root.add('interp_mda', Interpolation(na, ns, function = function_type, bias = bias_inter), promotes=['*'])\n",
    "\n",
    "    #Geometry and meshing Components\n",
    "    root.add('planform_geometry', PlanformGeometry(n_sec, b_sec), promotes=['*'])\n",
    "    root.add('aerodynamic_mesher', PanairMesher(n_sec, na, na_unique, network_info, ref_airfoil_file), promotes=['*'])\n",
    "    root.add('structure_mesher', StructureMesher(na_unique, node_id, node_id_all), promotes=['*'])\n",
    "\n",
    "    root.add('y_leading_edge', ExecComp(\n",
    "        'y_le = b/b_baseline*y_le_baseline', y_le=np.zeros(len(y_le_baseline), dtype=float), y_le_baseline=np.zeros(len(y_le_baseline), dtype=float)), promotes=['*'])\n",
    "\n",
    "    root.add('tc_ratio', ExecComp(\n",
    "        'tc = th/chords', tc=np.zeros(n_sec, dtype=float), th=np.zeros(n_sec, dtype=float), chords=np.zeros(n_sec, dtype=float)), promotes=['*'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explicit connections and group solvers definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "    #Explicit connections\n",
    "    root.connect('interp_struct_morph.H', 'G')\n",
    "    root.connect('xa_b', 'interp_struct_morph.node_coord')\n",
    "    root.connect('xs_b', 'interp_struct_morph.apoints_coord')\n",
    "\n",
    "    mda = Group()\n",
    "\n",
    "    #Add disciplines to the group\n",
    "    mda.add('displacement_transfer', DisplacementTransfer(na, ns), promotes=['*'])\n",
    "    mda.add('aerodynamics', Panair(na, network_info, case_name, sym_plane_index), promotes=['*'])\n",
    "    mda.add('load_transfer', LoadTransfer(na, ns), promotes=['*'])\n",
    "    mda.add('structures', NastranStatic(node_id, node_id_all, n_stress, tn, mn, sn, an, case_name), promotes=['*'])\n",
    "\n",
    "    #Define solver type and tolerance for MDA\n",
    "    mda.nl_solver = NLGaussSeidel()\n",
    "    # mda.nl_solver.options['maxiter'] = 0\n",
    "    mda.nl_solver.options['rutol'] = 1.e-2 #2 origin\n",
    "    mda.nl_solver.options['use_aitken'] = True\n",
    "    mda.nl_solver.options['aitken_alpha_min'] = 0.1\n",
    "    mda.nl_solver.options['aitken_alpha_max'] = 1.5\n",
    "\n",
    "    mda.ln_solver = ScipyGMRES()\n",
    "\n",
    "    root.add('mda_group', mda, promotes=['*'])\n",
    "\n",
    "    #Recorder\n",
    "    recorder = SqliteRecorder('mda.sqlite3')\n",
    "    recorder.options['record_metadata'] = False\n",
    "    # recorder.options['includes'] =\n",
    "    top.root.mda_group.nl_solver.add_recorder(recorder)\n",
    "\n",
    "    #Define solver type\n",
    "    root.ln_solver = ScipyGMRES()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Problem run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "##############################################\n",
      "Setup: Checking root problem for potential issues...\n",
      "\n",
      "The following parameters have no associated unknowns:\n",
      "a\n",
      "m\n",
      "\n",
      "The following components have no connections:\n",
      "Tensile_Yield_Strength\n",
      "airframe_mass\n",
      "factor_safety\n",
      "Group 'mda_group' has the following cycles: [['aerodynamics', 'load_transfer', 'structures', 'displacement_transfer']]\n",
      "\n",
      "The following params are connected to unknowns that are updated out of order, so their initial values may contain uninitialized unknown values: ['mda_group.aerodynamics.delta']\n",
      "mda_group.structures:\n",
      "The following input files are missing at setup  time: ['nastran_static_alpha_c.inp']\n",
      "\n",
      "\n",
      "Setup: Check of root problem complete.\n",
      "##############################################\n",
      "\n"
     ]
    }
   ],
   "source": [
    "    start1 = time.time()\n",
    "    top.setup()\n",
    "    end1 = time.time()\n",
    "    \n",
    "    start2 = time.time()\n",
    "    top.run()\n",
    "    end2 = time.time()\n",
    "    \n",
    "    print(\"Set up time = \" + str(end1 - start1))\n",
    "    print(\"Run time = \" + str(end2 - start2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Post Processing\n",
    "\n",
    "Once solved the problem, one can access to the post processing results for the structural deformation, stresses and pressure distribution.\n",
    "\n",
    "*Nastran_static_postproc.py* produce the .msh file containing the structure deformed and the elements Von Mises stress. At the bottom of the file one can find the maximum value of the stress, its position and the tip vertical displacement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run Nastran_static_postproc.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The postprocessing of the Panair output gives the pressure distribution. the drag and lift coefficient and deformed wing shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run wgs_to_gmsh.py\n",
    "%run panair_post.py"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}