AeroComBAT Tutorials

The following are tutorials intended convey how best to use AeroComBAT.

Tutorial 1 - Material Library and Classical Lamination Theory

This tutorial is intended to expose the user to the material library class as well as creating CLT laminates with the laminate class. Note that these laminate objects are integral to meshing the cross-section (XSect) objects.

# =============================================================================
# AEROCOMBAT TUTORIAL 1 - MATERIAL lIBRARY AND CLASSICAL LAMINATION THEORY
# =============================================================================

# IMPORT SYSTEM PACKAGES
# ======================
import sys
import os
# Append the root to the system path
sys.path.append(os.path.abspath('..'))

# IMPORT AEROCOMBAT CLASSES
# =========================
from AeroComBAT.Structures import MaterialLib, Laminate
from AeroComBAT.Utilities import RotationHelper
from AeroComBAT.tabulate import tabulate
import numpy as np


# MATERIAL lIBRARY VALIDATION
# ===========================
# Generate Empty Material Library
matlib = MaterialLib()
# Add a graphite orthotropic material
matlib.addMat(1, 'Graphite-Polymer Composite ortho', 'ortho', \
                [155.0e9, 12.1e9, 12.1e9, .458, .248, .248, 3.2e9, 4.4e9,\
                4.4e9, 1.7e3], .15e-3)
# Add a graphite transversely isotropic material
matlib.addMat(2, 'Graphite-Polymer Composite', 'trans_iso', \
                [155.0e9, 12.1e9, .458, .248, 4.4e9, 1.7e3], .15e-3)
# Add a glass transversely isotropic material
matlib.addMat(3, 'Glass-Polymer Composite', 'trans_iso', \
                [50.0e9, 15.2e9, .428, .254, 4.7e9, 1.2e3], .15e-3)
# Add a T300 transversely isotropic material
matlib.addMat(4, 'T300/5208', 'trans_iso', \
                [181.0e9, 10.3e9, .458, .28, 7.17e9, 1.8e3], .15e-3)
# Add a aluminum isotropic material
matlib.addMat(5, 'AL-2050', 'iso',[75.8, 0.33, 2.7e3], .15e-3)
# Add a rotated T300 transversely isotropic material
matlib.addMat(6, 'T300/5208', 'trans_iso', \
                [181.0e9, 10.3e9, .458, .28, 7.17e9, 1.8e3], .15e-3,th = [0.,45.,0.])
# Print a summary of the mat
matlib.printSummary()
# Get the material associated with MID 1
mat1 = matlib.getMat(1)
# Get the compliance matrix of the material mat1
Smat1 = mat1.Smat
# Get the stiffness matrix of the material mat1 and round for accuracy
Cmat1 = np.around(mat1.Cmat/1e9,decimals=2)
Cmat1[0,0] = np.around(Cmat1[0,0])
# The following matricies are the correct compliance and stiffness matricies
Smat1Test = np.array([[6.45e-12,-1.6e-12,-1.6e-12,0.,0.,0.],\
                      [-1.6e-12,82.6e-12,-37.9e-12,0.,0.,0.],\
                      [-1.6e-12,-37.9e-12,82.6e-12,0.,0.,0.],\
                      [0.,0.,0.,312e-12,0.,0.],\
                      [0.,0.,0.,0.,227e-12,0.],\
                      [0.,0.,0.,0.,0.,227e-12]])
Cmat1Test = np.array([[158e9,5.64e9,5.64e9,0.,0.,0.],\
                      [5.64e9,15.51e9,7.21e9,0.,0.,0.],\
                      [5.64e9,7.21e9,15.51e9,0.,0.,0.],\
                      [0.,0.,0.,3.2e9,0.,0.],\
                      [0.,0.,0.,0.,4.4e9,0.],\
                      [0.,0.,0.,0.,0.,4.4e9]])/1e9
# Check to make sure the calculated values are correct
np.testing.assert_array_almost_equal(Smat1,Smat1Test,decimal=12)
np.testing.assert_array_almost_equal(Cmat1,Cmat1Test,decimal=12)


# MATERIAL PROPERTY ROTATION HELPER VALIDATION
# ============================================
# Create a rotation helper object
rh = RotationHelper()
# Create an array of x-y-z rotations
th = [0.,45.,0.]
# Initialize a stiffness matrix
C = np.array([[1.8403e11,5.4101e9,5.4101e9,0.,0.,0.],\
              [5.4101e9,1.31931e10,6.12866e9,0.,0.,0.],\
              [5.4101e9,6.12866e9,1.31931e10,0.,0.,0.],\
              [0.,0.,0.,5.21455e9,0.,0.],\
              [0.,0.,0.,0.,7.17e9,0.],\
              [0.,0.,0.,0.,0.,7.17e9]])
# Convert it into a compliance matrix
S = np.linalg.inv(C)
# Rotate the compliance matrix
Sp = rh.transformCompl(S,th,xsect=True)
# Convert it back to a stiffness matrix
Cp = np.linalg.inv(Sp)
print('The rotated stiffness matrix:')
print(tabulate(np.around(Cp-C,decimals=3),tablefmt="fancy_grid"))

