ode45 System of Equations Difficulties - Launch Vehicle Simulation

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!

Answer 1

[Jim Riggs]

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.

[JT]

Thank you for the feedback on this, Jim. I did suspect potential issues could lie with the step functions, particularly around transitions. I will attempt the solution over discrete intervals bounded by the linear regions of the time-dependent profiles and provide an update.

[Jim Riggs]

The simplest approach would be to try the solver which is better suited for stiff systems, such as ode23s. This would not require you to break up your solution into multiple parts. Then, if the stiff solver works you are done. If the stiff solver doesn't work, there might be another issue.

[JT]

Neither ode23s or ode15s yielded improved results. Even on a mitigated time scale (tspan = [0 50]) the solution still appears to be approaching a singularity around 51s or so. ode15s allows for the 'NonNegative' specifier which yields function shapes for the 4th and 5th outputs that look correct but are not to scale - all other outputs are clearly incorrect.

[Jim Riggs]

I ran your code and I get the same result. With the ode45 solver I get the same error at 51 seconds, and I get a solution using "Nonnegative" but it doesn't make sense. There is no positive acceleration in the radial direction, so R remains constant and Vr is zero.

One thing I would do is check the application of alpha to make sure that it is being properly applied using units of degrees.

[Jim Riggs]

OK. The thrust terms are in millions of pounds. They need to be in pounds, so a factor of 1e6 is needed.

[JT]

To address the two points above, I did check alpha through several time steps before 50s in the debugger and it appears correct, declining in value from 90 degrees after 6s.

The units of the thrust profiles are a bit of a mystery. The document does not provide the 1e6 scale factor but earlier, I did think to add it. Doing so however, makes it so that the mass plot is no longer correct (output: y(5)) and the computation crashes with the following message:

 Warning: Failure at t=1.846013e-04.  Unable to meet integration 
 tolerances without reducing the step size below the smallest value
 allowed (4.336809e-19) at time t.

This makes me wonder if the thrust profile was left intentionally un-scaled in the document. I received the same crash and warning message between ode45 and ode15s solvers.

[Jim Riggs]

Yes, I am seeing the same behavior. But this does solve the problem of the initial Rdot is no longer negative.

[Jim Riggs]

Here is what I am finding;

There is a disparity between the thrust and the mass flow rate. As written, the thrust is too low to push the body and there is no radial acceleration. Increasing the thrust by a factor of 1e6 produces a radial acceleration that seems appropriate, but then the mass flow rate is way too high. So, I added a factor of (1/1e6) to the mass rate of change (dydt5) (just on the two thrust terms). This looks better to me. See what you think.

[JT]

With this correction, there is now pretty good agreement prior to ~480s. When the solution is calculated for the full 720s, wild deviations occur. Let me know if you're not seeing this, just want to ensure I'm not missing anything else. Using ode15s yields a more accurate mass profile to the document.

Two other things to note, the radial displacement value (y(1)) in the document is decreasing from the IC (hence h = IC-y(1) for the plot), instead, in our solution, it is increasing up until ~320s. The other point is that the plot for theta appears to start below 0, indicating a negative angular value for the IC. The document specifies 0 degrees for the IC however.

Thank you for all your time and hard work on this, Jim.

[JT]

Okay - by setting the non-negative specifier on argument 2 using ode15s, I am able to get results that are fairly consistent to the reference document. This parameter will prevent negative radial velocities and force stability across the time scale. The theta output argument (y(3)) will need to be converted to degrees for plotting, as well.

There are still scaling errors and dissimilarities but this is much closer than prior efforts. Let me know if you have a chance to implement this.

[Jim Riggs]

OK. I have one additional finding. I am looking at the ode45 solver with the "NonNegative" constraint. This solver sometimes is happy taking very large time steps (~10 seconds). When this happens, It can completely miss the SRB mass ejection transient, and the main fuel tank ejection transient. These are both very large mass ejections which occur in narrow, one second wide time windows, ( T=119.5 to 120.5, and t=479.5 to 480.5). If this happens, it results in the system mass being too large, consequently, the accelerations are too small, and the velocities are reduced. This can be mitigated by adding the "MaxStep" constraint to the solver options. This is something to be aware of.

If you find my input helpful, please accept my answer.

[Jim Riggs]

Here are some more observations: I continue to get a lower peak Vr velocity than is shown in the pdf paper. I get a peak of 3452 at t=120) and the pdf paper shows a higher peak, closer to 4000.

I am getting a final Radial position of 112.2 miles (this increases to 133.4 if I zero the drag). The pdf paper shows 177.4 miles. So my results show significantly lower velocity and altitude than the paper.

Now, note that the plot for the mass vs. time in the pdf paper (page 23) does not show the drop (vertical discontinuity) at t=480 where the main fuel tank is dropped. In fact, it shows a corner where the mass stays constant, but the corner is significantly past the 480 second (8 minute) mark. In addition, I get a corner of 53,770 at t=120, where the pdf paper shows a value closer to 40,000 at the 2 minute mark. The differences in this figure suggest that there are some internal inconsistencies in the paper, or that the graphs are not properly scaled, and maybe you will never get an exact match.

[Jim Riggs]

I found another significant difference.

In your original code, above, you have the following constants:

frMET = 117.9;
frSRB = 94.7;

However, the paper gives values (on page 10) of

ME rate = 94.7 slugs/s per Mlbs thrust

SRB rate = 105.4 slugs/s per Mlbs thrust

This makes more sense, since the main engines use liquid fuel and should be more efficient than the solid rockets. By making this change, I am now getting a better match with the Vr curve, and the radial position is increased up to 124 miles (147 with no drag). The mass vs. time curve looks closer to the paper as well. This change is a significant improvement.

[JT]

This is a good catch; the paper mixes up these rates and drops the 105.4 slugs/s rate altogether in the formula for mass rate, replacing it by 117.9 slugs/s. There are several unexplained inconsistencies throughout the document like this.

I took liberties to adjust the SRB rate and can get good agreement by setting it slightly higher. At a certain point, if this rate is too high (117.9 slugs/s for example), the radial velocity will not settle to 0 ft/s and the solution will diverge elsewhere. It is also sensitive to the solver step size as well; enforcing too much resolution through the 'MaxStep' parameter will prevent the radial velocity from settling. Despite solver behavior and document discrepancies, I am now happy with the result. Thank you for your hard work on this, Jim.

[Jim Riggs]

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.