# Comparing disturbance models

A customer recently asked us for help comparing the disturbance analysis available in the CubeSat Toolbox with the full model in SCT Professional, for a 3U CubeSat. That is, to compare CubeSatDisturbanceAnalysis to Disturbances. The CubeSat Toolbox uses a simplified model of the spacecraft geometry as a single set of fixed areas, nominally for a rectangular prism. The full model in SCT Pro allows for articulated and rotating bodies built of individual components. The CubeSat Toolbox has a subset of the environment and disturbance functions in Pro but includes drag and optical disturbances in Earth orbit. Given enough options for the two models it should be possible to get the exact same results. We will sketch out our process and discoveries in this post.

Creating a model of a 3U CubeSat in Pro was the easy part, as it is just a single box component, and verifying that the areas match those from CubeSatFaces is accomplished by a trivial command-line printout. Comparing the results of the optical and drag analysis is much more complex as there are so many variables:

• Atmosphere model: exponential, scale height, J70
• Solar flux model
• Earth albedo model
• Earth radiation model
• Satellite optical properties (specular, diffuse, absorptive)
• Attitude pointing (LVLH vs. ECI)

In order to compare the results, we had to call both disturbance models in a script and generate plots overlaying the resulting forces and torques. Here is the code defining the model for SCT Pro.

 m = CreateComponent( 'make', 'box', 'x', 0.1, 'y', 0.1, 'z', 0.3,... 'name', 'Core', 'body', 1, 'mass', mass, ... 'faceColor', 'gold foil', 'emissivity', thermal.emissivity,... 'absorptivity', thermal.absorptivity, 'sigmaT', optical.sigmaT,... 'sigmaA', optical.sigmaA, 'sigmaD', optical.sigmaD, 'sigmaS', optical.sigmaS,... 'inside', 0); BuildCADModel( 'add component', m ); 

You can see the thermal and optical properties that must be specified as well as the mass and dimensions. The spacecraft is inertially fixed and put into an equatorial orbit, so we would expect zero drag along the z axis and the x/y forces to oscillate at orbit rate. Then to call the disturbance model we generate a low Earth orbit, get the Earth environment and run the analysis, like so:

 [r,v] = RVOrbGen(el,t); e = EarthEnvironment( r, v, jD, d ); hD = Disturbances( 'init', g, e ); [fD,tD] = Disturbances( 'run', g, e, hD ); 

The EarthEnvironment function is where the guts of the space environment modeling occurs. This includes specifying albedo and radiation constants, calculating the atmospheric density over the orbit, computing the sun vector and solar flux magnitude, checking for eclipses, computing the Earth magnetic field, and correcting the inertial velocity for the rotation of the atmosphere for drag calculations. In the CubeSat toolbox, this data is computed inside the dynamical model, RHSCubeSat. The same steps of creating the model and calling the disturbance function are shown below.

 c = RHSCubeSat; c.mass = 3; c.inertia = InertiaCubeSat( '3U', c.mass ); [a,n,rho] = CubeSatFaces('3U',1); c.surfData.nFace = n; c.surfData.area = a; c.surfData.rFace = rho; for k = 1:6 % Radiation coefficients [absorbed; specular; diffuse] c.surfData.sigma(:,k) = [optical.sigmaA;optical.sigmaS;optical.sigmaD]; end c.atm = []; q = QZero*ones(1,size(r,2)); [tC, fC, h, hECI, fr, tq] = CubeSatDisturbanceAnalysis( c, q, r, v, jD ); 

Let’s look at drag first, as it proved to be the easiest to verify. The primary difference between the CubeSat model and full disturbance model for drag initially was the atmosphere model itself: CubeSat uses the Jacchia 1970 model by default, while EarthEnvironment specifies a scale height atmosphere. The Jacchia 1970 model accounts for changes in density with the time of day, resulting in an orbit rate periodicity; however it is computationally more intensive and not needed in preliminary analysis. The scale heights model depends only on altitude and is very quick. The CubeSat dynamic model already had an option to switch to the scale height atmosphere if desired, so we added that same option to the CubeSat disturbance analysis function. This promptly resulted in a close result between the models for the drag force.