# =============================================================================
# CLT VALIDATION
# =============================================================================
# Initialize the number of plies per each orientation
n_i = [1,1,1,1]
# Initialize the materials to be used at each orientation
m_i = [4,4,4,4]
# Initialize the angle orientations for the plies
th = [30,-30,0,45]
# Create a laminate with default orientations (for 4 orientations, this will
# default to th_defalt = [0,45,90,-45])
lam1 = Laminate(n_i,m_i,matlib)
# Print a summary of laminate 1
print('Laminate 1 summary:')
lam1.printSummary(decimals=3)
# Create a laminate with default orientations (for more or less than 4
# orientations, th_default = [0]*len(n_i))
lam2 = Laminate(n_i+n_i,m_i+m_i,matlib)
# Print summary of laminate 2
print('Laminate 2 summary:')
lam2.printSummary(decimals=3)
# Create a laminate using the above rotation orientations
lam3 = Laminate(n_i,m_i,matlib,th=th)
# Print Summary of laminate 3
print('Laminate 3 summary:')
lam3.printSummary(decimals=3)

Tutorial 2 - CQUADX and Airfoil Classes

This tutorial is intended to expose the user to the CQUADX and airfoil classes.

# =============================================================================
# AEROCOMBAT TUTORIAL 2 - CQUADX AND AIRFOIL
# =============================================================================

# IMPORT SYSTEM PACKAGES
# ======================
import sys
import os

sys.path.append(os.path.abspath('..'))

# IMPORT AEROCOMBAT CLASSES
# =========================
from AeroComBAT.Structures import Node, MaterialLib, CQUADX
from AeroComBAT.Aerodynamics import Airfoil

# IMPORT NUMPY MODULES
# ====================
import numpy as np
import matplotlib.pyplot as plt

# Material Info
mat_lib = MaterialLib()
# Add an aluminum isotropic material
mat_lib.addMat(1, 'AL-2050', 'iso',[75.8, 0.33, 2.7e3], .15e-3)
                

# CQUADX 2D ELEMENT CREATION
# ==========================
# Create a node 1 object
n1 = Node(1,[0.,0.,0.])
# Create a node 2 object
n2 = Node(2,[2.,0.,0.])
# Create a node 3 object
n3 = Node(3,[2.,3.,0.])
# Create a node 4 object
n4 = Node(4,[0.,5.,0.])
# Create a CQUADX element
elem1 = CQUADX(1,[n1,n2,n3,n4],1,mat_lib)
# Print a summary of the element
elem1.printSummary(nodes=True)

# AIRFOIL OUTER MOLD LINE VALIDATION
# ==================================
# Initialize a chord length of 1
c = 1.
# Create an airfoil object with a 'box' profile
af1 = Airfoil(c,name='box')
# Generate a set of non-dimensional x-coordinates
x = np.linspace(-.5,.5,50)
# Create the upper and lower box airfoil curves
xu,yu,xl,yl = af1.points(x)
# Create a matplotlib figure
plt.figure(num=1)
plt.plot(xu,yu)
plt.hold(True)
plt.plot(xl,yl)
plt.axes().set_aspect('equal', 'datalim')
plt.xlabel('x coordinate along the airfoil')
plt.ylabel('y coordinate along the airfoil')
plt.title('Box airfoil profile')
plt.hold(False)

# Create a NACA2412 airfoil profile
af2 = Airfoil(c,name='NACA2412')
# Generate a set of non-dimensional x-coordinates
x = np.linspace(0,1.,500)
# Create the upper and lower airfoil curves
xu,yu,xl,yl = af2.points(x)
# Create a matplotlib figure
plt.figure(num=2)
plt.plot(xu,yu)
plt.hold(True)
plt.plot(xl,yl)
plt.hold(False)
plt.axes().set_aspect('equal', 'datalim')

Tutorial 3 - Cross-Section Meshing and Analysis

This is the first extensive tutorial. It exposes the user to meshing several different cross-section types, all with varying complexities.

# =============================================================================
# AEROCOMBAT TUTORIAL 3 - Using XSect Objects
# =============================================================================

# IMPORT SYSTEM PACKAGES
# ======================
import sys
import os

sys.path.append(os.path.abspath('..'))

# IMPORT AEROCOMBAT CLASSES
# =========================
from AeroComBAT.Structures import MaterialLib, Laminate, XSect
from AeroComBAT.Aerodynamics import Airfoil

# IMPORT NUMPY MODULES
# ====================
import numpy as np


# ADD MATERIALS TO THE MATERIAL LIBRARY
# =====================================
# Create a material library object
matLib = MaterialLib()
# Add material property from Hodges 1999 Asymptotically correct anisotropic
# beam theory (Imperial)
matLib.addMat(1,'AS43501-6','trans_iso',[20.6e6,1.42e6,.34,.3,.87e6,.1],.005)
# Add material property from Hodges 1999 Asymptotically correct anisotropic
# beam theory (Imperial)
matLib.addMat(2,'AS43501-6*','trans_iso',[20.6e6,1.42e6,.34,.42,.87e6,.1],.005)
# Add an aluminum material (SI)
matLib.addMat(3,'AL','iso',[71.7e9,.33,2810],.005)

