#### Abstract

CubeSats, and small satellites in general, being small and relatively light, are sensitive to disturbance torques in the orbital environment. We developed a simulation tool that includes models of the major environmental torques and small satellite experiences in low Earth orbit, which allows users to study the attitude response for a given spacecraft and assist in the design of attitude control systems, such as selecting the magnet strength when using passive magnetic stabilization or designing the shape of the spacecraft when using aerodynamic attitude stabilization. The simulation tool named the Smart Nanosatellite Attitude Propagator (SNAP) has been public in precompiled form and widely used since 2010; this paper accompanies the release of SNAP’s source code with the inclusion of new models for aerodynamic torque and other new features. Details on internal models are described, including the models for orbit propagation, Earth’s magnetic field, gravity gradient torque, spacecraft shape modelling and aerodynamic torque, permanent magnetic dipole torque, and magnetic hysteresis. A discussion is presented on the significance of aerodynamic torque and magnetic hysteresis on a magnetically stabilized 3-unit CubeSat in the orbit of the International Space Station, from which many small satellites are deployed.

#### 1. Introduction

Small satellites have demonstrated significant utility in recent years at a significantly lower cost and risk tolerance than larger single spacecraft missions. Several factors challenge attitude control systems for a small spacecraft: (i) Being smaller and less massive make them more susceptible to environmental disturbance torques, which result in a large response to the satellite’s angular rate. (ii) Because of the larger number of small satellites deployed, space debris concerns favour deployment into low orbits, where orbit lifetime can be limited. At these altitudes, environmental disturbance torques such as gravity gradient and aerodynamic torques are more significant. Table 1 summarizes the general regions of influence for each of the environmental torques that can affect a small satellite. In orbits between 300 km and 650 km, in which many of CubeSat and small satellites are deployed, aerodynamic and gravity gradient torques can be expected to be significant and comparable in magnitude. (iii) With an increasing number of missions for small satellites venturing beyond low Earth orbit [1], solar pressure and gravitational torques can dominate attitude dynamics.

Passive satellite stabilization techniques are popular in the small-satellite community when only coarse pointing is required, because they require no power and no active control system while being fairly simple to implement in small satellites. Several methods of passive satellite stabilization are common [2]. Passive magnetic stabilization is achieved by using a magnetic dipole, typically achieved using permanent magnets, to track the Earth’s magnetic field in the same manner as a compass needle. Gravity gradient stabilization can be achieved by deploying a gravity gradient boom or designing the satellite with a highly asymmetric mass distribution such that the axis of minimum inertia aligns with the nadir vector. Finally, aerodynamic stabilization has been shown to be feasible for altitudes below 500 km, where velocity vector alignment can be obtained [3, 4].

The performance of a certain design is a function of its attitude dynamics under environmental torques, which depend on the expected orbit and altitude and the satellite’s geometry and mass properties. In order to design and quantify the performance of a certain satellite, a high-fidelity simulation is often implemented of the satellite parameters and the disturbance torques affecting it.

In the case of active control actuators such as reaction wheels, the design choice could be determined by the order of magnitude of the worst case expected torque on orbit, the minimum required slew times, and the desired pointing accuracy. A simulation to propagate the satellite in orbit may not be necessary, since the reaction wheels can be chosen to overcome the worst expected disturbance torques. Simpler simulations or calculations can be done on these special cases to quantify the drift and errors due to actuator resolution in order to quantify the pointing accuracy. In passive techniques, however, a simulator that models the orbit and environmental torques is very convenient to study passive stability system designs. Passively, stability is often achieved on only two of three rotation axes. Since rotation around the magnet axis in magnetic stabilization is uncontrolled, as well as roll in aerodynamic stability and rotations about the gravity gradient boom axis, it is difficult to predict the behaviour about these uncontrolled axes analytically. This motivates the development of a high-fidelity simulation to propagate the attitude in six degrees of freedom (DOF).

Notable related work from recent years includes the development of the NASA Jet Propulsion Laboratory (JPL) small-satellite dynamic testbed simulation [5] that includes numerous models for controllers, actuators, sensors, and dynamics. Simulations can address limitations in testbeds that typically model reduced degrees of freedom (e.g., using air bearings) and allow the modelling of entire systems including sensors, control algorithms, and actuation, in addition to attitude dynamics. The major torques affecting small satellites in LEO are gravity gradient, aerodynamic drag, and torques induced by the Earth’s magnetic field. To achieve passive stabilization, one of these environmental effects can be intentionally utilized for the satellite design to be greater than the other environmental torques. This concept is the essence of most passive attitude stabilization techniques, in which case the remaining torques act as disturbances and are the main cause of the pointing error. An attitude propagator would simulate the counteracting environmental torques, given a detailed satellite description, and can be used to tune the design of the stability system as well as to evaluate the amount of pointing accuracy that is achievable under the environmental disturbance torques.

