ode45 System of Equations Difficulties - Launch Vehicle Simulation

Question

JT 2018-6-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/406621-ode45-system-of-equations-difficulties-launch-vehicle-simulation

编辑： Jim Riggs 2018-6-26

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Jim Riggs 2018-6-20

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/406621-ode45-system-of-equations-difficulties-launch-vehicle-simulation#answer_325567

编辑：Jim Riggs 2018-6-20

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 个评论
显示 13更早的评论隐藏 13更早的评论

Jim Riggs 2018-6-25

编辑：Jim Riggs 2018-6-25

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.

JT 2018-6-26

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 2018-6-26

编辑：Jim Riggs 2018-6-26

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.

请先登录，再进行评论。