# CREATE A LAMINATE CROSS-SECTION
# ===============================
# Create a box airfoil object. Note that when creating an airfoil object, only
# the chord length is used. As such, it doesn't truly matter if the airfoil
# has an airfoil profile or a box profile. In this case we will just give it a
# box profile.
# Initialize the chord length
c1 = 1.
# Initialize the non-dimensional starting and stopping points of the cross-
# section. These bounds when dimensionalized will determine the overall
# dimesions of the cross-section. Therefore the total width of the laminate is:
# xdim[1]*c1-xdim[0]*x. In this case, the total width is 2!
xdim1 = [-1.,1.]
af1 = Airfoil(c1,name='box')
# Create a layup schedule for the laminate. In this case, we will select a
# layup schedule of [0_2/45/90/3]_s
th_1 = [0,45,90]
n_1 = [2,1,3]
m_1 = [1,1,1]
# Notice how the orientations are stored in the 'th_1' array, the subscripts are
# stored in the 'n_1' array, and the material information is held in 'm_1'.
# Create the laminate object:
lam1 = Laminate(n_1,m_1,matLib,th=th_1,sym=True)
# In order to make a cross-section, we must add all of the laminates to be used
# to an array:
laminates1 = [lam1]
# We now have all the information necessary to make a laminate beam cross-
# section:
xsect1 = XSect(1,af1,xdim1,laminates1,matLib,typeXSect='laminate',meshSize=2)
# With the cross-section object initialized, let's run the cross-sectional
# analysis to get cross-section stiffnesses, etc.
xsect1.xSectionAnalysis()
# Let's see what our rigid cross-section looks like when plotted in 3D:
xsect1.plotRigid(mesh=True)
# Note that while it might look like the cross-section is made of triangular
# elements, it's actually made of quadrilaterals. This is an artifact of how
# the visualizer mayavi works. Let's get a summary of the cross-section's
# stiffnesses, ect.
xsect1.printSummary(stiffMat=True)
# Notice that from the command line output, all of the important cross-
# sectional geometric properties are located at the origin. By observing the
# cross-section stiffness matrix, it can be seen by the 1,3 entry that there
# is shear-axial coupling. From the non-zero 4,6 entry, we can also tell that
# the cross-section has bending-torsion coupling.
# We can apply a force to the face of this cross-section at the reference axis
# (which in this case is at x,y = 0,0) and see what the stresses look like. In
# this case we'll apply [Fx,Fy,Fz,Mx,My,Mz]=[0.,0.,0.,100.,0.,0.] as if the
# beam belonging to this cross-section were in pure bending.
force1 = np.array([0.,100.,0.,10.,0.,0.])
xsect1.calcWarpEffects(force=force1)
# Having applied the force, let's see what the sigma_zz (normal stresses of the
# beam) look like
xsect1.plotWarped(figName='Laminate Sigma_33 Stress',warpScale=10,\
    contour='sig_33',colorbar=True)
# Let's look at the sigma_13 stress state now since we know there is torsion
# coupling:
xsect1.plotWarped(figName='Laminate Sigma_13 Stress',warpScale=10,\
    contour='sig_13',colorbar=True)
# Notice the increased stress in two of the plies? Recall which ones those are?
# Those plies are the 45 degree plies which are currently taking the shear!

# CREATE A LAMIANTE CROSS-SECTION WITH A DIFFERENT REFERENCE AXIS
# ===============================================================
# The cross-section we just made happened to have all of it's geometrical
# locations (mass center, shear center, tension center) at the origin, which
# is where we applied our force resultant. Suppose we wanted to give the cross-
# section a different reference axis. We can do this by executing the
# xSectionAnalysis method again:
ref_ax = [.5,0.]
# This will move the reference axis (location where we apply forces and
# moments) to x=0.5, y=0. Carrying out the cross-sectional analysis again, we
# get:
xsect1.xSectionAnalysis(ref_ax=ref_ax)
# Let's see how the cross-section's stiffness matrix has changed:
xsect1.printSummary(stiffMat=True)

# Let's not apply the same force resultant and see what the stresses look like:
xsect1.calcWarpEffects(force=force1)
xsect1.plotWarped(figName='Laminate Sigma_33 Stress New Reference Axis',\
    warpScale=10,contour='sig_33',colorbar=True)
xsect1.plotWarped(figName='Laminate Sigma_13 Stress New Reference Axis',\
    warpScale=10,contour='sig_13',colorbar=True)
# Notice how the stress resultants are fairly different once we moved the
# reference axis.

# CREATE A WRAPPED RECTANGULAR BOX-BEAM CROSS-SECTION
# ===================================================
# Layup 1 Box beam (0.5 x 0.923 in^2 box with laminate schedule [0]_6)

