Gravity Gradient Torque

[Arrowstar]

Hi guys!

I'm working on a little project to try and simulate gravity gradient torque for a generic spacecraft. I'm doing my work right now in MATLAB to try and develop the mathematics needed to get all of this right. However, I've run into a little snag and I was hoping someone can show me where I might be going wrong.

If we assume the vehicle is a rigid body and I set up a body-fixed coordinate system along principle axis, then I have Euler's Equations of Motion:

M_x = I_xx * omega_x_dot + (Izz - Iyy)*omega_y*omega_z
M_y = I_yy * omega_y_dot + (Ixx - Izz)*omega_z*omega_x
M_z = I_zz * omega_z_dot + (Iyy - Ixx)*omega_x*omega_y

...where "M_n" are the torques around x, y, z axis, "Inn" are principle moments of inertia, and "omega_n" are angular rates about x, y, and z axis. Note that I'm using lowercase x, y, and z for the body axis.

Now, if I am just looking at gravity gradient torques, then I can describe the left-hand side of Euler's Equations by the following (which I take from Curtis' "Orbital Mechanics for Engineering Students"):

M_x_g=((3*muEarth*Ry*Rz)/R^5) * (Izz - Iyy)
M_y_g=((3*muEarth*Rx*Rz)/R^5) * (Ixx - Izz)
M_z_g=((3*muEarth*Rx*Ry)/R^5) * (Iyy - Ixx)

Here, "muEarth" is the gravitational parameter of the Earth (but it could be whatever body I'm interested in), "R" is the radial distance from the CoM of the Earth to the CoM of the spacecraft, "Rn" are the x, y, z components of that R, and "Inn" is as above.

So, when I take these moment equations and plug them into Euler's equations, I get first order differential equations I can numerically integrate over time. Great. (From here I can go on to the kinematic equations and get the actual vehicle attitude, but bear with me.) However, I would expect that in a system where I only have gravity gradient torque, the torques would damp out eventually as the vehicle settled into an equilibrium position. That's not what I'm seeing, however: my calculations show that the angular rates of my vehicle actually grow over the course of a week or so! This can't be right.

Let me pose two questions:

1) I am not taking orbital motion into account, and thus I am not taking into account the attitude changes one would get as Earth "falls under" the spacecraft while in orbit. Is this important?

2) The distance vectors I'm using for the M_n_g gravity gradient torques are body frame vectors. Is this correct?

I can post some output if anyone would like. Can someone help me?

Thanks!

---------- Post added at 01:27 PM ---------- Previous post was at 11:34 AM ----------

Say, for example, I use the following initial conditions:

omega_x = omega_y = omega_z = 0
phi_x=pi/6
phi_y = phi_z = 0

...where phi_n are the euler angles.

In particular, here's what I see when I use an RK4 integration scheme on the above equations (with a 312 Euler angle kinematic equation thrown in):

I am perplexed, to say the least.

 

[Wishbone]

Aren't you better off by using implicit methods? Attitude equations are considered stiff.

A qualification: my diff.eq. knowledge mostly comes from reading books, is far from satisfactory, and some folks around here would consider self-educated people "half-wits".

 

[Arrowstar]

I did actually experiment with MATLAB's stiff solvers, but none of them really give me any better results. I did also think about trying the implicit solver, but that's a bit more work, so I haven't implemented anything like that yet.

By your comment, are you suggesting that my problem might be numerical and not in the solution of the equations themselves?

 

[tblaxland]

However, I would expect that in a system where I only have gravity gradient torque, the torques would damp out eventually as the vehicle settled into an equilibrium position.

1) I am not taking orbital motion into account, and thus I am not taking into account the attitude changes one would get as Earth "falls under" the spacecraft while in orbit. Is this important?

If you do not take the orbital motion into account, then your spacecraft will just behave like a pendulum - the angular rate will vary periodically, which looks like what you are getting. The reason that the angular motion dampens out in orbit is because, as the orbiting body moves around the parent, the direction of the gravity gradient changes in the inertial frame. For your current model, the direction of the gravity gradient in the inertial frame is constant and, in a frame of reference rotating with the rigid body axes, the direction of the gravity gradient does not change any faster or slower than the angular rate of the orbiting body. This means that there is no mechanism for exchanging the angular momentum of the orbiting body with its orbital angular momentum, and therefore no damping mechanism.

2) The distance vectors I'm using for the M_n_g gravity gradient torques are body frame vectors. Is this correct?

That looks correct to me.

 

[Arrowstar]

If you do not take the orbital motion into account, then your spacecraft will just behave like a pendulum - the angular rate will vary periodically, which looks like what you are getting. The reason that the angular motion dampens out in orbit is because, as the orbiting body moves around the parent, the direction of the gravity gradient changes in the inertial frame. For your current model, the direction of the gravity gradient in the inertial frame is constant and, in a frame of reference rotating with the rigid body axes, the direction of the gravity gradient does not change any faster or slower than the angular rate of the orbiting body. This means that there is no mechanism for exchanging the angular momentum of the orbiting body with its orbital angular momentum, and therefore no damping mechanism.

Hmm, okay. That does make physical sense, which is why I wondered about it when I started putting this together. Figuring out the angular rate I need to add to the vehicle to get "real life" behavior is going to be a bit of trick I think. I'll have to play around with some sketches and see if I can't decipher that.

Thanks for the help, tblaxland! :thumbup:

 

[tblaxland]

Figuring out the angular rate I need to add to the vehicle to get "real life" behavior is going to be a bit of trick I think.

I think you just need to take the vehicle's angular velocity vector and subtract the orbital angular velocity vector. Thinking that through, if the vehicle was tidally locked, both angular velocity vectors would be the same resulting in a zero net angular velocity of the parent body w.r.t. to the vehicle.

 

[Arrowstar]

I think you just need to take the vehicle's angular velocity vector and subtract the orbital angular velocity vector. Thinking that through, if the vehicle was tidally locked, both angular velocity vectors would be the same resulting in a zero net angular velocity of the parent body w.r.t. to the vehicle.

Well, that's looks like it's... kinda close. The results I get are a bit more stable, but still not settling down at the correct equilibrium position. Instead, the omega vector appears to settle down to a range of values and the Euler angle that I'm interested in (call it pitch, for lack of a better description) goes off to negative infinity (basically it's tumbling around that axis).

 

[tblaxland]

How does the magnitude of omega_vehicle - omega_orbit trend over time?

 

[n72.75]

Dormand–Prince is better then Runga-Kutta.

 

[tblaxland]

Dormand–Prince is better then Runga-Kutta.

I'm no expert but:

[Dormand–Prince] method is a member of the Runge–Kutta family of ODE solvers

So "better" doesn't really work there. Perhaps you mean it is better than other members of the R-K family? Also, I'm curious, better for what? Speed, accuracy, easier to code... :shrug:

 

[JEL]

@Arrowstar; Can this be of any help to you? :

http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19660013620_1966013620.pdf

http://www.eageseg.org/data/egm2010/Sessione B/Poster Papers/B_PP_02.pdf