Comparison of the drag force between the CubeSat Toolbox and SCT Pro disturbance models for a 3U CubeSat

A slight variation remains due to a difference in the transformation between the inertial frame and Earth-fixed frame between the two models. This transformation is used to account for the rotation of Earth’s atmosphere, as drag depends on the relative velocity between the surface and the air. CubeSat uses a fast almanac function, ECIToEF, to compute this matrix for a given Julian date. This model accounts for nutation but not as accurately as TruEarth does in Pro. The EarthEnvironment function in Pro, however, uses a simpler transformation using Earth’s rotational rate about the inertial z-axis, ignoring nutation. This accounts for the nonzero Z force in the CubeSat output, which can be seen to be four orders of magnitude less than the X/Y forces. Both approaches are valid for a preliminary analysis so we accept this small remaining difference.

Producing an equally close comparison for the optical forces unearthed a few bugs in the CubeSat version as well as differences in fidelity that are intentional. First, recall that there are three main contributions to optical forces in Earth orbit: solar flux; Earth albedo – that is, reflected flux; and Earth infrared radiation. The fluxes can be modeled simply as constants, or at higher fidelity by accounting for distance from the radiating body. The full disturbance model accounts for the change in solar flux over the year as the Earth moves in its orbit, which amounts to about 100 W/m2 or 7% of the average flux. The CubeSat environment model was not doing this, but since it was already calling the sun vector function which calculates the needed data, we decided to add it. The sun vector itself can be modeled a number of ways, with CubeSat providing a low fidelity almanac version and Pro a higher fidelity almanac option as well as an option for JPL ephemerides.

Making a temporary change in CubeSat to use the higher fidelity sun almanac provided closer results, but there were still differences in the optical forces. A check on the optical coefficients revealed that Disturbances assumed 100% diffuse reflection for planetary infrared radiation while the CubeSat version assumed 100% absorption. This was result of a misunderstanding of the model when the CubeSat version was created. The intention of the model is to assume 100% absorption, but the radiation has to be reemitted or the temperature of the body would increase to infinity. Hence the diffuse coefficient in the Pro model. So, the CubeSat version’s coefficients were corrected to match. This improved the agreement between the models but still not completely.

Further investigation uncovered a bug in the calculation of the planetary radiation flux and Earth albedo flux. We call this a copy/paste bug, as the distance scaling factor was in the code, but applied to the wrong variable when the CubeSat version was created from the Pro version – so the scaling was unitized out. Correctly applying the scale factor to the flux input (and not its unit vector) to the base solar force function, SolarF, resulted in exact agreement between the two models. The figure below shows the final differences with the CubeSat sun almanac function restored to its default value.

Comparison of the optical force between the simplified and full disturbance models for a 3U CubeSat

To summarize, we identified the following subtleties in making a direct comparison between these two disturbance models:

• Selecting the same atmosphere model, AtmDens2, for both cases
• Understanding the sun vector model in each case (SunV1 vs. SunV2)
• Adding the solar flux scaling with the Earth’s distance from the sun to the CubeSat model (1 line)
• Updating the optical coefficients for planetary radiation in the CubeSat model (to be diffuse)
• Correcting a bug in the scaling of the Earth albedo and radiation fluxes with altitude in CubeSat
• Accounting for Earth’s nutation in transforming velocity to the Earth-fixed frame

This provides a great illustration of how careful one needs to be in performing any disturbance analysis, as there are so many subleties in even simple models that must be recorded if results are to be replicated. Every constant, every scale factor, every coefficient, every source of ephemeris needs to be understood.

We were happy to provide the results of our analysis to our customer so that he could decide which models he wants to use in his analysis. This is an example of the technical support we can provide to all our SCT customers; if you have a question, just ask!

# Reaction Wheel Friction Models

Reaction wheels are used in many spacecraft for attitude control. A reaction wheel is a momentum exchange device because it controls the spacecraft by exchanging momentum with the rest of the spacecraft. Momentum is exchanged via a motor that is fixed to the spacecraft. As with all rotating parts it is subject to friction. Friction needs to be modeled as part of the design process.