# Let's make a slightly more complex cross-section. Using the 'rectBox' key
# word and four laminates, we can make a cross-section box-beam. In the next
# case, we are going to make a cross-section from Hodges 1999 Asymptotically
# correct anisotropic beam theory paper. We will do the most simple box-beam
# in this case, which is the "Layup 1" case:
# Establish the chord length
c2 = 0.53
# Establish the non-dimensional starting and stopping points of the cross-
# section.
xdim2 = [-.953/(c2*2),.953/(c2*2)]
# This can be confirmed by plotting, but it should be noted that for the
# 'rectBox' routine, the mesh y-coordinates will go from -c2/2 -> c2/2, and the
# x-coordinates will go from xdim2[0]*c -> xdim2[1]*c. Therefore the box's
# overall dimesions should be from 0.53 in x 0.953 in. Next we will generate
# the airfoil box:
af2 = Airfoil(c2,name='box')
# Now let's make all of the laminate objects we will need for the box beam. In
# this case it's 4:
n_i_Lay1 = [6]
m_i_Lay1 = [2]
th_Lay1 = [0.]
lam1_Lay1 = Laminate(n_i_Lay1, m_i_Lay1, matLib, th=th_Lay1)
lam2_Lay1 = Laminate(n_i_Lay1, m_i_Lay1, matLib, th=th_Lay1)
lam3_Lay1 = Laminate(n_i_Lay1, m_i_Lay1, matLib, th=th_Lay1)
lam4_Lay1 = Laminate(n_i_Lay1, m_i_Lay1, matLib, th=th_Lay1)
# Assemble the laminates into an array. Refer back to the documentation to
# remind yourself of stacking sequence direction, etc. It should be noted that
# the first laminate is at the top of the box cross-section. The next is the
# left laminate, and so on in a counter-clockwise direction.
laminates_Lay1 = [lam1_Lay1,lam2_Lay1,lam3_Lay1,lam4_Lay1]
# Create the cross-section vector:
xsect_Lay1 = XSect(2,af2,xdim2,laminates_Lay1,matLib,typeXSect='rectBox',\
    meshSize=2)
# Create the cross-section object. Note that since we aren't specifying the
# reference axis, it will automatically be set at the shear center. Since this
# cross-section is simple, this will still be at the origin.
xsect_Lay1.xSectionAnalysis()
# Let's look at the stiffness matrix:
xsect_Lay1.printSummary(stiffMat=True)
# Since this cross-section is simple, we can analytically calculate some of the
# simpler parameters. For example, the 3,3 entry is just E1*A. Similarly, the
# bending stiffnesses of the cross-section are just E1*I_xx and E1*I_yy. Try
# calculating them on your own to verify this! Note that this will only work
# since all of the fibers have a 0 degree orientation. Let's apply a load and
# see how it behaves!
# Force Resultant Vector:
force2 = np.array([0.,0.,0.,0.,0.,100.])
# Calculate the effects of the force resultant
xsect_Lay1.calcWarpEffects(force=force2)
# Plot the normal stress
xsect_Lay1.plotWarped(figName='Layup 1 Box Beam Sigma_33 Stress',\
    warpScale=100,contour='sig_33',colorbar=True)
# Now the shear 13 stress!
xsect_Lay1.plotWarped(figName='Layup 1 Box Beam Sigma_13 Stress',\
    warpScale=100,contour='sig_13',colorbar=True)

# Look at the differences in magnitudes of the stress between the two plots.
# Notice anything? There is virtually no normal stress, but A LOT of shear
# stress. This makes sense though since we only applied a torque to the cross-
# section. Note that this warping profile is vero common for any box type cross
# section. Let's try one more slightly more complex shape:

# NACA 2412 BOX BEAM
# ==================

# Now let's mesh a NACA2412 box beam. We will use the last of the supported
# meshing routines for this. This is the less restrictive than the 'rectBox'
# routine, and has different laminate mesh interfaces. This time we will also
# make a slightly more interesting mesh using unbalanced and unsymetric
# laminates. First let's initialize the airfoil shape:
# Initialize a chord length of four inches
c3 = 4.
# Initialize the non-dimesional locations for the airfoil points to be
# generated:
xdim3 = [.15,.7]
# Create the airfoil object:
af3 = Airfoil(c3,name='NACA2412')
# Create the laminates to make up the cross-section
n_i_1 = [1,1,1,1,1,1]
m_i_1 = [2,2,2,2,2,2]
th_1 = [-15,-15,-15,-15,-15,-15]
lam1 = Laminate(n_i_1, m_i_1, matLib, th=th_1)
n_i_2 = [1,1,1,1,1,1]
m_i_2 = [2,2,2,2,2,2]
th_2 = [15,-15,15,-15,15,-15]
lam2 = Laminate(n_i_2, m_i_2, matLib, th=th_2)
n_i_3 = [1,1,1,1,1,1]
m_i_3 = [2,2,2,2,2,2]
th_3 = [15,15,15,15,15,15]
lam3 = Laminate(n_i_3, m_i_3, matLib, th=th_3)
n_i_4 = [1,1,1,1,1,1]
m_i_4 = [2,2,2,2,2,2]
th_4 = [-15,15,-15,15,-15,15]
lam4 = Laminate(n_i_4, m_i_4, matLib, th=th_4)
# Organize the laminates into an array
laminates_Lay3 = [lam1,lam2,lam3,lam4]
# Create the cross-section object and mesh it
xsect_Lay3 = XSect(4,af3,xdim3,laminates_Lay3,matLib,typeXSect='box',meshSize=2)
# Run the cross-sectional analysis. Since this is an airfoil and for this,
# symmetric airfoils the AC is at the 1/c chord, we will put the reference axis
# here
xsect_Lay3.xSectionAnalysis(ref_ax=[0.25*c3,0.])
# Let's see what the rigid cross-section looks like:
xsect_Lay3.plotRigid()
# Print the stiffness matrix
xsect_Lay3.printSummary(stiffMat=True)
# Create an applied force vector. For a wing shape such as this, let's apply a
# semi-realistic set of loads:
force3 = np.array([10.,100.,0.,10.,1.,0.])
# Calculate the force resultant effects
xsect_Lay3.calcWarpEffects(force=force3)
# This time let's plot the max principle stress:
xsect_Lay3.plotWarped(figName='NACA2412 Box Beam Max Principle Stress',\
    warpScale=10,contour='MaxPrin',colorbar=True)

