Orbiter-Forum  

Go Back   Orbiter-Forum > Far Side of the Moon > Math & Physics
Register Blogs Orbinauts List Social Groups FAQ Projects Mark Forums Read

Math & Physics Mathematical and physical problems of space flight and astronomy.

Reply
 
Thread Tools
Old 05-15-2018, 03:44 PM   #1
Raccoon
Orbinaut
Default Yet another propagation problem (atmospheric drag mostly)

Hello! I’ve write a program for orbit propagation and choose numerical integration method runge-kutta 4. I wrote only 2 perturbation effect: 2-body and atmospheric drag. Here is the code of main cycle:

Code:
Satellite sat(QVector3D(7000,0,500), QVector3D(0,9,0),0.0025, 3725.0);
    QVector< QVector3D > satellite;
    double h = 1.0;
    QDateTime dt = QDateTime(QDate(1983,9,23),QTime(0,0,0,0));
    for (int j = 0.0; j < 100; j++) {
        double Mjd_UTC = MJD(dt);
        // 32.184 – UTC-TAI
        double Mjd_TT  = Mjd_UTC + 32.184/86400.0;
        RK4step(h, sat, Mjd_TT);
        dt = dt.addDays(j);
        satellite.push_back(sat.loc);
        satellite.push_back(sat.vel);
    }
My implementation of RK4:

Code:
void RK4step(const double& h, Satellite& y, const double &Mjd_TT) {
    QVector3D loc(y.loc);
    QVector3D vel(y.vel);
    QVector3D k1,k2,k3,k4, acceleration;

    // rotation matrix for satellite relative to Earth at time Mjd_TT
    // O.Montenbruck, E. Gill; Satellite Orbits - Models, Methods, and Applications
    QMatrix3x3 T = NutationMatrix(Mjd_TT) * PreccessionMatrix(MJD_J2000, Mjd_TT);

    k1 = partialStep(y,            T, y.size, y.mass);
    k2 = partialStep(y + k1 * h/2, T, y.size, y.mass);
    k3 = partialStep(y + k2 * h/2, T, y.size, y.mass);
    k4 = partialStep(y + k3 * h,   T, y.size, y.mass);

    acceleration += (k1 + (k2 + k3) * 2 + k4) / 6;

    y.vel += acceleration * h;
    y.loc += y.vel * h;
}
Partial calculation:

Code:
QVector3D partialStep(const Satellite &y, const QMatrix3x3 &T,
                      const double& satSize, const double& satMass ) {
QVector3D acceleration (0.0, 0.0, 0.0);
double r = Magnitude(y.loc);
// MU - Earth gravity = 398602 
    double temp = -MU / (r * r * r);
    QVector3D acc2Body = y.loc * temp;

    double Cd = 2.2;
    QVector3D accFromAir = AccelDrag(0.0, y.loc, y.vel, T, satSize, satMass, Cd);

    acceleration = acc2Body + accFromAir;

    return acceleration;
}
Implementation of atmospheric drag which I find in “O.Montenbruck, E. Gill; Satellite Orbits - Models, Methods, and Applications” :

Code:
QVector3D AccelDrag ( double Mjd_TT, const QVector3D &loc, const QVector3D &vel,
                   const QMatrix3x3& T, double Area, double mass, double CD )
{
    // W3 – Earth rotation = 0.0000729211 rad/sec
    const QVector3D omega (0.0, 0.0, W3);
    double vel_abs, dens;
    QVector3D loc_tod, vel_tod;
    QVector3D vel_rel, acc_tod;
    QMatrix3x3 T_trp;
    // for converting to ICRF/EME2000
    T_trp = T.transposed();
    
    loc_tod = T * loc;
    vel_tod = T * vel;
    // velocity relative to Earth rotation
    vel_rel = vel_tod - QVector3D::crossProduct(omega, loc_tod);
    vel_abs = Magnitude2(vel_rel);
    // Harris-Priester atmospheric density model
    dens = Density_HP (Mjd_TT, loc_tod);
    // atmospheric acceleration
    acc_tod = -0.5*CD*(Area/mass)*dens*vel_abs*vel_rel;
    return T_trp * acc_tod;
}
With solo 2-body perturbation effect orbit change correctly to check this I find other several implementations of it, for example matlab simple propagator of this guy: https://github.com/komrad36/SimpleOrbitalPropagator
I know that in matlab it’s dorman-prince method but difference in values only -+1-5, it’s okay for me.The problem is that with solo atmospheric drag satellite position change not quit correct, I guess. Following common sense satellite must go down sooner or later, but it fly away from Earth all the time position vector just increase. And I don’t find simple example to compare results with mine.
May be I am not right and propagate with only atmospheric drag is incorrect? Or it's error in code?
Raccoon is offline   Reply With Quote
Old 05-15-2018, 05:31 PM   #2
martins
Orbiter Founder
Default

What do you mean by incorrect? You are free to simulate whatever you like. And yes, adding drag should reduce the orbital energy and lead to a decaying orbit. If yours doesn't then there must be a bug.

Since you say your 2-body propagation works, the problem is likely in your drag term. Did you try to debug if it produces reasonable results? For example, is the vector returned by AccelDrag approx. opposite your y.vel value, as expected?

In the first instance, you can probably drop the second term in the calculation of your vel_rel value (i.e. non-rotating planet). It removes one potential source of error, and should still produce the expected result.

Are you sure that your definition of matrix T is as intended (i.e. transforms from global to local planet system and not the other way round)?
martins is offline   Reply With Quote
Old 05-15-2018, 10:48 PM   #3
cristiapi
Orbinaut
 
cristiapi's Avatar
Default

Quote:
Originally Posted by Raccoon View Post
 The problem is that with solo atmospheric drag satellite position change not quit correct, I guess. Following common sense satellite must go down sooner or later, but it fly away from Earth all the time position vector just increase.
It usually happens when the time step is too big.
I see:
for (int j = 0.0; j < 100; j++) {
...
dt = dt.addDays(j);
}
but I think that you should do:
for (int j = 0; j < 100; j++) {
...
dt = dt.addDays(t_step);
}
where t_step is a small constant time step.
cristiapi is offline   Reply With Quote
Reply

  Orbiter-Forum > Far Side of the Moon > Math & Physics

Tags
atmospheric drag, c++, propagation


Thread Tools

Posting Rules
BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts
Forum Jump


All times are GMT. The time now is 05:04 AM.

Quick Links Need Help?


About Us | Rules & Guidelines | TOS Policy | Privacy Policy

Orbiter-Forum is hosted at Orbithangar.com
Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2018, Jelsoft Enterprises Ltd.
Copyright ©2007 - 2017, Orbiter-Forum.com. All rights reserved.