Computational Engineering Analysis (AERO2463): MATLAB Project 2019
3.2 Space Science and Rocket Design: Orbital Manoeuvres
By changing the velocity of the lunar lander, it is possible to change the orbit. Changing the velocity allows the lunar lander to travel from low lunar orbit, down to the lunar surface, and back again. The transfer is carried out using an elliptical transfer orbit. We will choose the simplest possible transfer orbit, an elliptical orbit with the apogee at low lunar orbit, and the perigee at the radius of the Moon’s surface. This orbit will directly transfer from the low lunar orbit to the surface through half an orbit (see Figure 4 below).
Figure 4: Schematic of Low Lunar and Transfer Orbits
To use the rocket equation to determine the volume of fuel burnt, you will need to know the delta-V required to get onto the elliptical orbit and the final velocity of the orbit. The final velocity will be necessary to determine the last burn to stop at the Moon’s surface and return to the stable low lunar orbit at the end.
The first step is to determine the semi-major axis of the transfer orbit ();
2
and and the radii of the initial and final orbit, respectively.
The period of the elliptical transfer orbit () can is:
4
(2)
(3)
Where is the universal gravitational constant, and is the mass of the central object of the orbit, in
our case, the Moon.
Using the orbital period and the semi-major axis, the initial velocity of the transfer orbit ( ) is:
2 2 1
(4)
Knowing the initial transfer velocity, the delta-V (Δ) to reach the orbit can be found as the difference between the initial velocity on the transfer orbit, and the velocity before going to the transfer orbit.
(5)
For the first chronological burn, from low lunar orbit to the Moon, will be the velocity of the low lunar orbit. For the burn taking off from the Moon, 0. The rotation of the Moon is assumed to be negligible.
Finally, the velocity at the end of the elliptical orbit ( ) will be needed to determine the delta-V for landing on the Moon and returning to a stable low lunar orbit.
Page 5
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
2 2 1
(6)
To land on the Moon, delta-V is simply the final velocity of the transfer orbit, and the landing velocity is zero. For returning to low lunar orbit, delta-V is the difference between the velocity final velocity of the transfer orbit and velocity of the circular low lunar orbit.
3.3 Space Science and Rocket Design: Tsiolkovsky Rocket Equation
In its simplest terms, a rocket engine burns a mass of fuel, sending a stream of high-velocity gas in one direction and thrusting the rocket in the other direction, all in accordance to Newton’s Third Law and conservation of momentum. The mass of fuel (Δ) required to produce a specific delta-V or change in velocity for the rocket (Δ), can be directly related to the final mass of the rocket after the burn () and the effectiveness of the engine.
1
Here is the Earth’s gravitational acceleration. is chosen as a standard.
(7)
is called the specific impulse of the engine. The specific impulse quantifies the total impulse (or change in momentum) delivered per unit of propellant used. The study and optimisation of specific impulse is a
vast, and fascinating subject, and far beyond the requirements of this project. You simply need to understand it as an input to (7), and that bigger is better. Beware however than many high specific impulse engines are not suitable for space travel, or have only ever been tested in labs, simulations or Star Trek. When searching for a suitable engine, consider traditional chemical rocket engines such as the Raptor which has a specific impulse of 331 s in a vacuum (SpaceFlight101, 2019).
3.4 Space Science and Rocket Design: Mission
The whole mission for the lunar lander is shown chronologically below. You will need to work through the mission backwards, starting with only the dry mass of the ascent stage to determine the total fuel required.
3.5
1. 2. 3. 4.
5.
Start by orbiting the Moon at the low lunar orbit,
Fire the engines to reach a transfer orbit leading to the Moon’s surface,
Fire the engines to come to a stop at the Moon’s surface,
Fire the engines to take off from the Moon’s surface onto a transfer orbit leading back to the low lunar orbit,
Fire the engines to change speed to reach the correct insertion velocity for the low lunar orbit.
Finite Element Analysis: Material Selection
To determine whether your lander with a truss structure is safe under the loads, you will need to carry out static finite element analysis for each of the two load cases. For this analysis, a simple truss element analysis is sufficient to check if the stresses are below the allowable stress.
Your first decision will be the material you choose. We recommend using routine engineering materials. You should select and compare three materials, from these three categories, metals/alloys, polymers, and composites. You need to use some engineering judgement in your material choice.
You will need to find Young’s modulus and, density for the finite element analysis. The yield strength or ultimate strength depending on the requirements you choose for defining the failure of your truss. Ideally, you will select a lighter material, with the highest failure (yield/ultimate) strength.
Page 6
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
As a starting point, for metal, consider Mild steel with a density of 7850 kg/m3, Young’s Modulus of 200 GPa, Yield Strength of 200 MPa, and Ultimate Strength of 340 MPa (Gere and Goodno, 2009), or Aluminium with a density of 2700 kg/m3, Young’s Modulus of 70 GPa, Yield Strength of 270 MPa, and Ultimate Strength of 310 MPa (Gere and Goodno, 2009).
3.6 Finite Element Analysis: Node and Element Definitions
The truss you will use for your lander will be made up of discrete nodes connected by elements. A node is a distinct, infinitesimal point in 3D space. An element generally is a shape connecting two or more nodes. For the analysis in this assignment, we will use ‘truss’ elements which are straight lines joining two nodes. The same node can have multiple elements connected to it. A force or boundary condition can only apply at a node. Stresses and strains can only be computed for elements, as elements have finite lengths. We will provide you with all functions to calculate the stresses through finite element analysis. For details of the internal working process of FEA, see the appendix of this document and refer to the RMIT FEA course materials of AERO2357 or MIET1084.
4
Basic Overview and Steps
The “concept of operations” (overview) for the entire program is described below. We recommend you progressively build up the project one activity at a time, completing each function before using it in more complex steps. Activities to accomplish this are provided in Section 6.
(1) You are given a spreadsheet (excel file) < Sxxxxxxx_Syyyyyyy_Design_File.xlsx>. The main script template reads this file to get your inputs. Your main script will perform all of the following steps.
(2) Read the data from the spreadsheet (code supplied in RunAnalysis, marked “do not edit”).
(3) Starting with a truss design and chosen cross-sectional area of the truss members, complete and
use RunAnalysis to perform steps 6 and 7 below.
(4) To find the optimum area, define a broad range of cross-section areas to test to find the minimum cross-sectional area and hence minimum fuel consumption, with peak stress remaining with the allowable limit. This is the “coarse search”.
(5) Determine the length of each element using ElementLengths then combine and use to determine the total mass of the truss
(6) Calculate the fuel usage for the lunar landing mission using component A – Orbital Mechanics as described below.
(7) Determine maximum stress for the truss over all two load cases using component B – Implicit Finite Element Analysis (using provided function RunFEA).
(8) Plot the stress versus axial cross-section and required fuel mass versus axial cross-section for the rough estimates of cross-sectional area at which the maximum stress is nearly equal to the allowable stress.
(9) Use polyfit to determine a polynomial fit between the maximum stress and the cross- sectional area. You can fit a polynomial up to order 8, but use the lowest order that gives a reasonable fit (remember Assignment 4.1). Lower order means fewer data points are needed, so the computationally expensive FEA need be run fewer times.
(10) Identify the cross-section areas from your coarse search that are closest to the limit of maximum allowable stress. One is the smallest area under the maximum permissible stress, and one will be the nearest largest area above the maximum permissible stress.
Page 7
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
(11) Using the closest limits identified and the polynomial fit of the data, use the bisection method to determine the minimum allowable cross-section area (remember Assignment 1.2), and plot the progression of cross-section area vs maximum stress guesses. This is the “fine search”.
(12) Using supplied function DrawModel from within RunAnalysis, draw a representation of the truss model, e.g. Figure 5 below.
Figure 5: An example of the images expected from DrawModel with descent module fuel mass and ascent module mass
(13) Using this minimum allowable cross-section area, determine the dry mass, call RunAnalysis once more with plotting on to determine the fuel mass used and plot the stresses for each load case (using supplied function PostProcessFEA).
(14) Repeat your analysis for a total of three materials to determine which is most suitable for your truss design (by running your script once per material).
Component A – Orbital Mechanics (inside RunAnalysis function)
We will perform the orbital mechanics calculations backwards from the phase when the lander retuns to the low lunar orbital from the transfer orbit at the end of the mission.
(1) Determine the orbital velocity of the low lunar circular orbit with altitude (radius) defined in metres as the average of your student numbers divided by 10 (Orbital_Velocity)
(2) Use equations of orbital manoeuvres to determine the velocity of the ascent module at the end of the transfer orbit returning from the Moon’s surface to the low lunar orbit. Draw a plot of the orbital path, including the Moon’s surface and the path of the low lunar orbit (Manoeuvre).
(3) Calculate the fuel mass required for the ascent module to slow from the transfer orbit, back into the low lunar orbit using the rocket equation (Rocket_Eqn).
(4) Calculate the fuel mass required to produce enough velocity difference (delta-V) to reach the transfer orbit from the surface of the Moon for the ascent module using the rocket equation (Rocket_Eqn).
(5) Calculate the velocity at the end of the transfer orbit from the low lunar orbit to the Moon using orbital manoeuvre equations. Draw a plot of the orbital path, including the Moon’s surface and the path of the low lunar orbit (Manoeuvre).
Page 8
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
Figure 6: Example of the figures you should plot from the orbital mechanics calculations
(6) Calculate the fuel mass required to slow from the velocity at the end of the descent transfer orbit
to zero velocity at the surface of the Moon using the rocket equation (Rocket_Eqn).
(7) Calculate the fuel mass required to produce enough velocity change (delta-V) to enter the
descent transfer orbit from the low lunar orbit using the rocket equations (Rocket_Eqn).
Component B – Implicit Finite Element Analysis (inside RunAnalysis function)
(1) Run DrawModel to draw the launch phase lander truss model with all dry and fuel masses,
(2) Run RunFEA to determine the stress in each element for the two load cases with different fuel
masses determined using the orbital mechanics in component A.
(3) Determine the absolute maximum (peak) stress from the FEA and compare this to the maximum allowable stress for the selected materials.
Remember: This assignment is very open-ended. There are several choices you will have to make along the way and briefly justify in the report. Some of these choices may include:
Design of the truss structure configuration,
Material for the truss out of a metal/alloy, polymer or composite,
Cross-sectional areas of each element in the truss. Remember this is a truss analysis, so the actual cross-section shape or profile is not needed.
Allowable stress in the material, do you use yield strength or ultimate strength or failure strength, do you need to include a safety factor?
Type of engine for the lunar lander. Please choose a chemical engine that is reliable. We do not expect you to understand the limitations of some more efficient, but as yet un-proved engines.
Where are the locations of the fuel point masses?
Do you use a multi-stage lunar lander, leaving the descent module truss on the moon, or do you take
it back with you?
You may make use of previous assignment’s code for optimising the design or minimising the required fuel. However, be careful when repurposing old code as there may be different required inputs and outputs, or clashes with variable names.
Page 9
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
5 Load Cases
Two pseudo-static load cases need to be considered for your lunar lander. These are all included in Run FEA, and you will only need to interpret the results.
(1)
(2)
6
Launch from Earth:
Full fuel mass and both ascent and descent stages present,
You chose the attachment points of the ascent module by setting boundary conditions on nodes you
specify,
Peak global acceleration imparted by the launch vehicle is 40 ms-2 in the negative z-direction. This
implies if the rocket is accelerating at +40 ms-2, then in the pseudo-static load case you have to apply has an acceleration of -40 ms-2 on every element, and BCs where the truss is attached to the rocket.
Landed on the Moon:
Fuel mass is after landing burn,
Feet of the lander are supported on the Moon by assigning boundary conditions,
Global acceleration is the lunar gravity on the Moon’s surface, which is 1.62 ms-2 in the downward
(negative z-direction).
Development Process
Developing a large and complex project can be a difficult task. We recommend progressively developing one stage at a time, checking each behaves correctly before progressing onto the next one. The following activities are defined to assist you with this.
6.1
6.2
Activity 1 – Input data processing
(1) Use the provided spreadsheet (excel file) to develop a truss design by data input in the specified by the template.
(2) Read the data from the spreadsheet into Matlab using RunAnalysis.
(3) Use DrawModel to confirm that the truss is designed as you intended.
Activity 2 – Analysis of the truss structure for a given cross-sectional area (call RunAnalysis for a single case)
(1) Using the total length of the truss members, calculated from ElementLengths, determine the total dry mass of the descent module truss for a specific cross-sectional area.
(2) Calculate the fuel usage for the lunar landing mission using component A – Orbital Mechanics for a single cross-section area.
(3) Determine the maximum stress for the truss over both load cases using component B – Implicit Finite Element Analysis for a specific cross-section area.
6.3 Activity 3 – Analysis of the truss structure for multiple cross-sectional areas (call RunAnalysis for multiple cases)
(1) Define a coarse range of cross-sectional areas for investigating the optimal cross-section for the lunar lander truss members
(2) Call activity 2 – RunAnalysis for the various cross-sectional areas
(3) Plot the maximum stress against the cross-sectional area Page 10
6.4
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
(4) Plot the fuel mass against the cross-sectional area
Activity 4 – Derivation of a relationship between the peak stress and cross-
sectional area (using polynomial fit)
(1) Determine a polynomial fit to the coarse calculated relationship between maximum stress and cross-sectional area data from activity 3 (remember Tutorial and Assignment 4.1). Plot the polynomial fit compared to the coarse calculated values of maximum stress vs the cross- sectional area
6.5 Activity 5 – Optimisation of the truss cross-sectional area (using the bisection method)
6.6
Activity 6 – Analysis of the truss structure for different materials
(1) Repeat the full analysis for a total of three different materials, chosen from these three categories – metals/alloys, composites, and polymers. Composite properties in the fibre direction could be chosen as the relevant properties of the materials.
(1)
(2)
(3)
Add the ability to identify the cross-section areas from your coarse search in activity 3 that are closest to the limit of maximum allowable stress. One is the smallest area under the maximum permissible stress, and one will be the nearest largest area above the maximum permissible stress.
Using the closest limits identified and the polynomial fit from activity 4, use the bisection method to determine the minimum allowable cross-sectional area (remember Assignment 1.2) and plot the progression of the guesses of maximum stress vs cross-section area. This involves writing and calling the Bisection function (not supplied).
Using this minimum allowable cross-section area, determine the dry mass, re-run the orbital- mechanics calculations and FEA analysis to determine the fuel mass used and plot the stresses for each load case
Page 11
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
7 Programming
7.1 Overview
7.2
You must create a main program (script file) named Run_Matlab_Project.m
The marker assessing your program should be able to run your program with their own input
spreadsheet by changing the variable Filename.
The function header and comments provided below describe what each function must do. The requirements for the main program and the functions are also outlined below.
Main Program (Script)
The main program must be a Matlab script that you must program and write.
The main Matlab script program MUST BE NAMED Run_Matlab_Project.m
The user will edit the input at the start of the script (given in the template) before running the main program. The following variables are for you to define at the start of the script named Run_Matlab_Project.m
%% EXCEL FILE NAME from which to read Parameters, FE Model Definition, and Load
Cases:
Filename = ‘Sxxxxxxx_Sxxxxxxx_Design_File.xlsx’;
%% Material Properties
rho =; % Element Material Density (kg/m^3)
E =; % Element Material Young’s Modulus (N/m^2)
yield_strength =; % Element Material Yield/Ultimate Strength (N/m^2)
%% Parameters
AminCoarse =; % initial min cross-sectional area in m^2
AmaxCoarse =; % initial max cross-sectional area in m^2
nA =; % number of cross-sectional areas to check
n=; % Order of polynomial to fit to coarse data
safety_factor = ; % Factor by which to divide yield strength to have maximum
allowable stress
stress_allowed = yield_strength / safety_factor;
tolerance=; % Convergence tolerance for bisection method, in units of stress
(Pa)
maxiter=; % Maximum number of iterations for bisection method
NOTE: When submitted for marking, all variables should already be set and kept in the top section of the file.
Page 12
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
7.3 Matlab Functions
The main program should use the following functions that you must program.
These functions help break down the program into smaller, easier to program sections. You MUST ensure that:
1) The function names are EXACTLY the same as those specified below, and
2) Each function MUST use the inputs and outputs of the variables in the same order as defined.
You MUST follow this convention to allow assessment of the individual sections of your program.
The main program should use the following Matlab functions (all function templates are on Canvas).
Summary of Functions
Function Name
Function Description
Supplied?
Called By:
1
Orbital_Velocity
Calculates the velocity of a circular orbit.
No
RunAnalysis
2
Manoeuvre
Calculates the change in velocity required to complete an orbital transfer manoeuvre and the resultant velocity at the opposite side of the orbit, then plots the orbit (if plot_or_not==1).
No
RunAnalysis
3
Rocket_Eqn
Uses the Tsiolkovsky rocket equation to calculate the required fuel mass to perform a manoeuvre with a given delta-V.
No
RunAnalysis
4
RunAnalysis
This function brings the orbital mechanics and FEA together to determine required fuel and peak stress in a truss structure with a certain material and cross- sectional area. It will draw the structure and stress results if desired.
No
Run_Matlab_Project
5
DrawModel
Draws a representation of the truss in 3D space.
Yes
RunAnalysis
Page 13
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
6
ElementLengths
Determines the length of each element.
Yes
RunAnalysis
7
RunFEA
Calculates the displacement of nodes and stresses and strains in each element of the truss using static finite element analysis on truss elements.
Yes
RunAnalysis
8
PostProcessFEA
Creates a 3D representation of the stressed truss given FEA results for the load case.
Yes
RunAnalysis
9
Bisection
This function performs a bisection root-finding algorithm, outputting the history of the process, inclduing the converged root.
No
Run_Matlab_Project
Note: Many of these functions will require only a few lines of code to complete.
7.3.1 OrbitalVelocity
This function calculates the velocity of a circular orbit.
function V = Orbital_Velocity(M_Central,R)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Orbital_Velocity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This function calculates orbital velocity for a circular orbit at a given % radius from the center of a planetary body.
% Inputs:
% M_Central: Mass of the planetary body (kg)
% R: Radius of orbit from center of the body (m)
% Outputs:
% V: Orbital Velocity (m/s)
% Constants:
G = 6.67e-11; % Universal Gravitational Constant
7.3.2 Manoeuvre
This function calculates the change in velocity required to complete an orbital transfer and the resultant velocity at the opposite side of the orbit.
Page 14
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
function [Delta_V,V_final] = Manoeuvre(V_0, R_0, R_f,M,Phase,plot_or_not)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Manoeuvre
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the change in velocity required to complete an orbital transfer using a Hohmann transfer, and the resulting velocity at the opposite side of the new orbit.
% Then, if plot_or_not==1, polar plot the shape of the transfer orbit along
with the initial and final radii.
% Inputs:
% V_0:
% R_0:
% R_f:
% M:
% Phase:
Initial velocity
Initial transfer orbit radius
Final transfer orbit radius
Mass of central planetary body
Name of the mission phase when the manoeuvre occurs (for plot
title)
% plot_or_not: Logical 1/0 to enable/disable the orbit polar-plot
% Outputs:
% Delta_V: Change in velocity required to complete manoeuvre
% V_final: Velocity at the opposite side of the orbit resulting from the
manoeuvre
% Plot of the orbit transfer in a new figure window
% Constants:
G = 6.67e-11; % Universal Gravitational Constant
7.3.3 RocketEquation
The function uses the Tsiolkovsky rocket equation to calculate the change in fuel mass for a given velocity difference (delta-V).
function Delta_m = Rocket_Eqn(m_f, Delta_V, I_sp)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rocket_Eqn
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Using the Tsiolkovsky rocket equation, calculate the change in fuel mass % required to perform a manoeuvre with a certain change in velocity.
% Inputs:
% m_f: Final mass after manoeuvre (kg)
% Delta_V: Change in velocity of manoeuvre (m/s)
% I_sp: Specific impulse (s)
% Outputs:
% Delta_m: Change in mass (kg) (i.e. mass of fuel burned for manoeuvre)
% Constants:
g = 9.81; % Earth Surface Gravity in m/s^2 (used as constant for calculating
I_sp);
7.3.4 RunAnalysis
This function will bring the orbital mechanics and FEA together to determine the amount of fuel required for a given material and cross-sectional area.
Page 15
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
function [peak_stress,m_Descent_Module_Fuel] = RunAnalysis(Filename, A, E, rho, plot_or_not)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RunAnalysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function analyses the orbital motion and static structural response
% of a lunar lander, outputting peak stress over all load cases and the mass
% of fuel required in the descent module. Specifically it does the following:
% 1 – Reads parameters of a truss model from the excel file named by the
Filename input
% 2 – If plot_or_not is true, draws the truss model using DrawModel
% 3 – Calculates the total dry mass using ElementLengths
% 4 – Calculates fuel required for the ascent phase of the mission
% 5 – Calculates fuel required for the descent phase of the mission
% 6 – Performs static FEA for the LAUNCH load case using RunFEA (with fuel mass
and ascent module mass included)
% 7 – If plot_or_not is true, draws the stress results on the truss model for
the LAUNCH load case using PostProcessFEA
% 8 – Performs static FEA for the LANDED load case using RunFEA (with fuel mass
and ascent module mass included)
% 9 – If plot_or_not is true, draws the stress results on the truss model for
the LANDED load case using PostProcessFEA
% 10- Saves the highest peak stress over the two load cases to peak_stress
% Inputs:
% Filename:
Name of excel file to read truss model parameters from
Cross-Sectional Area to apply to all elements (m^2)
Young’s Modulus to apply to all elements (Pa)
Density to apply to all elements (kg/m^3)
% A:
% E:
% rho:
% plot_or_not: Input 1 to performing steps 7 and 9, input 0 to not
% Outputs:
% peak_stress: Maximum absolute stress over all elements and load cases
% m_Descent_Module_Fuel: Mass of fuel required in descent module
7.3.5 DrawModel
This function is supplied to you.
Creates a 3D representation of the truss in space. The function header is provided below.
function DrawModel(n_xyz, e_n, As, point_ms,Ne, Nn, Title)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DrawModel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function creates a 3D representation of the input truss structure
% defined by nodal coordinates and truss element connections.
% Elements are drawn with thickness proportional to their cross-sectional area. % Point masses are drawn as filled 2D circles with radius relative to their mass.
% Inputs:
% n_xyz: 3*Nn matrix containing xyz coodinates of each node in each
column
Page 16
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
% e_n: 3*Ne matrix containing each element’s number (row 1) and
numbers of it’s two connected nodes (rows 2 and 3)
% As: 1*Ne array containing cross-sectional area of each element
% point_ms: 1*Nn array containing point mass on each node (kg), e.g. [0 510
0] would be a 510 kg point mass on node 2
% Ne:
% Nn:
% Title:
number of elements
number of nodes
string containing load case title/name, e.g. ‘Launch’ or ‘No
Masses’
% Outputs: NONE
The output could look something like this:
Figure 7: Example of a DrawModel output with no fuel or ascent module point masses
7.3.6 Element Lengths
This function is supplied to you.
This function determines the length of each element or truss member. The function header is defined below
function Ls = ElementLengths(n_xyz,e_n,Ne)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ElementLengths
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function calculates the length of each element of a 3D truss.
% Inputs:
% n_xyz: 3*Nn matrix containing xyz coodinates of each node in each
column
% e_n: 3*Ne matrix containing each element’s number (row 1) and
numbers of its two connected nodes (rows 2 and 3)
% Ne: number of elements
% Outputs:
% Ls: 1*Ne vector containing length of each element
Page 17
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
7.3.7 RunFEA
This function is supplied to you.
Run FEA calculates the displacement of nodes and stresses and strains in each element of the truss using implicit finite element analysis on truss elements.
The function header is provided below.
function [U, Strain, Stress] = RunFEA(CaseName, Nn, Ne, n_xyz, e_n, rhos, As,
Es, point_ms, BCs, Loads, a)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RunFEA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Utilising local sub-functions, runs a static FE analysis for a particular
load case.
% 1. Call PreProcessFEA to prepare FE model from excel data
% 2. Call ApplyBCs to apply boundary conditions defined in excel
% 3. Call StaticFEA to perform static FEA, getting disps, strains, and stresses
% Inputs and Outputs to be passed in/out of the sub-functions:
% Inputs:
% CaseName
% Nn:
% Ne:
% n_xyz:
% rhos:
% As: (m^2)
% Es:
Name of the case to be processed, input as a string
number of nodes
number of elements
3*Nn matrix containing xyz coordinates of each node in each
3*Ne matrix containing each element’s number (row 1) and
numbers of it’s two connected nodes (rows 2 and 3)
column
% e_n:
% a:
or gravity) (m/s^2)
% Outputs: % U
% Strain
% Stress
globally applied acceleration (e.g. rigid body acceleration
3Nn*1 DoF-wise displacement vector
1*Ne Elemental Strain vector
1*Ne Elemental Stress vector
1*Ne array containing density of each element (kg/m^3)
1*Ne array containing cross-sectional area of each element
1*Ne array containing Young’s Modulus of each element (Pa)
1*Nn array containing point mass on each node (kg), e.g. [0
% point_ms:
510 0] would be a 510 kg point mass on node 2
% BCs: #*2 array read from excel file containing list of applied
boundary conditions (each row states a fixed DoF e.g. [node1,dof2] means node 1
is fixed in y-direction)
% Loads: #*3 array read from excel file containing list of applied
loads (each row states a loads position, direction, and magnitude e.g.
[node1,dof2,10 Newtons])
7.3.8 PostProcessFEA
This function is supplied to you.
Creates a 3D representation of the stressed truss using FEA results for the load case. The function header is provided below.
function PostProcessFEA(n_xyz, e_n, As, point_ms, Ne, Nn, U, Strain, Stress,
ScaleFactor, CaseName, t, FrameSkip, z_surface)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Page 18
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
% PostProcessFEA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function colour-plots stress on a 3D representation of the deformed
% truss structure defined by nodal coordinates and truss element connections,
% for input nodal displacement vector U and elemental stress vector S.
% Deformation is scaled by input scale factor, and the plot title includes
% the name of the load case and the value of the scale factor.
% Elements are drawn with thickness proportional to their cross-sectional
area.
% Point masses are drawn as filled 2D circles with radius relative to their
mass.
% NOTE: For static FEA results, input 1 for t, FrameSkip, and z_surface, so
they will be ignored.
% Inputs:
% n_xyz: 3*Nn matrix containing xyz coordinates of each node in each
column
% e_n: 3*Ne matrix containing each element’s number (row 1) and
numbers of it’s two connected nodes (rows 2 and 3)
% As: 1*Ne array containing cross-sectional area of each element
(m^2)
% point_ms: 1*Nn array containing point mass on each node (kg), e.g. [0
510 0] would be a 510 kg point mass on node 2
% Ne:
% Nn:
% U:
% Strain:
% Stress:
% ScaleFactor: Constant by which to scale displacements in the image
% CaseName: Title/name of the load case, entered as a string
number of elements
number of nodes
3Nn*1 column Displacement vector
1*Ne Elemental Strain vector
1*Ne Elemental Stress vector
% Explicit FEA Inputs (NOTE: input 1 for each to ignore when using for static
FEA results):
% t:
% FrameSkip:
Array of times (s) at each timestep (e.g. [0:dt:tf])
Number of timesteps skipped per frame
Anonymous function defining z-coordinate of surface wrt x
This function will implement a bisection root-finding algorithm. You should adapt the bisection method code you wrote before.
function [x_list,y_list] = Bisection(fun, a, b, maxiter, tolerance)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Bisection
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function uses the bisection method to find the root of an input
function, fun.
% z_surface:
and y (in your truss model’s coordinate system)
% Outputs: NONE
7.3.9 Bisection
% Inputs:
% fun:
% a:
% b:
% maxiter:
% tolerance:
Function of which to find roots (i.e. where f(inputs) = 0)
Initial lower limit for bisection method
Initial upper limit for bisection method
Maximum number of iterations the bisection can run
Convergence tolerance of bisection
Page 19
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
% Outputs:
% x_list: x-coordinate (“s”) of each point along the process (not
just the final one).
% y_list: y-coordinate (“fun(s)”) of each point along the process
(not just the final one).
Page 20
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
8 Testing
Testing of the Matlab program is highly recommended. Test each function individually using your own set of data. You must also test the assembled code as you integrate each function and each phase. We recommend comparing your results with the Apollo missions and using your engineering judgement to check if your results are realistic.
For example:
The Apollo missions used ~2000 kg of fuel in the ascent stage and ~8000 kg in the descent stage
(NASA report, 1969). If you require 2 kg or 2,000,000 kg of fuel, something is wrong.
A solid truss with a cross-sectional area of around 0.1 m2 is quite large if you use a standard metal. This is a “ballpark figure”; some truss designs will use a larger size, some will be much smaller, some will be variable, but if you need a cross-section area of 1.0 m2, then something may be wrong.
Velocity difference (Delta V) of order 100 – 10000 ms-1 could be considered reasonable, 3 10 ms-1 is possibly not.
Ensure your code is working at the end of each of the activities defined above before moving onto the next activity. This will assist in locating errors as simply as possible.
Page 21
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
9 Assessment
Your task is to implement the program design in Matlab (as outlined below) to calculate and optimise the weight (hence fuel consumption) for the lander, and ensuring the truss is safe (the maximum stress is below the allowable limit) under all of the load cases considered. You must present the results in a short report as detailed below. You will work in teams of two members for the project and submit one report and Matlab code for the team.
9.1 Report
You should prepare a short report (maximum 2000 words) including figures and appendix. The report needs to contain the following information. You should use the orbit radius assigned to you via your student numbers,
1. Three-dimensional plot showing the lunar lander geometry and dry point masses for one specific (example) material choice,
2. Description of each load case analysed for structural FEA,
3. Basic description and justification of fuel/engine type,
4. Two plots showing the transfer orbits from the low lunar orbit to the lunar surface and back again for one specific (example) material choice,
5. A table incorporating the velocities at each phase of the mission for a typical case of one specific (example) material choice,
6. Three-dimensional plots and charts showing the stress in the geometry for each load case for each material selected,
7. A plot of the relationship between the maximum stress and cross-sectional area for each material selected including the polynomial fits,
8. A plot from the bisection curve fit showing the progression of guesses of maximum stress versus cross-section area,
9. A plot of the relationship between fuel consumption and cross-sectional area for each material,
10. Three-dimensional plots and a table showing the peak stress in the lander geometry for the dynamic load case for each material choice,
11. Discussion of the materials chosen for the lunar lander truss structure including, Young’s modulus, density, and design criteria for assessing suitability,
12. A table incorporating the total lander mass at each phase of the mission for each material selected,
13. Total fuel mass required for each material choice,
14. List of references.
The target audience for the report is an engineer familiar with the design of a lunar lander and Matlab programming.
THE REPORT MUST BE SUBMITTED AS A PDF FILE
Any other file type (e.g. MS Word) is not accepted
NOTE: ALL files must be submitted inside a single Zip File
Page 22
9.2 Marking Criteria
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
Section
Weighting
Successfully designed a truss and displayed using DrawModel
10%
Implementation of 3 Orbital Motion functions
15% (5% × 3)
Run static FEA
10%
Run for multiple truss cross-sectional areas
10%
Produce polynomial fit
10%
Optimisation of the truss area using bisection function and calculation of minimum weight and fuel mass
15%
Implementation into main script
10%
Report
20%
Total
100%
Bonus 5% for the group that has a feasible design with the minimum fuel mass and meets the static load requirements for marking inputs
+ 5%
Note:
Functions that do not conform to the function headers given (correct inputs, outputs, and purpose) may receive mark deductions.
9.3 Result Examples
The following sample results and figures are examples of some of the images you could include in your report.
Figure 8: Truss model before launch with all fuel masses shown
Page 23
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
Figure 9: Launch load case – Maximum stresses on the truss during launch (in Pa)
Figure 10: Landed load case – Maximum stresses on the truss when landed on the Moon’s surface (in Pa)
Page 24
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
Figure 11: Plot of the orbital path from the low lunar orbit (LLO) to the Moon
Figure 12: Plot of the orbital path from the Moon to the low lunar orbit (LLO)
Page 25
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
Figure 13: Example of the polynomial fit to coarse data showing the upper and lower limits about the intersection of the absolute peak stress and the cross-sectional area
Figure 14: Example of course data showing the total fuel mass versus truss member cross-sectional area.
Page 26
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
Figure 15: Example output of the bisection method attempts to reach the optimum cross-sectional area
10 Submission
You must submit to the AERO2463 Canvas site by the due date ALL OF THE FOLLOWING ATTACHMENTS (ALL ZIPPED UP INTO A SINGLE FILE named
StudentID1_StudentID2.zip):
(The 7-Zip File Manager is installed on RMIT systems to create the Zip Package).
(1) Report, as an individual PDF File within the Zip package.
(2) Each of the Matlab functions, including the pre-built functions (supplied ones) as specified above. All of these functions must be submitted in the Zip package as individual m-files (corresponding to the function name).
(3) The main program must be submitted as a Matlab script, named Run_Matlab_Project.m and included within the Zip package.
(4) The spreadsheet (excel file) containing the truss structure inputs (e.g., the nodes and elements, material, boundary conditions, etc.) should be named Sxxxxxxx_Syyyyyyy_Design_File.xlsx and included within the Zip package.
Page 27
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
The DUE DATE is: Midnight Sunday, 27 October 2019 NOTE:
The report must be submitted as a pdf file (StudentID1_StudentID2.pdf) within a zip package at the AERO2463 Canvas site.
Please ensure the Matlab project zip file is named as
StudentID1_StudentID2.zip
before submitting
If you use any other file type (e.g. MS Word Doc), we will NOT be able to assess your submission.
11 References
Article: “Missions to the Moon”, The Planetary Society, Accessed 2019: http://www.planetary.org/explore/space-topics/space-missions/missions-to-the-moon.html.
Article: “SpaceX Raptor”, Spaceflight 101, Accessed 2019: http://spaceflight101.com/spx/spacex-raptor/. Gere, James M & Goodno, Barry J., (author.) 2013, Mechanics of materials, Eighth edition, Cengage
Learning ; [Andover : Cengage Learning, distributor], Stamford, Conn
Apollo 11 preliminary science report, NASA, SP-214, Wash, DC, 1969.
Article: “Apollo 11 Lunar Module”, NASA Space Science Data Coordinated Archive, Accessed 2019: https://nssdc.gsfc.nasa.gov/nmc/spacecraft/display.action?id=1969-059C.
Howard D. Curtis, Chapter 2 – The Two-Body Problem, In “Aerospace Engineering, Orbital Mechanics for Engineering Students” (Second Edition), Butterworth-Heinemann, 2010, Pages 61-153, ISBN 9780123747785, https://doi.org/10.1016/B978-0-12-374778-5.00002-7.
M.D. Griffin and J.R. French, “Space Vehicle Design”, 2004.
Page 28
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
12 Appendix: Theory of Implicit Finite Element Analysis
The following theory covers the work done by functions provided to you for solving static finite element analysis.
12.1 Implicit Finite Element Analysis: Truss elements
We will carry out finite element analysis using truss elements, the simplest type of element. Truss elements are only able to carry axial loads. A truss element can be equivalently be thought of as a spring with stiffness (). The stiffness is defined from the cross-sectional area (), Young’s modulus ( and the length ().
Figure 16 shows an example of a truss element.
Figure 16: A truss element The force at each end and is give as:
(8)
Or in matrix form
Or equivalently
.
(9)
(10)
(11)
(12)
is referred to as the displacement matrix, is the displacement matrix, and is the force matrix. In the absence of external forces, the force matrix can generally be obtained from the global accelerations () and the mass matrix () using Newton’s Second Law.
(13)
Page 29
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
(11) is only the stiffness matrix for one element, you will have numerous elements, and the stiffness matrices will need to be combined. Combining matrices is done by looking at the effect of each element on each of the nodes. Figure 17 shows a truss structure with two elements and three nodes.
Figure 17: A truss structure with two elements and three nodes
In this case, the stiffness matrix will need to take into account the effect of element 1 on node 1 and node 2, and the effect of element 2 on node 2 and node 3.
Node 1
Node 2
The assembled stiffness matrix is then:
(14)
In 3D, the process becomes more complex, as now has components in the x,y and z. The coordinate system then needs to be rotated in terms of the global coordinate system. The transformation uses variables , and , and the length of the element () to form a transformation matrix for each element (.
0
0
0 0 0 0 0 0
(11) can then be rotated into the components.
(15) (16) (17) (18)
(19)
, 0
0
0 0
0
0 Page 30
Node 3
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
, ,
(20) (21)
(22)
Finally, once you have the stiffness matrix, displacement matrix and force matrix for the elements, these can be assembled as shown in 1D in (14) to form a huge matrix, including the effect of each element on each other. This is where MatLab will be helpful. For each element, insert the global coordinate stiffness matrix , into a stiffness matrix of the entire assembly. Remember that each node has 3 degrees of freedom listed in the matrices, so for example if element contains nodes and ,
,
0 0 0
,
,
: 2, : 2 ,1: 3,1: 3 : 2, : 2 ,4: 6,4: 6
(23) (24)
Some loads will have boundary conditions as well as loads. To apply boundary conditions, you remove degrees for freedom from the global stiffness matrix. You do this by identifying the rows and columns of the global stiffness matrix corresponding to displacement which are set to zero due to the boundary condition, and then removing these rows. Removing these rows and columns produces the reduced global stiffness matrix.
The reduced global matrix system can be solved for the displacements using .
(25)
Ensure that when you create the force and displacement matrices, you enforce the boundary conditions correctly, choosing known forces and removing elements where the displacement is already known.
Finally, Once the global displacement vector is calculated, it can be used to find strain and stress in each element. For each element connected to nodes and , convert back to local coordinates.
, ,@
(26) (27)
0 0 0
Then use the initial length to calculate strain .
Page 31
% Inputs:
% n_xyz:
coordinates of the nodes (m)
element number and corresponding nodes
densities of each element (kg/m^3)
cross-sectional areas of each element (m^2)
Young’s modulus of each element (Pa)
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
(28)
(29)
The following headers are for the internal subfunctions of RunFEA. You should not need to interact with these subfunctions.
function [M, K, Ts, Ls] = PreProcessFEA(n_xyz, e_n, rhos, As, Es, point_ms, Ne,
Nn)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PreProcessFEA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function “pre-processes” the FEA model in preparation for static or
modal FEA to be performed.
% Specifically it must:
% 1. Check the truss is stable using the stability condition m+6 >= 3j where m
is number of members (Ne) and j is number of joints (Nn).
% NOTE: There are cases where this is true but the truss is still not stable
e.g. if all nodes lie in a plane
% 2. Calculate the length of each element (Ls)
% 3. Calculate each element’s coordinate transformation matrix (Ts)
% 4. Assemble the Global Stiffness Matrix (K)
% 5. Assemble the Global Mass Matrix (M)
From the strain, the stress can be calculated from the Young’s Modulus using Hooke’s Law.
12.2 Interface for internal subfunctions
% e_n:
% rhos:
% As:
% Es:
% point_ms: Value of point masses on the each node (kg)
% Ne: number of elements
% Nn: number of nodes
% Outputs:
% M: Global Mass matrix
% K: Global Stiffness matrix
% Ts: Coordinate Transformation Matrix of each element (2*6*Ne matrix)
% Ls: Lengths of each element (m)
function [K_BC, M_BC, freeDOFs] = ApplyBCs(K, M, BCs, Nn)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ApplyBCs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function removes fixed degrees of freedom from input global mass and
% stiffness matricies, producing the trimmed global matricies.
% For input simplicity, it is implemented so that fixed nodes are fixed
% Inputs: % K:
% M: % BCs:
Global Stiffness Matrix
Global Mass Matrix
nBC*2 Matrix of boundary conditions with each fixed degree of
freedom on each row: [node number, fixed nodal DOF number], as
Page 32
% Nn:
% Outputs:
% K_BC:
Reduced Mass Matrix (DoFs with BCs removed)
Computational Engineering Analysis (AERO2463): MATLAB Project 2019
defined in the excel file for each load case. Nodal DOF number
gives direction in which it is fixed where x=1, y=2, z=3.
Number of nodes
Reduced Stiffness Matrix (DoFs with BCs removed)
% M_BC:
% freeDOFs: Matrix containing list free degrees of freedom (i.e.
row/column indexes of the global matricies without BCs
applied)
function [U, Strain, Stress] = StaticFEA(K_bc,M, Loads, a, freeDOFs, Nn, Ne,
e_n, Ts, Ls, Es)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% StaticFEA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Create reduced force vector (F_bc) from Loads and rigid body acceleration
(a)
% 2. Calculate global displacement from U_bc = K_bc \ F_bc
% 3. Calculate each element’s strain
% 4. Calculate each element’s stress from Strain * Es
% Inputs:
% K_bc:
% M:
% Loads:
Reduced Stiffness matrix
Global mass matrix (for converting acceleration to force on
each node)
nLoads*3 matrix containing following data for an applied load
on each row: [node number, DoF number, Load magnitude (N)]
Globally applied acceleration
% a:
% freeDOFs: List of all free (unfixed) Degrees of Freedom
% Nn:
% Ne:
% e_n:
% Ts:
% Ls:
% Es:
% Outputs:
% U:
Number of nodes
Number of elements
Element connectivity matrix
List of each element’s Coordinate Transformation Matrix
Length of each element (m)
Young’s Modulus of each element (Pa)
Global displacement vector
% Strain:
% Stress: Stress of each element
Strain of each element
Page 33