# Gravity turn solver only works for a gamma of 90 degrees

31 ビュー (過去 30 日間)
Jordan Booth 2023 年 2 月 1 日
コメント済み: William Rose 2023 年 2 月 7 日
I am trying to plot a gravtiy turn trajectory for a rocket given an inital pitch over angle using the following equations of motion: When the intial gamma is set to 90 there will be no pitch over and the solution seems correct so the v_dot equation should be correct. As soon as some angle less than 90 is used it doesnt seem to work.
main:
clc, clear
rE = 6371008.8; % radius of earth in m
t0 = 0;
tfinal = 1500; % total time of model
% starting values
p0 = [rE; 0; 0.001; 90]; % r, theta, velocity, gamma
[t,p] = ode45(@ode,[t0 tfinal],p0); % ode solver
fig1 = figure(1); % velocity plot
plot(t,p(:,3))
xlabel('time (s)')
ylabel('velocity (m/s)')
fig2 = figure(2); % altitude plot
plot(t,p(:,1)-rE) % height above surface of earth
xlabel('time (s)')
ylabel('height (m)')
axis([0 1500 0 3e6])
fig3 = figure(3); % ground relative plot
x = p(:,1).*cosd(p(:,2)); % polar to cartesian conversion
y = p(:,1).*sind(p(:,2));
plot(x,y)
hold on;
viscircles([0 0],rE);
axis([0 10e6 0 10e6])
axis square
hold off
ode function:
function dpdt = ode(t,p)
rE = 6371008.8; % radius of earth in m
dSL = 1.225; % sea level density in kg/m^3
h0 = 10400; % scale alt in m
A = 43.0084; % cross-sectional area in m^2
Cd = 0.237; % drag co-efficent
u = 3.986004418e14; % standard gravitational parameter
gA = 9.80665; % gravtiational acceleration at sea level in m/s^2
pm = 20000; % payload mass
% second stage
smp = 92670; % propellant mass in kg
sIsp = 348; % specific impulse in secs
sF = 981000; % thrust
sm0 = smp+3900+pm; % total mass
stf = (smp*gA*sIsp)/sF; % second stage burn time in secs
% first stage
fmp = 395700; % propellant mass in kg
fIsp = 283; % specific impulse in secs
fF = 7607000; % thrust
fm0 = fmp+25600+sm0; % total mass
ftf = (fmp*gA*fIsp)/fF; % first stage burn time in secs
aD = dSL*exp(-(p(1)-rE)/h0); % air density function
T = thrust(t,ftf,fF,stf,sF); % thrust function
m = mass(t,ftf,stf,sm0,smp,gA,sIsp,fm0,fIsp,T); % mass function
dpdt = [p(3)*sind(p(4)); % r
(p(3)*cosd(p(4)))/p(1); % theta
(-u*cosd(p(4)))/(p(3)*p(1)^2)+(p(3)*cosd(p(4)))/p(1)]; % gamma
end
outputs for inital gamma of 90 degrees:   outputs for initial gamma of 89 degrees:   The relative graph should show some pitch over if the function is working correctly
thrust:
function T = thrust(t,ftf,fF,stf,sF)
if t < ftf
T = fF;
elseif t < ftf + 12
T = 0;
elseif t < stf + 12 + ftf
T = sF;
else
T = 0;
end
end
mass:
function m = mass(t,ftf,stf,sm0,smp,gA,sIsp,fm0,fIsp,T)
if t >= stf + 12 + ftf
m = sm0-smp;
elseif t > ftf + 12
m = sm0-((T*(t-ftf-12))/(gA*sIsp));
elseif t > ftf
m = sm0;
else
m = fm0-((T*t)/(gA*fIsp));
end
end
The thrust and mass functions perform as expected:  I have been trying for a while to figure it out and have looked online for solutions but cant find anything similar that works. I would be grateful for any help, thanks :)

サインインしてコメントする。

### 採用された回答

William Rose 2023 年 2 月 2 日
This looks like a nice model!
Your equation for = d(theta)/dt = dp(2)/dt is
dpdt = [...
(p(3)*cosd(p(4)))/p(1); % theta...
I think you want to do
dpdt = [...
(180/pi)*p(3)*cosd(p(4))/p(1); % theta in degrees...
because you compute x(t),y(t) with cosd(p(:,2)) and sind(p(:,2)), i.e. you assume theta is in degrees.
Alternatively, you could compute x and y with cos(p(:,2)) and sin(p(:,2)), but then you would have theta in radians and gamma in degrees, which would be kind of wierd.
##### 17 件のコメント表示非表示 16 件の古いコメント
William Rose 2023 年 2 月 7 日
@Jordan Booth, well done! You're welcome.

サインインしてコメントする。

### カテゴリ

Find more on Geology in Help Center and File Exchange

R2021b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!