Tutorial 4 - Using Superbeams to Conduct Simple Beam Analysis

This tutorial exposes users to conducting beam analysis using the SuperBeam class. It also exposes users to the Model class for the first time.

# =============================================================================
# AEROCOMBAT TUTORIAL 4 - Using Super Beams to Conduct Analysis
# =============================================================================

# IMPORT SYSTEM PACKAGES
# ======================
import sys
import os

sys.path.append(os.path.abspath('..'))

# IMPORT AEROCOMBAT CLASSES
# =========================
from AeroComBAT.Structures import MaterialLib, Laminate, XSect, SuperBeam
from AeroComBAT.Aerodynamics import Airfoil
from AeroComBAT.FEM import Model

# IMPORT NUMPY MODULES
# ====================
import numpy as np
import mayavi.mlab as mlab

# ADD MATERIALS TO THE MATERIAL LIBRARY
# =====================================
# Create a material library object
matLib = MaterialLib()
# Add material property from Hodges 1999 Asymptotically correct anisotropic
# beam theory (Imperial)
matLib.addMat(1,'AS43501-6','trans_iso',[20.6e6,1.42e6,.34,.3,.87e6,.1],.005)
# Add material property from Hodges 1999 Asymptotically correct anisotropic
# beam theory (Imperial)
matLib.addMat(2,'AS43501-6*','trans_iso',[20.6e6,1.42e6,.34,.42,.87e6,.1],.005)
# Add an aluminum material (SI)
matLib.addMat(3,'AL','iso',[71.7e9,.33,2810],.005)

# CREATE A WRAPPED RECTANGULAR BOX-BEAM CROSS-SECTION
# ===================================================
# Layup 1 Box beam (0.5 x 0.923 in^2 box with laminate schedule [0]_6)

# Before we make a beam, we must first make the cross-section of that beam. We
# are going to start with a cross-section we used in the third tutorial.
c2 = 0.53
# Establish the non-dimensional starting and stopping points of the cross-
# section.
xdim2 = [-.953/(c2*2),.953/(c2*2)]
# Generate the airfoil box:
af2 = Airfoil(c2,name='box')
# Now let's make all of the laminate objects we will need for the box beam. In
# this case it's 4:
n_i_Lay1 = [6]
m_i_Lay1 = [2]
th_Lay1 = [0.]
lam1_Lay1 = Laminate(n_i_Lay1, m_i_Lay1, matLib, th=th_Lay1)
lam2_Lay1 = Laminate(n_i_Lay1, m_i_Lay1, matLib, th=th_Lay1)
lam3_Lay1 = Laminate(n_i_Lay1, m_i_Lay1, matLib, th=th_Lay1)
lam4_Lay1 = Laminate(n_i_Lay1, m_i_Lay1, matLib, th=th_Lay1)
# Assemble the laminates into an array.
laminates_Lay1 = [lam1_Lay1,lam2_Lay1,lam3_Lay1,lam4_Lay1]
# Create the cross-section vector:
xsect_Lay1 = XSect(2,af2,xdim2,laminates_Lay1,matLib,typeXSect='rectBox',\
    meshSize=2)
