Simulink: spacecraft motion integration using spherical harmonic gravity model problem

Hi, I am trying to simulate a spacecraft motion under a gravity field comprising of spherical harmonics. I wrote a simple Simulink model that implements the equation of motion:

dr/dt = v; dv/dt = g;

where r,v, and g are the radius, velocity, and gravity acceleration vectors, respectively. the gravity acceleration vector is calculated twice, once using the Simulink block "Zonal Harmonic Gravity Model", and once using the Simulink block "Spherical Harmonic Gravity Model". The produced acceleration vectors from both blocks are compared and used separately as inputs for the equations of motion.

Here we see the utilized Simulink diagram:

The Simulink diagram includes 5 sections:

1. Gravity calculation in ECEF (Earth Centered Earth Fixed) frame using 21x21 spherical harmonics model
2. Gravity calculation in ECEF frame using J2 - J4 zonal harmonics
3. Setting the initial conditions
4. Integrating the equations of motion in ECI (Earth Centered Inertial) frame.
5. calculating the rotation matrix from ECEF to ECI and back. Notice that this transformation is an approximation because it does not include precession, nutation, and polar motion.

The code in the included embedded Matlab function "calc ECI2ECEF" is:

 function D_I2E = fcn(GMST)
cGST = cos(GMST);
sGST = sin(GMST);
D_I2E =  [cGST, sGST, 0;
         -sGST, cGST, 0;
           0  ,   0 , 1];

I am using a variable step solver with 1e-9 relative tolerance. All the other solver definitions are set as "auto".

Results

• When I feed the equations of motion with the"Zonal Harmonic Gravity Model" gravity acceleration, I get the followingresults (The "Spherical Harmonic Gravity Model" is calculatedfor reference, but not used):

The attached plots include:

1. Left-up plot the radius components in inertial frame.
2. Right-up plot the velocity components in inertial frame.
3. Left-down plot the "Spherical Harmonic Gravity Model"gravity acceleration components in inertial frame.
4. Right-down plot the "Zonal Harmonic Gravity Model"gravity acceleration components in inertial frame.

We see that when the "Zonal Harmonic Gravity Model" block is used, we get bounded radius and velocity. This is logical because the gravity is a conservative force. We see that the "Spherical Harmonic Gravity Model" calculated gravity acceleration is the same as the gravity from the block "Zonal Harmonic Gravity Model" (up to small variations which are not visible in the plotted scale).

• When I feed the equations of motion with the"Spherical Harmonic Gravity Model" gravity acceleration, I get the followingresults (The "Zonal Harmonic Gravity Model" is calculatedfor reference, but not used):

We see that both "Spherical Harmonic Gravity Model" and "Zonal Harmonic Gravity Model" produced similar results (up to small variations which are not visible in the plotted scale). Furthermore, we see that the radius diverges and consequently, the gravity magnitude decreases. The reason for this divergence is not clear to me. the radius and velocity should have been maintained in the same magnitude, as was previously demonstrated.

Does anyone have any idea why the radius diverged in the last scenario, as opposed to the previous scenario?

P.S.

I attached the Simulink file to this question.

Answer 1

[Sergio Tamayo]

It seems that the problem lies in the solver configuration. In essence, the step that is automatically selected by the variable step solver is too big, causing the divergence of the solution. A workaround for this is to set the maximum step size for the simulation to a fixed value like 1. A different approach, is to set the solver to a fixed step solver and set the step size to similar values. Once this is done, the position for the spherical harmonic does not diverge anymore.