The standard way to model friction is with three terms. One is damping which is proportional to wheel speed. The faster the wheel spins the more friction torque is produced. Ultimately, this limits the net control torque. At some speed the motor is just balancing the friction torque. The second component is Coulomb friction that is constant but flips signs when the wheel speed changes sign. the third is static friction. It is like Coulomb friction but only exists at zero speed.

An alternative friction model is known as the bristle friction model. This models friction as bristles that bend. It also has the same friction components described above but they are applied though the bristle state.

Both models can be made to produce similar results as shown in the following figure.

The static friction is clearly seen. The wheel speeds are nearly the same. The middle plot is of the bristle state. The problem with these models is when the torque is low and the wheel speed passes through zero. We then get limit cycling as shown below.

This is due to numerical error.

We can eliminate the limit cycling by using a very small integration time step with the bristle friction model. An alternative approach is to use the first model and multiply the sum of the static and Coulomb friction with a sigmoid, or s function.

The coefficient of the sigmoid function is beta. Very small betas remove the static friction, and all Coulomb friction, near zero speed. The large betas retain the form of the friction and eliminate the limit cycling!

These models can be found in the Spacecraft Control Toolbox 2015.1 . This particular script will be available in 2015.2 which will be released in July.

# Attitude Maneuvers with the CubeSat Control System

The CubeSat control system is designed to work with either thrusters or reaction wheels. It has a number of handy built in maneuver modes such as pointing at the sun, nadir pointing or pointing at a specific latitude and longitude on the ground. Here is the spacecraft shown in the VisualCommander interface.

The movie in the link below shows attitude maneuvers in the VisualCommander interface. The interface has pages for the various subsystems and attitude control system functions. We start by seeing the spacecraft in a polar orbit on the Summary page. The solar arrays are reorienting themselves so that their cell faces are pointed at the sun. We switch the 3D display to look along the boresight of the telescope. We then go to the ACS page and select a sun pointing maneuver. We go back to the Summary page and see that the sun appears in the display. We then return to the ACS page and command nadir pointing. The remainder of the movie shows the reorientation maneuver to nadir pointing.

For more information on our simulation frameworks including our real-time control system framework, ControlDeck, go to Simulation Framework page.

For more information on VisualCommander go to VisualCommander page.

You can also send us an email to find out more about our CubeSatControl System. All of these products are available now.

# Optical Navigation for Geosynchronous Transfer Orbits

Optical navigation, using the Earth’s chord width and angles between nadir and stars, is an alternative to GPS based navigation for autonomous spacecraft. The Optical Navigation System, developed under a NASA contract, is well suited for this application.

In the Spacecraft Control Toolbox, we provide an easy-to-use demo script in the 2014.1 release that shows you how to implement optical navigation. The system uses Unscented Kalman Filters (also known as sigma point filters) with non-linear dynamics and measurement models.

This demo uses our new UKF functions shown below:
 ukf.t = t; ukf = UKFPredict( ukf ); ukf = UKFUpdate( ukf ); 
ukf is a data structure that includes all filter information. Measurements are passed as data structures to UKFUpdate which have pointers to the measurement functions. In this way, any type of measurement can be used in the filter and introduced at any time.

The following plots show some results from the script. The first shows the orbit and the estimated orbit which are essentially the same.

The second shows the position errors. Of course, actual errors would depend on the accuracy of the sensors, particularly the Earth sensor. Great care needs to be taken when setting up the UKF parameters. As you can see, the largest errors are at perigee and if the UKF parameters are not set properly, the filter might think a hyperbolic orbit was a valid solution!

Check out our Spacecraft Control Toolbox page for more information on the 2014.1 release! More information about optical navigation can be found on our Deep Space Navigation page.

# Simulating Magnetic Hysteresis Damping