# Create the cross-section object.
xsect_Lay1.xSectionAnalysis()
# Having created the cross-section, we can now generate a superbeam. A
# superbeam is just a collection of beam elements. In other words, a superbeam
# is just there to fascilitate beam meshing and other pre/post processing
# benefits. In order to make a superbeam, we need to initialize a few things.
# First, let's initialize the starting and stopping location of the beam:
x1 = np.array([0,0,0])
x2 = np.array([0,0,4])
# Initialize a superbeam ID
SBID = 1
# Next we need to initialize the number of elements the superbeam should mesh:
noe = 20
# Now let's make the superbeam
sbeam1 = SuperBeam(SBID,x1,x2,xsect_Lay1,noe)
# In order to analyze this beam, we'll need to add it to a finite element
# model. First let's make a finite element model!
model = Model()
# Easy right? Now let's add the superbeam to the model.
model.addElements([sbeam1])
# Now that our beam is loaded into the FEM, let's visualize it!
model.plotRigidModel(numXSects=8)
# First let's constrain the beam at it's root. Since we only added one
# superbeam and we never specified a starting node ID, we know that the first
# node ID is actually 1! So when we constrain the model, we can just select:
model.applyConstraints(1,'fix')
# There are two supported keywords for constraints, either 'fix' or 'pinned'.
# If the user wanted to apply a different constraint, they can just enter the
# degrees of freedom to constrain on the model. This can be done by supplying
# an array such as: [1,2,3,5]. Now let's apply a load. We will make two load
# cases. In the first case, we are going to apply a simple tip load:
load1 = {21:np.array([100.,100.,0.,0.,0.,100.])}
# We can also create a function for a distributed load:
def load2(x):
   vx = (1/10)*10*x[2]**2-7*x[2]-2.1
   vy = 10*x[2]**2-7*x[2]
   pz = 0
   mx = 0
   my = 0
   tz = (10*x[2]**2-7*x[2])/10+3*x[0]**2
   return np.array([vx,vy,pz,mx,my,tz])
# Ok now let's add these loads to the model:
model.applyLoads(1,F=load1)
# Notice that when I applied a tip load, I did it using the argument 'F'. When
# we apply a distributed load function, we use the argument 'f' instead.
model.applyLoads(2,f=load2,allElems=True)
# Now with constraints and loads, we can run a static analysis! Let's run the
# first load case.
model.staticAnalysis(1,analysis_name='tip load')
# Let's see what results we get:
model.plotDeformedModel(analysis_name='tip load',figName='Tip Load Analysis',\
    numXSects=8,contour='VonMis',contLim=[0,1e5],warpScale=50,displScale=10)
# Now let's try analyzing the distributed load:
model.staticAnalysis(2,analysis_name='distributed load')
# Let's see what results we get for the distributed load:
model.plotDeformedModel(analysis_name='distributed load',\
    figName='Distributed Load Analysis',numXSects=8,contour='VonMis',\
    contLim=[0,1e5],warpScale=50,displScale=10)
# Really quickly, let's discuss some of the keywords used in the analysis and
# plotting. The analysis_name designates the name where the results should be
# stored for the beam. Therefore if you want to keep multiple results stored
# at once you can give it a name. Otherwise this will always result to the
# default name. figName is just the name of the MayaVi figure, and numXSects
# designates how many cross-sections should be plotted. Note that the more
# cross-sections are plotted, the slower the plotting process will be. Both
# contour, contLim and warpScale were discussed in Tutorial 3. displScale is
# the scalling factor applied to the beam displacements and rotations.
# We can also run a normal modes analysis on the beam as well: having already
# applied the constraints, we can run:
model.normalModesAnalysis(analysis_name='normal modes')
# For now, the frequencies can be accessed via the model attribute 'freqs'
Frequencies = model.freqs
# Let's plot the first mode:
model.plotDeformedModel(analysis_name='normal modes',\
    figName='Normal Modes 1',numXSects=10,contour='none',\
    warpScale=1,displScale=2,mode=1)
# How about the second?
model.plotDeformedModel(analysis_name='normal modes',\
    figName='Normal Modes 2',numXSects=10,contour='none',\
    warpScale=1,displScale=5,mode=2)
# How about the third? It happens to be the torsional mode.
model.plotDeformedModel(analysis_name='normal modes',\
    figName='Normal Modes 3',numXSects=10,contour='none',\
    warpScale=1,displScale=1,mode=3)

Tutorial 5 - Creating a Wing and Conducting a Flutter Analysis

This tutorial exposes users to creating a wing and conducting flutter analysis on it. Note that every analysis a superbeam can conduct can also be conducted on a wing.

# =============================================================================
# AEROCOMBAT TUTORIAL 4 - Using Super Beams to Conduct Analysis
# =============================================================================

# IMPORT SYSTEM PACKAGES
# ======================
import sys
import os
sys.path.append(os.path.abspath('..'))

# IMPORT NUMPY PACKAGES
# =====================
import numpy as np
import matplotlib.pyplot as plt

# IMPORT AEROCOMBAT CLASSES
# =========================
from AeroComBAT.Structures import MaterialLib
from AeroComBAT.AircraftParts import Wing
from AeroComBAT.FEM import Model

# ADD MATERIALS TO THE MATERIAL LIBRARY
# =====================================
# Create a material library object
matLib = MaterialLib()
# Add an aluminum material (SI)
matLib.addMat(1,'AL','iso',[68.9e9,.33,2700*2],.00025)
# Add an soft material material (SI)
matLib.addMat(2,'Weak_mat','iso',[100,.33,10],.00025)
# Add material property from Hodges 1999 Asymptotically correct anisotropic
# beam theory (SI)
matLib.addMat(3,'AS43501-6*','trans_iso',[142e9,9.8e9,.34,.42,6e9,20000],0.0005)