Solar pressure is neglected here, for LEO, as it is typically at least one order of magnitude smaller than any of the other torques because the surface area of small satellites is typically small, as expected from Table 1. Beyond LEO, solar radiation pressure can become significant, causing forces that can be used for orbit manoeuvres [6] and torques that affect attitude [7].

#### 2. Smart Nanosatellite Attitude Propagator

This section describes the attitude propagator implemented in Simulink® that includes the orbital and attitude dynamic components. Simulink® is a graphical modelling environment that is a convenient design and simulation tool for time-varying dynamic systems. Several differential equation solvers are available within Simulink with varying parameters and capabilities, which makes it a convenient environment in which to spot and mitigate numerical instabilities or solver-specific issues and also to gain confidence in the simulation results by seeking consistent results across different differential equation solvers.

Figure 1 shows the high-level view of the Simulink implementation of our simulator, named the Smart Nanosatellite Attitude Propagator (SNAP). The satellite 6DOF body dynamics are implemented in the centre block, which has translational forces and rotational moments as inputs. The translational force is found using a two-body gravitational model to simulate orbital motion. Rotational moments are a sum of the environmental effects, namely, gravity gradient, aerodynamics, and magnetic effects. Feedback elements can be observed in the figure; the value of these forces and moments at each time step is a function of the satellite’s position in orbit and its attitude at the previous time step. Simulink’s solvers propagate the satellite’s state with time, given this description of the dynamics. Figure 2 is a screenshot of the graphical user interface, which outlines all simulation parameters and includes several import and export tools and convenient postsimulation plotting tools.

##### 2.1. Body Kinematics and Dynamics

Orbital dynamics are not the focus of this simulation environment, because attitude studies range over hours while orbital decay analysis spans years. Therefore, a simple two-body model is used to represent Earth’s gravity to keep the satellite in its orbit. The position in orbit is initialized using Keplerian orbital elements and thereon propagated using the two-body model. The satellite’s attitude is propagated using models for aerodynamic effects, gravity gradient, and magnetic effects.

The 3DOF rotational dynamics of the satellite body are modelled using a set of rotational dynamic equations in body-frame coordinates [7]: where is the satellite inertia matrix, is the body angular rate vector in body-frame, is the total external torque vector in body-frame, and is the quaternion attitude vector in body-frame. The quaternion attitude representation is used to describe the kinematics. Quaternions offer a compact representation within four parameters while avoiding singularities (division by zero for certain attitude angles), which, for example, the Euler angle representation suffers from.

The total external torque is found as the combination of the gravity gradient, aerodynamic, magnetic, and magnetic hysteresis models of the Orbital Environment Simulator:

These models calculate the torque components due to the respective environmental effects at a certain point in orbit as a function of the satellite mass and magnetic properties, the attitude at that point, the position in orbit, and the velocity in orbit. The individual torques are discussed in detail in the following sections.

##### 2.2. Gravity Gradient Torque

Gravity gradient torque is a significant source of angular moments for small satellites in low Earth orbit. The gravity gradient torque for an Earth-orbiting satellite is caused by differences in the distance to Earth across the satellite body; mass that is closer to Earth experiences higher gravitational attraction. For a given satellite geometry, the torque profile due to the gravity gradient is a function of attitude. An asymmetric body in a gravitational field will experience a torque tending to align the axis of the least inertia with the field direction [8]. The mass distribution of the satellite is adequately described in the inertia matrix. The torque to gravity gradient is modelled in the attitude propagator as [9, 10] where is the gravity gradient torque, is the unit vector toward nadir, is the distance from the centre of Earth to the satellite, is the inertia matrix, and is the geocentric gravitational constant. This equation is modelled in Simulink to calculate the gravity gradient torque at each time step, given the position in orbit, which defines the distance , and the current attitude, which is used to find the nadir vector expressed in body-frame coordinates.

##### 2.3. Aerodynamic Torque

