ode45 System of Equations Difficulties - Launch Vehicle Simulation
11 views (last 30 days)
Show older comments
I am attempting to solve a system of ODEs that describe a shuttle launch. The calculation incorporates time-dependent unit step functions with differential equations to describe launch kinematics with realistic engine burn and thrust vectoring profiles.
I am using this document for reference, trying to recreate the results presented on pages 24 - 25 ('Requirement 15') of the PDF which were originally computed in Mathematica using the NDSolve[] functionality. I have successfully recreated the results from Part I and Part II of the document but have encountered a lot of difficulty in Part III. For this system of equations, the output from 'ode45' will output the following warning if a 'NonNegative' argument is not included in 'odeset':
'Warning: Failure at t=5.106812e+01. Unable to meet integration tolerances without reducing the step size below the smallest value
allowed (1.136868e-13) at time t.'
If a 'NonNegative' argument is provided, the solution covers the given 'tspan' but is incorrect according to the reference document. The 4th and 5th index output values have the correct shape over time but incorrect values. I am to the point where I cannot determine the origin of the errors; I believe the system of equations and the time-dependent unit step functions are defined correctly within the context of 'ode45' but I cannot justify the rationale behind the use of the 'NonNegative' specifier. When plotted using a time vector the non-DE unit step functions ('SRBT' and 'MET') look correct compared to the reference PDF. My code is as follows:
%%Requirement 15
R = 20.9e6; %earth radius in ft
g = 32.174; %gravitational acceleration in ft/s^2
tspan = [0 720];
y0 = [R 0 0 1530 4.5e6/g];
opts = odeset('NonNegative',[2 3 4 5]);
[T,Y] = ode45(@(t,y) myode(t, y), tspan, y0, opts);
function dydt = myode(t, y)
R = 20.9e6; %earth radius in ft
g = 32.174; %gravitational acceleration in ft/s^2
% Define Heaviside step function
hvsd = @(x) [1*(x == 0) + (x > 0)];
%Fuel rate - Main Engines
frMET = 117.9; %slugs/s
%Fuel rate - Solid Rocket Boosters
frSRB = 94.7; %slugs/s
% Mass of SRBs
mSRB = 2*6250; %slugs
% External tank mass
mEXT = 2062; %slugs
% SRB Thrust Schedule
SRBT = ((6.6 - 2.2.*(t/50)).*(hvsd(t) - hvsd(t-50)) +...
4.4.*(hvsd(t-50) - hvsd(t-120)));
% ME Thrust Schedule
MET = (1.125.*(hvsd(t) - hvsd(t-26) + 0.95.*(hvsd(t-26) - hvsd(t-60))) +...
1.125.*(hvsd(t-60) - hvsd(t-460)) + 0.731.*(hvsd(t-460) - hvsd(t-480)));
% Thrust Vector Control
alpha = 90.*(hvsd(t) - hvsd(t - 6)) +...
(90-(12.*((t - 6)./20))).*(hvsd(t - 6)- hvsd(t - 26))+...
(78-(57.*((t - 26)./214))).*(hvsd(t - 26)- hvsd(t - 240)) + atand(y(2)/y(4))*hvsd(t - 240);
% Drag Coefficient
D = 430*(2.33e-3)*exp(-(y(1) - R)/28276)*(y(2)^2 + (y(4) - 1530)^2);
dy1 = y(2);
dy2 = y(4)^2/y(1) - g*(R/y(1))^2 + (1/y(5))*((SRBT + MET)*(1.25 - 0.25*exp(-(y(1) - R)/28276)) - D)*sind(alpha);
dy3 = y(4)/y(1);
dy4 = (y(2)*y(4))/y(1) + (1/y(5))*((SRBT + MET)*(1.25 - 0.25*exp(-(y(1) - R)/28276)) - D)*cosd(alpha);
dy5 = -frMET.*MET - frSRB.*SRBT - mSRB.*(hvsd(t - 119.5) - hvsd(t - 120.5)) - ...
mEXT.*(hvsd(t - 479.5) - hvsd(t - 480.5));
dydt = [dy1 dy2 dy3 dy4 dy5]';
end
Any input regarding the mechanics of the 'NonNegative' specifier or correct implementation of time-dependent unit step functions would be appreciated. My objective is to achieve matching results with the PDF document.
Thanks!
0 Comments
Accepted Answer
Jim Riggs
on 20 Jun 2018
Edited: Jim Riggs
on 20 Jun 2018
The ODE solver has a hard time dealing with discontinuities. If any part of the equation set contains discontinuities, or even very rapid slope changes, you will need to break up the solution into segments based on the location of the discontinuities.
I see from the pdf paper that there is a discontinuity in the thrust profile at 50 seconds. This could be the cause of the error you are getting at 51 seconds.
Another option is to use a different solver with less stringent error tolerances. ODE45 is about the most stringent.
You might want to try ode23s or ode15s which are designed for stiff systems.
15 Comments
Jim Riggs
on 26 Jun 2018
Edited: Jim Riggs
on 26 Jun 2018
My pleasure. This is the kind of problem that I enjoy.
By the way, normally the thrust of a rocket is computed by specifying the specific impule, Isp, for the motor, then the thrust is computed from
Thrust = Isp * mdot ( mdot being the mass flow rate in pounds-mass per sec.)
an alternate form is
Thrust = Isp * Mdot * gravity (where mdot is in slugs/sec)
Isp has units of seconds. The total mass of fuel is specified (as an initial condition) and is subtracted from the system mass as the fuel burns. This way, there is a specified amout of fuel mass which is burned and the burn rate (mdot) determines how long it will burn. When the mass of fuel is used up, no more thrust. The Isp is the measure of efficiency of the motor.
The way that the thrust and mass are related in this paper, you have varying amounts of fuel mass burned based on the 'fr' factor. This results in a change in the system mass if/when the thrust level changes. The thrust curve is specified and the "fr" factor determines how much mass is used to produce the thrust. Not physically correct.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!