# CREATE THE WING
# ===============
# Define the chord length of the model
c = .076
# Define the chord length (this will be the hight of the box beam)
ctip = 0.0076+.001
# Since the wing isn't tapered
croot = 0.0076+.001
# Define the non-dimensional starting point of the cross-section
x1 = -0.039/croot/2
# Define the non-dimensional ending point of the cross-section
x2 = 0.039/croot/2
# Define the span of the beam
span = 0.76
# Define the starting and stopping point of the wing structure
p1 = np.array([c/2,0.,0.])
p2 = np.array([c/2,span,0.])
# Define the non-dimesional locations of the ribs in the wing.
Y_rib = np.linspace(0.,1.,2)

# Initilize the layup schedule for the cross-section.
n_ply = [4,4,4,4]
# Initilize the material ID corresponding to an orientation in the
# cross-section
m_ply = [1,1,1,1]
# Initilize the orientations used in the box beam.
th_ply = [0,0,0,0]

# Define the number of orientations used in each laminate. In this case, we'll
# just use one.
n_orients = 1
# Define the number of laminates per cross-section. In this case since we are
# using a box beam, it will be 4.
n_lams = 4
# Define the type of cross-section. In this case it'll be the 'rectBox' mesh
typeXSect = 'rectBox'
# Define the number of elements per unit length are to be used. The structure
# will have then 120*span beam elements.
noe_dens = 120
# Define the aditional vector required to define the orientation of the beam.
# In this case, we'll have it pointing down the length of the wing.
chordVec=np.array([1.,0.,0.])
# Create the wing object. For more information about some of the input
# parameters see the AeroComBAT documentation
wing1 = Wing(1,p1,p2,croot,ctip,x1,x2,Y_rib,n_ply,m_ply,matLib,name='box',\
    noe=noe_dens,chordVec=chordVec,ref_ax='shearCntr',n_orients=n_orients,\
    n_lams=n_lams,typeXSect=typeXSect,meshSize=2,th_ply=th_ply)
# This is an optional step for ease of programming. Since only one wing section
# was created and the wing isn't tapered, there is only one superbeam which
# contains all of the beam elements in the wing.
sbeam1 = wing1.wingSects[0].SuperBeams[0]

# ADD A LIFTING SURFCE TO THE WING
# ================================
# Define the root leading edge location
x1 = np.array([0,0.,0.])
# Define the root trailing edge location
x2 = np.array([c,0.,0.])
# Define the tip trailing edge location
x3 = np.array([c,span,0.])
# Define the tip leading edge location
x4 = np.array([0,span,0.])
# Determine the number of boxes to be used in the spanwise direction
nspan = 36*2
# Determine the number of boxes to be used in the chordwise direction
nchord = 10
# Add the lifting surface to the wing.
wing1.addLiftingSurface(1,x1,x2,x3,x4,nspan,nchord)

# MAKE THE FINITE ELEMENT MODEL (FEM)
# ===================================
model  = Model()
# Add the aircraft wing to the model
model.addAircraftParts([wing1])
# Apply the constraint for the model
model.applyConstraints(0,'fix')
# Plot the rigid wing in the finite elment model
model.plotRigidModel(numXSects=10)
# Conduct a normal modes analysis. Since the normal mode shapes are used in the
# flutter analysis it is good to run this ahead of time to make sure you are
# selecting enough mode shapes to include any relevant torsional and bending
# mode shapes.
model.normalModesAnalysis()
# Save the frequencies of the modal analysis
freqs = model.freqs

# IMPORT NASTRAN RESULTS
# ======================
# The same analysis was conducted on a plate model in NX NASTRAN to verify the
# results produced by AeroComBAT
NASTRAN = np.genfromtxt('NASTRANFlutterResults.csv', delimiter=',')
UNAST = NASTRAN[:,0]
Damp1 = NASTRAN[:,1]
Freq1 = NASTRAN[:,2]
Damp2 = NASTRAN[:,3]
Freq2 = NASTRAN[:,4]
Damp3 = NASTRAN[:,5]
Freq3 = NASTRAN[:,6]
Damp4 = NASTRAN[:,7]
Freq4 = NASTRAN[:,8]
Damp5 = NASTRAN[:,9]
Freq5 = NASTRAN[:,10]
Damp6 = NASTRAN[:,11]
Freq6 = NASTRAN[:,12]

# CONDUCT FLUTTER ANALYSIS
# ========================
# Whenever a flutter analysis is conducted, several quantities need to be
# defined. The first is an array of free-stream airspeeds:
U_vec = np.linspace(1,342,100)
# In addition, a vector of trial reduced frequencies must be initialized. Keep
# in mind that the because the non-looping pk method is used, a wide range of
# reduced frequencies may need to be used.
kr_vec = np.array([0.,1e-06,1e-04,.001,.01,.05,.1,.5,1.,5.,10.,50])*10
# A vector of mach numbers must also be used. These should be kept close to
# The suspected flutter mach number. If mach numbers greater than 0.7 are used,
# is is likely doublet lattice method is no-longer valid.
M_vec = [0.]*len(kr_vec)
# Initialize the sea level density
rho_0 = 1.225
# Determine the number of modes to be used
nmodes = 6
# Run the flutter analysis. Depending on how many reduced frequencies and
# velocities sampled, this could take anywhere from a 30 seconds to 2 minutes.
model.flutterAnalysis(U_vec,kr_vec,M_vec,c,rho_0,nmodes,symxz=True,g=0.0)

