PDA

View Full Version : Modelling a rocket

Glom
2006-Feb-12, 07:28 PM
I'm writing a program to model the launch of a rocket. Here's the particular function so far.

//1. Load the position and speed of the rocket on the pad
double lslat = ( latdegbox->value() + ( latminbox->value() + latsecbox->value()/60) / 60 ) * M_PI / 180;
if ( southbutton->set() ) {
lslat *= -1;
}
double lslong = ( longdegbox->value() + ( longminbox->value() + longsecbox->value()/60) / 60 ) * M_PI / 180;
if ( westbutton->set() ) {
lslong *= -1;
}
double lselev = elevbox->value() * 0.3048;
double t0 = ( utchourbox->value() + ( utcminbox->value() + utcsecbox->value()/60) / 60 ) * M_PI / 180;
double x = (lselev + georad) * cos (lslat) * cos ( lslong + t0);
double y = (lselev + georad) * cos (lslat) * sin ( lslong + t0);
double z = (lselev + georad) * sin (lslat);
vec r(x,y,z);
x = - ( lselev + georad ) * cos (lslat) * georate * sin (lslong);
y = ( lselev + georad ) * cos (lslat) * georate * cos (lslong);
z = 0;
vec v(x,y,z);
lsspeedbox->value(v.mag());
vec att = r.unit();
double t=0;
double tstep = 0.001;

//2. Load the parameters of the launch
double oew = oewbox->value();
double zfw = oew + payload;
double fuel = fuelbox->value();
double mass = zfw + fuel;
double mdot = mdotbox->value();
double vex = vexbox->value();
double azimuth = azimuthbox->value() * M_PI / 180;
double pitchtime = pitchtimebox->value() * 0.3048;
double pitchrate = pitchratebox->value();
double pitchterm = pitchtermbox->value();
double pitch = 0;
vec g = r*(-mu/pow(r.mag(),3));
vec mg = g*mass;
towbox->value(mg.mag());
thrustbox->value(mdot*vex);

do {
//mdv = cdm + mgdt
vec g = r*(-mu/pow(r.mag(),3));
v = v + att*vex*mdot*tstep/mass + g*tstep;
vec gforce = att*vex*mdot/mass;
r = r + v*tstep;
mass -= mdot*tstep;
t += tstep;
wtf << t << " " << r.mag() - georad << " " << v.mag() << " " << gforce.mag()/9.81 << endl;
}
while ( r.mag() - georad < pitchtime && mass > zfw);

The problem bit is that upon reaching the altitude of pitchtime, the rocket is then supposed to pitch over at an inputted rate toward a final pitch at an inputted azimuth. The describe the attitude of the rocket, I created the vec att, which is a unit vector pointing in the direction of the rocket (vec is a class defined to act like a 3D physics vector). The tricky part is that I'm not sure how to alter the value of the vec att through the pitch rate and the azimuth.

Launch window
2006-Feb-14, 09:38 AM
sorry, I can't help

Nicolas
2006-Feb-14, 10:57 AM
you have a timestep "tstep" and a pitch rate "prate". That gives an angle "deltapitch" = prate*tstep per time interval. You turn the vector att around it's local pitch axis with this angle. To find the components of the vector, keep track of the roll, pitch and yaw angles:

pitch = pitch + detapitch

Now apply Euler angles to get the parameters of vec:

Euler angles (http://mathworld.wolfram.com/EulerAngles.html)

You can calculate the matrices B,C and D. From that, you can calculate matrix A.

vec = A.vec gives the new vector. (vec = [x y z]' in my case)

Unless I'm mistaken here, I'm writing this on the fly :).

I explained things for pitch; keep track of roll and yaw in the same way and take care of order of angles and matices.

Glom
2006-Feb-14, 05:14 PM
Thanks, Nicolas. I'll have to read through that multiple times in order to grasp it.

For the moment though, I cheated a little and eliminated azimuth control. The rocket can only be launched easterly for the moment. That was a simplification that made the maths easier. What I did was I worked out a local North unit vector. The North vector would lie in the plane of r and K, but would be perpendicular to r. The solution was obviously cross products. N = r × ( K × r ). Then I apply my trusty unit vector operator to give it unit magnitude.

The next step was to get the delta att vector. The simplest way to find it would be to imagine that delta att was perpendicular to att. That way, the magnitude would be tan(delta pitch). If I limited the direction to easterly (perpendicular to N), then the direction would be given by N × att. Since those are both unit vectors, the magnitude is unity before scaling it with the tan(delta pitch). Then I add delta att to att to get the new att, which has the right direction, but unfortunately, not unit magnitude, so I simply apply the unit operator.

Easy!

Here's the code I have for the burn of the first stage (the rocket is now three stage).

do {
//mdv = cdm + mgdt
vec g = r*(-mu/pow(r.mag(),3));
vec N = (r.cross( K.cross(r) )).unit();
string phase = "Straight ascent";
if ( (r.mag() - georad) > pitchtime && att.arg(r) < pitchterm ) {
att = (att + N.cross(att)*tan(pitchrate*tstep)).unit();
phase = "Pitching";
}
v = v + att*vex1*mdot1*tstep/mass1 + g*tstep;
vec gforce = att*vex1*mdot1/mass1;
r = r + v*tstep;
mass1 -= mdot1*tstep;
t += tstep;
telapsebox->value(t);