1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

[MATLAB] Modeling gravity using polar coordinates! Argh!

  1. Dec 12, 2014 #1
    I've been using MATLAB (ode45) to simulate the mechanics of a rocket under the forces of gravity, drag, and internal thrust.y I've recently refactored my simulation to include 2d space, orientation of the rocket, etc. (So I can try to make it orbit, finding optimal ascent profiles, etc)

    Because I want to retain the simplicity of my 1d model, I've attempted to use polar coordinates for defining the rockets position, velocity, and acceleration. I've implemented it, but when I run the simulation and plot height vs theta... I get a disturbing picture.

    Gravity does not affect the orbital velocity in any way. It always accelerates the height of the rocket downward. And yet for it to orbit, it needs to have a cyclic, non-zero height. Something is very wrong here, and I can't tell what it is.

    The line defining gravitational acceleration is as follows. I define it as a vector, with the first component being orbital acceleration (how much theta's rate of change changes due to gravity), and the second component being height's acceleration. (height is y, here)
    Code (Text):
    grav = [0; -S/(R+y)^2];
    Since S (the standard gravitation parameter) is positive, the acceleration of height will always be negative. The height of the rocket must inevitably go to zero. My intuition of gravity no longer makes sense.

    This is bothering me on a visceral level. I can't tell what I'm doing wrong. I now suspect the Earth will fall into Sol in the next few hours.
     
  2. jcsd
  3. Dec 12, 2014 #2

    jedishrfu

    Staff: Mentor

    It would help if you posted the whole MATLAB script as some folks have MATLAB and could run it to see whats happening.
     
  4. Dec 12, 2014 #3
    All right, I will post the whole thing. It's large though. There's four files, p_step.m, simul.m, eventy.m, and designer.m

    p_step.m is the function ode45 calls, and is basically a system of nonlinear differential equations.
    simul.m is the main script which calls ode45 and plots the data.
    eventy.m defines events (out_of_fuel, and ground_collision)
    designer.m specifies rocket properties (thrust, efficiency, etc)

    p_step.m:
    Code (Text):

    function dx = p_step(~, Z)

      global T II IF;
      global TH;
     

      S = 3.5316E12;
      R = 6E5;
      D = 9.784758844E-4;
      H = 5E3;
      G = 9.82;
     
      x = Z(1);
      y = Z(2);
      TH = y*3.14/200000;
     
      vx = Z(3);
      vy = Z(4);
      th = Z(5);
      m = Z(6);
     
      atm = exp(-y/H);
     
      grav = [0; -S/(R+y)^2];
      drag = -D*atm*[vx; vy].^2;
      thrust = T/m*[sin(TH); cos(TH)];
     
      a = thrust + drag + grav;
      ax = a(1);
      ay = a(2);
     
      dm = T/(((IF-II)*atm-IF)*G);
      dth = 0;
     
      dx = [vx; vy; ax; ay; dth; dm];
     
    end
     
    simul.m:
    Code (Text):


      global T II IF ME MF;

      type = 8;

      T = engine_thrusts(type);
      II = engine_atm_isps(type);
      IF = engine_vac_isps(type);
      ME = engine_masses(type);
     
      G = 9.82;
      R = 1.5;
     
      MI = T/(G*R);
      MF = (MI+7*ME)/8;
     
     
      opt = odeset('Events',@eventy);
     
      [t,y,~,ey,ei] = ode45(@p_step,[0 1000],[0, 0, 0, 0, 0, MI], opt);
      plot(y(:,1), y(:,2));
     
      if (numel(ei)~=0)
      switch ei(1)
      case 1
      max_dv = ey(3);
      case 2
      max_dv = -1;
      end
      else
      fprintf('Flight timed out\n');
      max_dv = 0;
      end
     
      fprintf('max-velocity:\t%.2f\n',max_dv);
     
    eventy.m:
    Code (Text):

    function [zero, terminal, negative] = eventy(~,x)

      global MF;
     
      zero = [x(6) - MF; sign(x(2))+1];
      terminal = [1; 1];
      negative = [0; 0];
     
    end
     
    designer.m:
    Code (Text):


    global engine_names engine_thrusts engine_masses engine_atm_isps engine_vac_isps;

    engine_names = ['LV-1 Liquid Fuel Engine'... % 1
      'Rockomax 48-7S'... % 2
      'LV-909 Liquid Fuel Engine'... % 3
      'LV-N Atomic Rocket Motor'... % 4
      'R.A.P.I.E.R. Engine'... % 5
      'Toroidal Aerospike Rocket'... % 6
      'LV-T45 Liquid Fuel Engine'... % 7
      'LV-T30 Liquid Fuel Engine'... % 8
      'Rockomax "Poodle" Liquid Engine'... % 9
      'Rockomax "Skipper" Liquid Engine'... % 10
      'Rockomax "Mainsail" Liquid Engine'... % 11
      'Kerbodyne KR-2L Advanced Engine'... % 12
      'S3 KS-25x4 Engine Cluster']; % 13

    engine_thrusts = [4 30 50 60 175 175 200 215 220 650 1500 2500 3200];

    engine_masses = [.03 .1 .5 2.25 1.2 1.5 1.5 1.25 2 3 6 6.5 9.75];

    engine_atm_isps = [220 300 300 220 320 388 320 320 270 320 320 280 320];

    engine_vac_isps = [290 350 390 800 360 390 370 370 390 370 360 380 360];

    fuel_names = ['Kerbodyne S3-14400 Tank'...
      'Kerbodyne S3-7200 Tank'...
      'Rockomax Jumbo-64 Fuel Tank'...
      'Kerbodyne S3-3600 Tank'...
      'Rockomax X200-32 Fuel Tank'...
      'Rockomax X200-16 Fuel Tank'...
      'FL-T800 Fuel Tank'...
      'Rockomax X200-8 Fuel Tank'...
      'FL-T400 Fuel Tank'...
      'FL-T200 Fuel Tank'...
      'FL-T100 Fuel Tank'...
      'ROUND-8 Toroidal Fuel Tank'...
      'Oscar-B Fuel Tank'];

    fuel_full_mass = [82 41 36 20.5 18 9 4.5 4.5 2.25 1.125 0.5625 0.136 0.078675];

    fuel_empty_mass = [10 5 4 2.5 2 1 0.5 0.5 0.25 0.125 0.0625 0.025 0.015];

     

    I don't think the problem is in the MATLAB code, though. I think it's a conceptual problem I'm having. If acceleration due to gravity is perpendicular (has no affect on) orbital velocity, why do orbiting objects get faster near the gravitational body?
     
  5. Dec 12, 2014 #4

    jedishrfu

    Staff: Mentor

    Thats large? Medium for programming is 10K lines and large is 100K lines...
     
  6. Dec 12, 2014 #5
    My program is plenty large! D:<

    Can you help me understand orbits jedish, or are you satisfied with insulting the length and girth of my program? ;)
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: [MATLAB] Modeling gravity using polar coordinates! Argh!
Loading...