# POST-PROCESS THE FLUTTER RESULTS
# ================================
# Note that in these figures, the dashed lines are the NASTRAN results, whereas
# the solid lines are the AeroComBAT results.
cvec = ['b','g','r','c','m','y']
plt.figure(1)
plt.hold(True)
for PID, point in model.flutterPoints.iteritems():
    plt.plot(U_vec,point.gamma,color=cvec[PID],label='Mode '+str(PID+1))
plt.legend(loc=2)
plt.ylim([-1,1])
plt.plot(UNAST,Damp1,str(cvec[0])+'--',label='NASTRAN Mode 1')
plt.plot(UNAST,Damp2,str(cvec[1])+'--',label='NASTRAN Mode 2')
plt.plot(UNAST,Damp3,str(cvec[2])+'--',label='NASTRAN Mode 3')
plt.plot(UNAST,Damp4,str(cvec[3])+'--',label='NASTRAN Mode 4')
plt.plot(UNAST,Damp5,str(cvec[4])+'--',label='NASTRAN Mode 5')
plt.plot(UNAST,Damp6,str(cvec[5])+'--',label='NASTRAN Mode 6')
plt.title('Damping of the Wing Modes')
plt.xlabel('Free-stream airspeed, m/s')
plt.ylabel('Damping, g')
plt.grid(True)
plt.hold(False)

plt.figure(2)
plt.hold(True)
for PID, point in model.flutterPoints.iteritems():
    plt.plot(U_vec,point.omega,color = cvec[PID],label='Mode '+str(PID+1))
plt.legend(loc=1)
plt.plot(UNAST,Freq1,str(cvec[0])+'--',label='NASTRAN Mode 1')
plt.plot(UNAST,Freq2,str(cvec[1])+'--',label='NASTRAN Mode 2')
plt.plot(UNAST,Freq3,str(cvec[2])+'--',label='NASTRAN Mode 3')
plt.plot(UNAST,Freq4,str(cvec[3])+'--',label='NASTRAN Mode 4')
plt.plot(UNAST,Freq5,str(cvec[4])+'--',label='NASTRAN Mode 5')
plt.plot(UNAST,Freq6,str(cvec[5])+'--',label='NASTRAN Mode 6')
plt.title('Frequency of the Wing Modes')
plt.xlabel('Free-stream airspeed, m/s')
plt.ylabel('Mode Frequency, Hz')
plt.grid(True)
plt.hold(False)

# The following figures demonstrate how the damping and frequencies are
# interpolated.
Uind = 80
point1 = model.flutterPoints[0]
point2 = model.flutterPoints[1]
point3 = model.flutterPoints[2]
point4 = model.flutterPoints[3]
omegaAeros = point1.omegaAeroDict[U_vec[Uind]]
omegaRoots1 = point1.omegaRootDict[U_vec[Uind]]
omegaRoots2 = point2.omegaRootDict[U_vec[Uind]]
omegaRoots3 = point3.omegaRootDict[U_vec[Uind]]
omegaRoots4 = point4.omegaRootDict[U_vec[Uind]]
gammas1 = point1.gammaDict[U_vec[Uind]]
gammas2 = point2.gammaDict[U_vec[Uind]]
gammas3 = point3.gammaDict[U_vec[Uind]]
gammas4 = point4.gammaDict[U_vec[Uind]]

plt.figure(3)
plt.hold(True)
plt.plot(omegaAeros,omegaAeros,str(cvec[0])+'o-',label='omega_aero')
plt.plot(omegaAeros,omegaRoots1,str(cvec[1])+'o-',label='omega_root_1')
plt.plot(omegaAeros,omegaRoots2,str(cvec[2])+'o-',label='omega_root_2')
plt.plot(omegaAeros,omegaRoots3,str(cvec[3])+'o-',label='omega_root_3')
plt.plot(omegaAeros,omegaRoots4,str(cvec[4])+'o-',label='omega_root_4')
plt.legend(loc=2)
plt.ylim([0,20000])
plt.xlim([0,25000])
plt.xlabel('Aerodynamic frequency, rad')
plt.ylabel('Root frequency, rad')
plt.title('Interpolation of Root Requencies at V=%4.2f m/s'%(U_vec[Uind]))
plt.grid(True)
plt.hold(False)

plt.figure(4)
plt.hold(True)
plt.plot(omegaAeros,gammas1,str(cvec[1])+'o-',label='gamma_root_1')
plt.plot(omegaAeros,gammas2,str(cvec[2])+'o-',label='gamma_root_2')
plt.plot(omegaAeros,gammas3,str(cvec[3])+'o-',label='gamma_root_3')
plt.plot(omegaAeros,gammas4,str(cvec[4])+'o-',label='gamma_root_4')
plt.legend(loc=3)
plt.ylim([-1,.1])
plt.xlim([0,25000])
plt.xlabel('Aerodynamic frequency, rad')
plt.ylabel('Damping (g)')
plt.title('Interpolation of Root Damping at V=%4.2f m/s'%(U_vec[Uind]))
plt.grid(True)
plt.hold(False)