Atmospheric drag for small satellites becomes a prominent source of disturbance and angular moments in the low part of LEO, at altitudes of 500 km and below. The atmospheric density decreases exponentially as a function of altitude, with atmospheric drag effects becoming minimal at higher altitudes. The aerodynamic torque for a certain area element can be calculated by [10] where is the aerodynamic torque, is the unit velocity vector, is the vector from the centre of pressure to the centre of mass, is the atmospheric density, is the satellite velocity, is the drag coefficient, and is the affected area.

The centre of pressure is not readily known. The aerodynamic torque for a certain attitude is a function of the area facing the velocity vector that is not shadowed by any other parts of the spacecraft body. Taking torque due to aerodynamics into account requires a method of representing the spacecraft geometry. Then an approach is needed to calculate the torque spacecraft experiences, given the geometric representation and the attitude of the satellite relative to the velocity vector.

We discretize the geometry of the satellite into volumetric elements, as shown for a hypothetical satellite in Figure 3 to illustrate the geometric representation. This particular design is a 3U CubeSat with side panels that deploy into a shuttlecock configuration, which can aerodynamically stably track the velocity vector, as shown in prior work [4].

A lookup table is generated that maps the attitude relative to the velocity vector to the amount of torque the satellite experiences at that orientation. This torque profile is generated before the simulation runtime to reduce the amount of computations and minimize the simulation duration. At runtime, the satellite’s angle to the velocity vector (the incoming wind) is computed from the current attitude, and the lookup table returns the torque associated with that deflection angle. The mapping table is generated across a full range of satellite rotations by considering elements directly facing the wind vector that are not shadowed by other satellite components. A form of numerical integration is performed by summing up the torque contributions of all the satellite elements to find the total torque that affects the satellite at a given attitude. Although shadowing is often ignored in the literature when the main body of the satellite is small relative to the dimensions of the fins, shadowing is an important factor to consider when finding a general solution and studying arbitrarily shaped satellites with unpredictable responses. This geometric representation becomes a very convenient tool for solving this type of problem.

The lookup table is used at runtime to obtain a torque factor given the satellite attitude at any given time step. That value is then scaled by the atmospheric density at that altitude and the satellite’s linear velocity to find the final torque affecting the satellite at that time step, according to the aerodynamic torque equation shown above.

##### 2.4. Magnetic Torque

A magnetic dipole in a magnetic field experiences an angular moment that aligns the dipole with the magnetic field lines, like a compass needle pointing north. Magnetic dipoles occur in spacecraft transiently due to the on-board electronics as well as in the structure of the spacecraft, which may contain a residual dipole; both of these effects can be a source of unwanted angular moment disturbance. Magnetic effects can also be used to control and stabilize the attitude of a satellite. Passive magnetic stabilization utilizes a set of permanent magnets to align the satellite with the Earth’s magnetic field and prevent random tumbling. The torque affecting a satellite due to a magnetic dipole interacting with the Earth’s magnetic field is modelled as [8]
where is the magnetic torque vector in body-frame, is the magnetic dipole moment vector in A·m^{2} in body-frame, and is the Earth magnetic flux density vector in body-frame.

The Earth’s magnetic field is modelled as a simple dipole (L-Shell Model) [9]. Given the position in orbit at a given simulation step, the local magnetic field from the Earth can be found using the dipole model and rotated to body-frame coordinates given the attitude of the satellite. At runtime, given the Earth’s magnetic field and the satellite’s orientation, the magnetic torque due to the magnetic dipole is found. Simulink’s Aerospace Blockset includes the more accurate World Magnetic Model [11], which can be used at the cost of reduced simulation speed.

##### 2.5. Magnetic Hysteresis Damping

Angular rate damping is a major problem in satellite attitude dynamics that is normally dealt with using active systems such as reaction wheels or magnetic torque rods. The nature of the space environment is such that there is no friction (or damping), i.e., torques proportional to and opposing the angular rate of the satellite are practically nonexistent. In a system sense, passive stabilization behaves as a second-order system with a very low damping factor, often resulting in unstable responses [9]. Magnetic hysteresis material is a completely passive solution for angular rate damping; it is, however, nontrivial to study and predict, motivating the implementation of a simulation environment [12].

Magnetic hysteresis material, which is magnetically “soft” material of low coercivity, can be magnetized by the Earth’s magnetic field and follows hysteresis patterns as it cycles in a magnetic field. This makes it suitable as a means for angular rate damping for small satellites in orbits with a significant magnetic field. Magnetic hysteresis is a physical property of ferromagnetic material. The material becomes magnetized when an external magnetic field is applied, forcing the magnetic domains on the atomic level to polarize. Depending on the magnetic remanence () of the material, it will retain a magnetic dipole of some strength when the external magnetic field is removed. The magnetic coercivity () of the material is the intensity of the external magnetic field required to diminish the magnetization from saturation to zero when applied against the polarity of the material. The lag (or “hysteresis”) in tracking the externally applied magnetic field caused by the coercivity and remanence of the ferromagnetic material results in energy lost as heat in the material. The phenomenon can be thought of as the magnetic dipoles having “friction” when their orientation changes.

