MAE 8 MATLAB Project Report

Man-Yeung Tsay
A11313074
MAE 8 MATLAB Project Report
Overview
This project models projectile motion of a ball on a 3-dimensional terrain. The objective was to
identify the location of where the ball would hit the terrain. Given an initial position (Xo, Yo,
and Zo in meters) and initial velocity (Umag in m/s) with launch directions (theta and phi in
degrees), the project computes the time (in seconds), 3 components of position (X, Y, Z in m),
and 3 components of velocity (U, V, W in m/s) along the trajectory. Theta is the angle from the
positive x-axis in the x-y plane. Phi is the angle from the positive z-axis to the vector Umag. In
all tasks, 9.81 m/s^2 is used as the value for gravity. Except for calculating the max height and
final position, the time step (dt) of 0.001s was used.
Task1
In the first task, air resistance was neglected, and projectile motion was modeled using Euler's
method. These equations were used to solve for projectile motion:
𝑑𝑋
=π‘ˆ
𝑑𝑑
π‘‘π‘ˆ
=0
𝑑𝑑
π‘‘π‘Œ
=𝑉
𝑑𝑑
𝑑𝑉
=0
𝑑𝑑
𝑑𝑍
=π‘Š
𝑑𝑑
Using Euler's method, the equations are:
𝑋(𝑛 + 1) = 𝑋(𝑛) + π‘ˆ(𝑛) βˆ— 𝑑𝑑
π‘Œ(𝑛 + 1) = π‘Œ(𝑛) + 𝑉(𝑛) βˆ— 𝑑𝑑
𝑍(𝑛 + 1) = 𝑍(𝑛) + π‘Š(𝑛) βˆ— 𝑑𝑑
π‘‘π‘Š
= βˆ’π‘”
𝑑𝑑
π‘ˆ(𝑛 + 1) = π‘ˆ(𝑛)
𝑉(𝑛 + 1) = 𝑉(𝑛)
π‘Š(𝑛 + 1) = π‘Š(𝑛) βˆ’ 𝑔 βˆ— 𝑑𝑑
A function was written in file "projectile3d.m" to calculate the equations to solve for the
positions and velocities along the trajectories of the ball. The inputs were the 3 components of
initial position (Xo, Yo, Zo), the initial velocity (Umag), and the launch directions (theta and
phi). A while loop was used to loop through the equations. The condition used for the while loop
was W(i)>=0. To find the max height position at where "W" would change signs, an "if"
statement was used. The condition for the "if" statement was "if W(i+1)<0 && W(i)>0," and the
outputs would be calculated with a new time step. To find where the ball lands on the ground, an
"if" statement of "if Z(i)<0" was used to calculate the outputs with a different time step. To find
the new time step for max height position and final position, the difference between the previous
Z value and current value was divided by the previous W value. To extract the max height values
from the outputs, the built-in command "max" was used to find the maximum Z value. A
variable was set to find the index for the maximum Z value. The remaining corresponding
outputs were called using this index.
Man-Yeung Tsay
A11313074
Task 2
In the second task, air resistance was included. The process for the second task mirrors the first
task, except for including air resistance in the equations. A separate file "projectiletask2.m" was
created. This new file mirrors "projectile3d.m", except for "C" (coefficient of friction) as an
additional input and the inclusion of air resistance in the equations. Four trajectories (with same
initial conditions of Xo=Yo=Zo=0, Umag=40m/s, theta=phi=45 degrees) were launched with
different coefficients of friction: C=0, 0.2, 0.4, and 0.6. The equations used are:
𝑑𝑋
π‘‘π‘ˆ
𝜌𝐴
=π‘ˆ
= βˆ’πΆ
π‘ˆβˆšπ‘ˆ 2 + 𝑉 2 + π‘Š 2
𝑑𝑑
𝑑𝑑
2π‘š
π‘‘π‘Œ
𝑑𝑉
𝜌𝐴
=𝑉
= βˆ’πΆ
𝑉 βˆšπ‘ˆ 2 + 𝑉 2 + π‘Š 2
𝑑𝑑
𝑑𝑑
2π‘š
𝑑𝑍
π‘‘π‘Š
𝜌𝐴
=π‘Š
= βˆ’π‘” βˆ’ 𝐢
π‘Š βˆšπ‘ˆ 2 + 𝑉 2 + π‘Š 2
𝑑𝑑
𝑑𝑑
2π‘š
C=coefficient of friction; A=Ο€ r2 where r=0.04 m; m=0.15kg, ρ=1.2kg/m3 (air)
Using Euler's method, the equations are:
𝑋(𝑛 + 1) = 𝑋(𝑛) + π‘ˆ(𝑛) βˆ— 𝑑𝑑
π‘Œ(𝑛 + 1) = π‘Œ(𝑛) + 𝑉(𝑛) βˆ— 𝑑𝑑
𝑍(𝑛 + 1) = 𝑍(𝑛) + π‘Š(𝑛) βˆ— 𝑑𝑑
𝜌𝐴 √ 2
π‘ˆ π‘ˆ + 𝑉2 + π‘Š2 ) βˆ— 𝑑𝑑
2π‘š
𝜌𝐴 √ 2
𝑉(𝑛 + 1) = 𝑉(𝑛) βˆ’ (𝐢
𝑉 π‘ˆ + 𝑉2 + π‘Š2 ) βˆ— 𝑑𝑑
2π‘š
𝜌𝐴 √ 2
π‘Š(𝑛 + 1) = π‘Š(𝑛) + (βˆ’π‘” βˆ’ 𝐢
π‘Š π‘ˆ + 𝑉2 + π‘Š2 ) βˆ— 𝑑𝑑
2π‘š
π‘ˆ(𝑛 + 1) = π‘ˆ(𝑛) βˆ’ (𝐢
Task 3
The goal of the third to task was to accomplish the main objective of this project. A new file
"trackprojectile.m" was created to incorporate the terrain and locate the final positions of the ball on the
terrain. This new file mirrors the file from the second task "projectiletask2.m" with some exceptions. The
built-in MATLAB command "interp2" was used for interpolation. The condition for the while loop was
"Z(i) >= interp2(terrain.x, terrain.y, terrain.z, X(i), Y(i)." Then, a variable "surface" was created by
setting it equal to the given value from "terrain.mat" after the while loop. To find the final position of the
ball, an "if" statement of "if Z(i)<surface" was created, and a new time step was used to calculate the
outputs for the final position. Similarly to the previous tasks, this new time step was calculated by
finding the difference between the previous Z value and current value and dividing it by the
previous W value.
Man-Yeung Tsay
A11313074
Figures and Tables
Task 1
Initial Conditions for Task 1
Results for Task 1
Man-Yeung Tsay
A11313074
Task 2 - Initial conditions of Xo=Yo=Zo=0, Umag=40m/s, theta=phi=45 degrees
with C=0, 0.2, 0.4, 0.6
Results for Task 2
Man-Yeung Tsay
A11313074
Task 3
Initial Conditions for Task 3
Results for Task 3
Man-Yeung Tsay
A11313074
Flow Chart of Programs
Start with Initial Inputs
(Xo,Yo,Zo,Umag,theta,phi, C (if there is air resistance))
projectile3d.m (for no air resistance)
projectiletask2.m (with air resistance)
calculates equations of motion
using Euler's method until ball has
landed
when ball has landed, a new time
step is calculated then used to find
the final outputs
for max height, a new time step is
calculated when W changes sign,
and then this time step is used to
find the position, velocities, and
time during max height
End with Outputs
(T, X, Y, Z, U, V, W)
Man-Yeung Tsay
A11313074
Start with Initial Inputs
(Xo,Yo,Zo,Umag,theta,phi, C)
trackprojectile.m
calculate until ball has landed
(using interp2 for interpolation)
find where ball has landed
find new time step to calculate
final outputs
End with Outputs
(T, X, Y, Z, U, V, W)
Man-Yeung Tsay
A11313074
Appendix
project.m - runs all the commands for tasks
projectile3d.m - solves for positions and velocities along trajectories of ball with no air
resistance
Inputs
Xo, Yo, Zo (in meters) are components of initial position
Umag (in meters/sec) is initial velocity
theta and phi (in degrees) are components of launch direction
Outputs
X, Y, Z (in meters) are components of position
U, V, W (in meters/sec) are components of velocity
T=time
projectiletask2.m -a modified version of projectile3d.m that solves for positions and velocities
along trajectories of ball with air resistance
Inputs
Xo, Yo, Zo (in meters) are components of initial position
Umag (in meters/sec) is initial velocity
theta and phi (in degrees) are components of launch direction
C is coefficient of friction
Outputs
X, Y, Z (in meters) are components of position
U, V, W (in meters/sec) are components of velocity
T=time
trackprojectile.m - tracks motion of ball as it lands on terrain and calculates positions and
velocities of trajectory
Inputs
Xo, Yo, Zo (in meters) are components of initial position
Umag (in meters/sec) is initial velocity
theta and phi (in degrees) are components of launch direction
C is coefficient of friction
Outputs
X, Y, Z (in meters) are components of position
U, V, W (in meters/sec) are components of velocity
T=time