CubeSats have caused a renewed interest in magnetic control of satellites, and passive hysteresis damping in particular. Modeling actual hysteresis rods on a satellite is not trivial, and generally requires empirical data on the properties of the rods selected. Our newest CubeSat simulation demonstrates damping using rods in LEO. A permanent magnet is modeled using a constant dipole moment, and we expect the satellite to align with the magnetic field and damp. We evaluate the results by plotting the angle between the dipole and the Earth’s magnetic field and the body rates.

First, let’s verify the magnetic hysteresis model in the toolbox using the bulk material properties in orbit. We use a dipole model of the Earth’s magnetic field. The nice hysteresis curves below confirms that we are computing the derivatives of the magnetic field correctly in the body frame, which requires careful accounting of rotating coordinates. Also we stay within the saturation limits which means our magnetic flux derivatives are correct too.

Hysteresis curves from simulating magnetic hysteresis in orbit

We will assume the rods are 1 mm radius and 95 mm length, with rods placed perpendicular to each other and the permanent magnet. Three rods are used per axis. The apparent rod parameters are taken from the literature. The actual rods will not reach saturation while in orbit, so we will see a minor loop.

Minor loops from damping rods using apparent properties

The rods produce only a small amount of damping per orbit, so we have to run for many orbits or days to see significant damping – in some passive satellites, the total time allotted for stabilization is two months! In this case we test the rods’ ability to damp the torque induced by turning on a torque rod with a dipole of 1 AM2 and allowing the CubeSat to align itself with the magnetic field, starting from LVLH pointing.

Damping in LEO using hysteresis rods

Simulating the rods is time-intensive, with a timestep of about 4 seconds required – which makes a simulation of several days on orbit take several minutes of computation. Once performance of the rods has been verified, a simple damping factor can be substituted.

This new simulation along with the functions for hysteresis rod dynamics will be in the new version of our CubeSat Toolbox, due for release in June!

References:

1. F. Santoni and M. Zelli, “Passive magnetic attitude stabilization of the UNISAT-4 micro satellite”, Acta Astronautica,65 (2009) pp. 792-803
2. J. Tellinen, “A Simple Scalar Model for Magnetic Hysteresis”, IEEE Transactions on Magnetics, Vol. 34, No. 4, July 1998
3. T. Flatley and D. Henretty, “A Magnetic Hysteresis Model”, N95-27801 (NASA Technical Repoets Server), 1995

# Kepler Telescope Reorientation Maneuver

The Kepler telescope has suffered the loss of two reaction wheels. This means that it cannot use the wheels to control orientation about all three axes.

One option is to use thrusters and reaction wheels at the same time as actuators. Princeton Satellite Systems Core GN&C Bundle does just that.

We’ve simulated the system for the Kepler spacecraft

You can see a movie of a reorientation here:

# New Attitude Profiling Functions and Visualization in SCT v11

Creating attitude profiles just got easier! Satellites typically have multiple antennas and sensors that must be pointed in different directions at various times. We often want to point a sensor payload or a directional antenna at a certain location on Earth, while keeping the solar panels aligned with the sun and aiming the star camera away from bright areas of the sky. The group of new attitude profile functions in SCT v11 allow high-level directives to be defined, and facilitate the automatic computation of an attitude profile that meets the target alignment objectives while satisfying all pointing constraints. Detailed time-history plots and 3D visualization with playback enable you to explore and understand the attitude profile in depth.

The 2D plot below shows a time history of the rotation angle about the primary body axis. The primary body axis is aligned with the primary target. We can then rotate about this axis to align a secondary body axis as closely as possible with a secondary target. At the same, we have one or more pointing constraints which impose time-varying bounds on the rotation angle. The dark gray regions illustrate how these bounds change over time.

The 3D view below shows the orbital path (cyan) of the satellite about the Earth, with a CAD model at the current orbit location in the center of the figure. The sun vector is shown (yellow) and the Earth lighting is based on the sun location. The primary alignment vector (green) is directed towards a coordinate on the Earth, and the secondary alignment is pointed in the orbit-normal direction. Constraint directions are shown in red with angular sweeps to show their size. The sensor cone is a star camera that has to keep the sun, Earth and moon out of its field of view.