Transcript Converting 2nd order ODEs to 2 1st order ODEs Example  MEAM 211

MEAM 211
Converting 2nd order ODEs to 2 1st order ODEs

Example
MEAM 211
Converting 2nd order ODEs to 2 1st order ODEs

Example

Define the state vector
MEAM 211
Converting 2nd order ODEs to 2 1st order ODEs

Example

Define the state vector

Write state equations
MEAM 211
Converting 2nd order ODEs to 2 1st order ODEs

Example

Define the state vector

Write state equations

With initial conditions
MEAM 211
Solving Equations of Motions for Particles
Example
 Particle in 3-D subject to thrust, gravity and drag

State vector is 6x1

6 State Equations
MEAM 211
MATLAB Example
A ball is thrown upward against gravitational attraction and air resistance
with an initial velocity of 30 meters/second. The air resistance opposes the
velocity and is proportional to the square of the velocity. The acceleration is:
a = - g – cv2 sign(v)
where g = 9.81 meter/sec2 and c = 0.001 1/meter. Solve for the position and
velocity of the particle as a function of time through a six second time interval.
main.m
timeInterval=[0, 6];
x0=[0; 30];
% interval for integration
% initial position = 0, velocity = 30
intFn ('vertical', timeInterval, x0);
function xdot=vertical(t, x);
c=0.001;
g=9.81;
% c=0.001 1/m
% g=9.81 m/sec/sec
xdot1 = x(2);
% derivative of position = velocity
xdot2 = -g - c*x(2)*abs(x(2));
% acceleration
vertical.m xdot=[xdot1;
xdot2];
MEAM 211
Example
Figure SA2.2.1 (p. 91)
Schematic of ejection seat test device.
The seat is moving at 600 knots when it is ejected with
a specified acceleration a-bt2 for 800 msec. A chute is
released to slow the chair down until the chair reaches
a speed of 100 knots. Then a parachute is deployed.
[Assume the parachute drag is 10 times the drag of the
first (drogue) chute.]
MEAM 211
Ejection Seat
T= 0: pilot pulls handle, T=0.15 secs: seat clears cockpit;
T=0.5-0.8 secs: seat/man separator fires; T = 200-400 seconds;
main parachute deploys
MEAM 211
if t < 0.8 % this is phase I
acceleration=a-b*t*t;
xdot(3) = -acceleration*sin(pi/9); % sin(20 deg) = sin
% This function computes the rate of change
(pi/9 rads)
% of position and velocity for the ejection seat
xdot(4) = acceleration*cos(pi/9); % cos(20 deg) = cos
(pi/9 rads)
%
else % phase II or III
% x(1) = position (y);
if speed > v_crit % this is phase II
% x(2) = velocity (ydot);
xdot(3)=-c*speed*xdot(1);
%
xdot(4)=-c*speed*xdot(2)-g;
else % we are in phase III
k1=1.1508122; % conversion from knots to miles per hour
if x(2) > 0 % the ejection seat is in the air
k2=1.4666249; % conversion from miles per hour to ft/sec.
xdot(3)=-10*c*speed*xdot(1);
k=k1*k2;
% conversion from knots to feet/second
xdot(4)=-10*c*speed*xdot(2)-g;
v_crit=k*100; % speed at which Phase III starts
else % the ejection seat has hit the ground
xdot(3)=0;
c=0.003;
% c=0.003 1/ft
xdot(4)=0;
g=32.2;
% g=32.2 ft/sec/sec
xdot(1)=0;
a=16*g;
% a=16g
xdot(2)=0;
b=9*g;
% b= 9g
end
end
end
xdot(1) = x(3); % derivative of x position = x velocity
xdot(2) = x(4); % derivative of y position = y velocity
xdot=xdot'; % need to transpose the vector to get
speed=sqrt(x(3)^2+x(4)^2);
% xdot to be a column vector
function xdot=ejectionSeat(t, x);