Figure 4 shows a sample magnetization curve generated using the mathematical model of the hysteresis material for an oscillating satellite. The mathematical recipe was developed to simulate the NASA Passive Aerodynamically Stabilized Magnetically Damped Satellite (PAMS) satellite, launched in 1996, and is a set of first-order differential relationships that we implemented and introduced to SNAP [13]. The model is an improvement on previous parallelogram approximations used in previous versions of SNAP [2, 14], where the continuous nonswitching nature of the curve adds fidelity to the simulations.

When magnetic hysteresis material with low enough coercivity is selected, the Earth’s magnetic field is sufficient to magnetize and demagnetize it. This material can then provide an effective angular rate damping method for lightweight satellites. It is also a simple and completely passive solution; the only requirement is a calculation of the amount of hysteresis material required to achieve adequate damping. However, quantifying the amount of hysteresis material to include in a satellite design is challenging. The amount of damping caused by hysteresis material is not a fixed or calculated amount; it is a result of the behaviour of the hysteresis material interacting and cycling through the Earth’s magnetic field. Modelling and simulation is a convenient and effective way to study hysteresis material.

##### 2.6. Validation

SNAP has been validated in our prior work [14] by comparing simulated results to on-orbit data recovered. We validated using the design parameters and on-orbit observations by (i) NASA Passive Aerodynamically Stabilized Magnetically Damped Satellite (PAMS) from 1996 [15], (ii) Delfi-C3 from 2008 [16], and (iii) QuakeSat from 2003 [17]. Simulation results and on-orbit observations were consistent, which increased our confidence in the fidelity of SNAP’s simulations.

#### 3. Effect of Aerodynamic Torques on Passive Magnetic Stabilization

Many CubeSats are deployed at low altitudes (500 km and below) in order to limit orbit lifetime and reduce space debris. For example, several launches by the NASA Educational Launch of Nanosatellites (ELaNa) program were from ISS resupply missions or using the International Space Station’s nanosatellite deployer. At the ISS altitude (apogee 406 km and perigee 403 km, as of the time of this writing), aerodynamic torques become significant for asymmetric spacecraft whose geometric centre and centre of aerodynamic pressure are at an offset. Aerodynamic torque can be utilized to achieve alignment with the velocity vector, resembling a badminton shuttlecock, as presented in prior work based on SNAP [4].

At the low altitudes described above, aerodynamic torques act as a significant disturbance. Because of challenges in aerodynamic modelling, aerodynamic torques have often been neglected in prior work. In this section, we present a case study to illustrate the effect of aerodynamic torques on a passive magnetic attitude control system for a CubeSat-class satellite at that altitude.

##### 3.1. 3U CubeSat—Neglecting Aerodynamics

A 3U CubeSat deployed from the ISS is a compelling case. 3U CubeSats suffer the most significant gravity gradient and aerodynamic effects, being the least symmetric compared to smaller CubeSats. The ISS altitude is compelling, being relatively low and of stronger gravity and larger atmospheric density and also being a common orbit for CubeSat deployments.

Table 2 outlines the satellite description and simulation parameters used to illustrate the effect of neglecting aerodynamic effects, as is commonly done due to the complexity of modelling aerodynamics. The simulation parameters are based on the CubeSat Specification Revision 13 [18], the latest as of this writing. The satellite is designed to be magnetically stabilized using a permanent magnet with magnetic hysteresis material for damping. This form of stabilization is also very common for CubeSats with only coarse attitude control requirements, such as pointing a wide-beam antenna and preventing uncontrolled tumble.

Figure 5 shows the magnitude of the environmental torques the satellite experiences. In Figure 6, we observe that the magnetic field is successfully tracked by the magnet axes (the -axis) within approximately 30 degrees at a steady state, as the magnetic torque is two orders of magnitude larger than the strongest disturbance, i.e., the gravity gradient torque. Results are falsely optimistic when neglecting aerodynamics, as discussed next.

##### 3.2. 3U CubeSat—Including Aerodynamics

Next, we utilize the same simulation parameters as in the previous section and outlined in Table 2; however, we incorporate aerodynamics. As previously discussed, in aerodynamic modelling, the 3U CubeSat is modelled as a point cloud as shown in Figure 7. Rotating this point cloud through a range of attitude angles, a lookup table of aerodynamic torque values as a function of attitude is created, which is used at runtime by SNAP. The torque profile, shown in Figure 8, is normalized by velocity and atmospheric density, which are a function of the satellite’s position in orbit that SNAP propagates.

When aerodynamics is not neglected, we observe a significantly different attitude response. The magnetic stabilization system should be designed to overcome the strongest disturbance torque, which is the aerodynamic torque at this altitude, as shown in Figure 9. In this scenario at the steady state, the satellite tracks the magnetic field with oscillations within 45°, as shown in Figure 10. This is in contrast to 30° maximum steady-state error when aerodynamic torques are neglected, as presented in the previous section. Figure 11 shows the steady-state error plots for the two cases for clear side-to-side comparison.

##### 3.3. Hysteresis Material Underdamped Behaviour

Another case study we present is underdamped behaviour due to insufficient magnetic hysteresis damping. A large number of CubeSats use passive magnetic stabilization as the primary attitude control technique, which is based on including a permanent magnet in the spacecraft design to give it a significant magnetic dipole, which will align it with the Earth’s magnetic field similar to a compass needle. Damping is a critical component in order for the spacecraft to lose momentum and for oscillations to dampen before the spacecraft achieves stability in tracking the magnetic field. In passive magnetic attitude control, damping is typically achieved with magnetic hysteresis material, which is magnetically soft material that is magnetized by the Earth’s magnetic field and creates a dipole that lags the Earth’s field as the satellite rotates, creating a drag effect.

Selection of the correct amount of magnetic hysteresis material is critical to achieve a small steady-state error. From prior experience [2, 4], we note that using too little or too much hysteresis material has a significant effect on the settling time and steady-state error. To illustrate this, Figure 12 shows the attitude response for the same simulation described in Figure 10, with, however, half the amount of hysteresis material. A total of 0.05 cm^{3} of HyMu80 is used (0.025 cm3 along - and -axes).

#### 4. Conclusions

The modelling and implementation of an attitude propagator for Earth-orbiting small satellites are described. SNAP incorporates the major environmental torques from a small-satellite perspective (gravity gradient, magnetic, and aerodynamic), as well as magnetic hysteresis material, which is a passive solution to angular rate damping. The Simulink model propagates the satellite’s attitude given its mass properties and orbit parameters. At each time step, the various environmental torques are calculated given the magnetic field at that point, the velocity, the position in orbit, and the satellite orientation. The satellite position and orientation are modelled by a 6DOF body model. Simulink offers a variety of differential equation solvers to propagate the models and obtain attitude reports for analysis and animation.

A set of case studies of a 3U CubeSat deployed in the ISS orbit, for example, from the ISS nanosatellite deployer or from a resupply mission, are presented. For such a CubeSat, which is stabilized using permanent magnets in order to track the magnetic field, it is important to consider aerodynamic torques, especially at altitudes where they may be the most significant disturbance torque. Also, simulations show that the amount of magnetic hysteresis material used can significantly affect the attitude response, where too little results in underdamped behaviour and too much results in unresponsiveness and large steady-state error.

This paper accompanies the release of the source code of SNAP. In the current version, models for the major environmental torques and passive control methods, namely, magnetic and magnetic hysteresis, are included. Using magnetic hysteresis for damping, it is possible to investigate and design small satellites that are magnetically stable, as shown in this paper and aerodynamically stable [4]. Additionally, all models needed for gravity gradient stabilization are included as well. SNAP is designed to be expandable, using the MATLAB and Simulink programming environments, to include additional forces such as thrusters or solar pressure and additional torques such as those from active attitude control actuators, for example, reaction wheels and magnetic torque rods.

#### Data Availability

The release of the SNAP 3.0 models and source code and an initialization file with the simulation parameters outlined in Table 2 can be used to reproduce the results and figures presented in this paper. This manuscript accompanies the release of SNAP’s source code available at http://sar-lab.net/snap/.

#### Conflicts of Interest

The author declares that there is no conflict of interest regarding the publication of this paper.

#### Acknowledgments

SNAP extends from simulations developed between 2009 and 2012 as graduate research at the University of Kentucky advised by Dr. James E. Lumpp, Jr. Funding for this work and publication were covered by internal resources at the University of Michigan-Dearborn.