Design Module for a Spacecraft Attitude Control System

Design Module for a Spacecraft Attitude Control System Software Designed to Help Estimate Spacecraft Design Requirements 16.851 Satellite Engineering Motivation Spacecraft often have pointing requirements. Satellite antennas and optics will generally require that the spacecraft remain pointed at the desired target within a certain pointing tolerance. There are many environmental effects that disturb the spacecraft, and the attitude control systems (ACS) must be designed to account for this. The ACS system for this project will be designed for a three-axis stabilized spacecraft. This module will be useful for next semester when the x-ray telescope being designed will require three-axis stabilization. Problem Statement Design a tool that sizes attitude control system actuators for a three-axis stabilized spacecraft given disturbance torques created by various environmental effects. Environmental effects include gravity gradient, solar radiation, magnetic field, and aerodynamic forces. ACS actuators being investigated include momentum wheels and thrusters. Introduction MATLAB is used to evaluate the disturbance torques on the spacecraft and select and size appropriate ACS actuators to meet the pointing requirements of the spacecraft. For spacecraft without large slewing requirements, momentum wheels are primarily used to counteract disturbance torques. Momentum wheels and the propulsion system to dump the excess momentum will be the primary focus of the software module. The user will input information about the satellite including mass properties, physical dimensions, and orbit information (altitude, inclination, and eccentricity). The tool then examines the relevant environmental disturbance torques and finds the worstcase torque conditions for the specific spacecraft and orbit. Next, the tool will size the momentum wheels required to overcome worst-case disturbance torques. The masses of the momentum wheels are obtained from a database of available momentum wheel specifications. In addition, the module calculates the mass of the propulsion system required to dump momentum from the momentum wheels when they become saturated over the lifetime of the spacecraft. Finally, the tool is executed for a test case in order to check the validity of the module. In addition, other test cases are run in order to construct parametric design figures that show trends of the output involving disturbance torques and ACS mass as a function of spacecraft size and orbit parameters. Software Module ACS Environmental Disturbance Torque Tool Description of Code This MATLAB software tool sizes ACS angular momentum storage and dumping devices for a specified vehicle in any orbit. The tool only addresses disturbance torques from the main environmental sources: atmospheric drag, solar radiation, magnetic field, and gravity gradient. The ACS sizing is only meant to account for disturbance torques and does not address the design needs for vehicle slewing. The tool requires input data structures to describe the vehicle, orbit, and planet (Table 1). momentum. This ratio was chosen because cyclic momentum loads are known to drive the ACS design and this ratio matched well with guidance provided in SMAD. Table 1 Software tool inputs Parameter Description veh.dim [m], vehicle dimensions (x,y,z) veh.CG [m], vehicle CG offset from geometric center veh.mass [kg], vehicle mass veh.mat vehicle surface material code (see Appendix) veh.life [years], vehicle design life OE.a [m], orbit semi-major axis OE.e orbit eccentricity OE.i [deg], orbit inclination OE.Om [deg], argument of periapsis OE.om [deg], longitude of the ascending node planet.mu [m /s ], earth gravity constant planet.r_pol [m], earth polar radius planet.r_equ [m], earth equatorial radius 3 The cyclical angular momentum drives the momentum wheel design. The tool selects a momentum wheel system that can store the calculated cyclical momentum and then estimates the mass of the wheel. Next, the tool selects the thruster system to dump the accumulating secular torque for the duration of the vehicle life. This determines the thruster system mass and thrust capability. Finally, the tool outputs these calculated ACS specs for the worst case environmental disturbance torques. Table 2 shows the software tool output. 2 Table 2 Software tool output Parameter The vehicle is modeled as a rectangular prism and described by the three edge lengths, mass, center of gravity, and exterior surface material. Using the dimensions and mass, the vehicle’s moment of inertia can be determined using1: ⎡ m( y + z ) ⎢ 0 I=⎢ ⎢ 0 ⎣ 1 12 2 2 1 12 0 m( x 2 + z 2 ) 0 ⎤ ⎥ ⎥ 2 2 1 ⎥ 12 m ( x + y ) ⎦ 0 0 Units Cyclical angular momentum (per orbit) [Nms] Secular angular momentum (per orbit) [Nms] ACS Momentum wheel mass [kg] ACS Thruster system mass [kg] ACS Thruster thrust [N] (1) Program Execution The ACS ENV tool is executed using MATLAB. First, various constant parameters are specified by the user within the “ACS_env_ini.m” file (see Table 1). Next, the tool is executed by entering at the Matlab command prompt: Using the orbital elements, the orbital period and inertial position and velocity are calculated2. Throughout the orbit, worst case environmental disturbance torques are calculated for the four environmental sources. Details of the torque calculations will be addressed in the following sections. >> [h_cyc,h_sec,wheel_mass,thr_mass,thr_force] = ACS_env(OE,veh,planet); The five program outputs are specified in Table 2. Copies of the initialization file and all program modules can be found in the appendix. The total disturbance torque is then integrated for one complete orbit. This provides the worst case angular momentum imposed on the vehicle for one orbital period. Integrating the orbital torque instead of assuming an average from the maximum torque is advantageous because it incorporates the varying effects resulting from different orbit types. Two very different orbits might have the same maximum torque, but would in reality experience very different disturbance effects (i.e. circular orbit versus highly elliptic orbit). Aerodynamic Torque Module Requirements This MATLAB module, torque_aero.m, calculates the torque due to aerodynamic drag on the spacecraft. This drag is caused by the vehicle flying through the Earth’s atmosphere during close approaches to the planet. An atmospheric model is used to approximate the density of the atmosphere at each position of the spacecraft. Combined with geometric and aerodynamic information, the torque caused by drag is determined. For the angular momentum experienced by the vehicle, it is assumed that 80% is cyclic in nature and therefore only requires temporary storage during the orbit. The remaining 20% is secular and will continue to accumulate until the momentum dumping system applies an external torque and eliminates this 2 of this vehicle. Also, Chobotov recommends that a drag coefficient value of 2.2 be used as a conservative estimate.4 Description of Code The code uses an atmospheric model for the “upper atmosphere” from NASA. This model provides an approximation of the density of Earth’s atmosphere for use in drag calculations on the spacecraft.3 Inputs size (m): This input was a vector containing the three edge lengths of the rectangular parallelepiped-shaped spacecraft. The code also assumes the spacecraft is shaped as a rectangular parallelepiped of edge dimensions X, Y, and Z. The origin of the spacecraft is defined to be the geometric center of the 3-D shape. The coordinate system definition for the assumed spacecraft in the MATLAB module is shown below in Figure 1. V (m/s): This input is the instantaneous velocity of the spacecraft at any given time when the module is called by the main program. h (m): This input is the altitude of the spacecraft at any time in which the module is called by the main program. CG (m): This input is a vector of the center-of-gravity offset from the geometric center of the spacecraft. This is used to determine the worst-case moment arm on which the drag force is acting. Outputs Ta (N*m): This output is the disturbance torque on the spacecraft due to atmospheric drag in the spacecraft. This is sent back to the main program when this subroutine is called. Figure 1 S/C coordinate system definition Theory & Equations In order to estimate the atmospheric density at any given time during the orbit of the spacecraft, a density model provided by NASA was used. The model used was for the “upper atmosphere.” This meant that it was acceptable for use above 25km altitude. The equations below were used for this approximation.5 In addition, the center of pressure, the location at which the aerodynamic drag is acting, is assumed to be at the center of one face of the spacecraft. It is assumed that the worst case drag the spacecraft experiences is the case in which the vehicle has one face directly facing the direction of travel though the atmosphere. In addition, it is assumed that the face on which the drag is acting has the largest surface area on the spacecraft. This results in the worst-case drag possible on the vehicle. T = −131.21 + .00299h ⎡ T + 273.1⎤ p = 2.488⎢ ⎣ 216.6 ⎥⎦ p ρ= .2869[T + 273.1] The worst-case of the aerodynamic drag torque also results from the worst-case moment arm of the force acting on the vehicle. This worst-case moment arm is determined from the absolute value of the maximum CG offset from the geometric center. The geometric center is assumed to be the coordinate system origin for the spacecraft. Since the center of pressure acting on each face is at the center of the face, the largest moment arm acting around the center of gravity is the maximum distance in one axis of the center of gravity from the geometric center. (2a) −11.388 (2b) (2c) In the above equations, T is the atmospheric temperature (degrees Celsius), p is the atmospheric pressure (KPa), and ρ is the atmospheric density (kg/m3). The force acting on the spacecraft due to drag was determined by the following equation.6 F = 12 ρC D AV 2 Constants The only constant used in this module was the coefficient of drag, CD, which is assumed to be 2.2. This assumption was made due to a lack of experimental data about the drag and lift characteristics (3) The A in the above equation is the surface area of the face of the spacecraft normal to the airflow hitting the vehicle. 3 from the local vertical. To calculate the worst case scenario, θ is assumed to be 45o. The torque due to aerodynamic drag on the spacecraft is calculated using the following equation. Ta = F (max( Di )), where i = x, y, z Magnetic Field Disturbance Torques The magnetic field of the earth can cause torques on a satellite. The strength of the magnetic field torque on a satellite is a function of the satellite’s position. The magnetic field intensity B is found by the equation: (4) The torque, Ta, is calculated by multiplying the drag force on the spacecraft, F, by the maximum distance of the center of gravity from the geometric center out of the x, y, or z directions, Di. This would allow for the largest moment arm and therefore would result in the highest disturbance torque due to aerodynamic drag. B= Fs As (1 + q ) cos i(cps − cg ) c (7) where B0 is the magnetic field at the equator at the Earth’s surface, 0.3 gauss, R is the radial distance from the Earth, and λ is the magnetic latitude. The worse case disturbance torque from the magnetic field on the satellite can then be found by multiplying B by the residual dipole of the satellite, D: Torque from Solar Radiation Solar radiation pressure if found from the equation: Tsp = Bo 1+ sin 2 λ 3 R (5) Tm = DB where Fs is the solar constant, 1367 W/m2, c is the speed of light, 3x108 m/s, As is the surface area of the satellite facing the sun, q is the reflectance factor, i is the angle of incidence to the sun, cps is the location of the center of solar pressure, and cg is the center of gravity of the satellite. The units of solar radiation pressure are Newton-meters. (8) The vehicle residual magnetic dipole (D) is assumed to be 1 Am2 for this problem. Mass in Reaction Wheels Using the cyclic angular momentum calculated by the tool, the reaction wheel mass can be estimated. To determine mass of a reaction wheel, data was collected on many existing, commercial reaction wheels. A fourth-order polynomial curve was fit to this data using MATLAB, to compare momentum storage in reaction wheels and their weight in kilograms. From the curve, the approximate mass of a momentum wheel given its momentum storage can be determined. The figure below illustrates the data collected (red points) as well as the polynomial curve fit used (shown as the blue line). The worst case is for an incidence angle of zero degrees, so i=0 was always used in calculations. Also in a worst case situation, As would be the surface area of the largest face of the satellite. The value for q is based on the material of the satellite, which is specified on input. Gravity-gradient Disturbance Torques The gravity field of the Earth can place unwanted torques on a satellite. When a satellite is affected by the Earth’s gravity field, the longitudinal axis becomes aligned toward the center of the Earth. The strength of the torque is a function of the distance of the satellite from the Earth. Also, the torque is symmetric about the nadir vector of the satellite. The equation for finding the worst case gravity-gradient torque is: Tg = 3µ I z − I y cos( 2θ ) 2R3 (6) where µ is the gravitational constant of the Earth, 3.986*1014 m3/s2, R is the orbit radius, Iz and Iy are the moments of inertial about the z and y axes of the satellite, and θ is the maximum deviation of the z-axis Figure 2 Mass vs. angular momentum capability for momentum wheels 4 the mass of the propulsion system is 85% propellant.19 This, by sizing the propellant, the total system mass can be determined by dividing by 0.85. Sources used for collecting data on reaction wheels are included in the references section of this paper.7,8,9,10,11,12,13,14,15,16,17 It should also be noted that it is assumed in this module that the time for each thruster firing is one second. ACS Propulsion System Design Module Constants The only constant used in this module is the acceleration due to gravity at the Earth’s surface, g. This is assumed to be 9.8 m/s2. Requirements This MATLAB module, called prop_system.m, calculated the thrust required to dump momentum at each momentum dump thruster firing as well as the total ACS propulsion system mass required for momentum dumping. Inputs size (m): This input was a vector containing the three edge lengths of the rectangular parallelepiped-shaped spacecraft. Description of Code This code calculates the worst-case moment arm of an ACS thruster about a spin axis of the spacecraft and combines that moment arm with the required thrust for dumping momentum to size the propulsion subsystem. CG (m): This input is a vector of the center-of-gravity offset from the geometric center of the spacecraft. This is used to determine the worst-case moment arm on which the drag force is acting. The code assumes the spacecraft is shaped as a rectangular parallelepiped of edge dimensions X, Y, and Z, with thrusters arranged at locations shown in the figure below. The thrusters are the red triangles in the figure below. lifetime (years): This input is the lifetime of the satellite. It is used to determine how many thruster firings to dump momentum will be needed throughout the life of the spacecraft. H (N*m*s): This input is the stored maximum momentum in any one momentum wheel. It is the momentum which will need to be dumped by using the thruster firing of the ACS system. sat_rate (days/saturation): This input is the time it takes for the wheels of the spacecraft to become saturated with angular momentum. It is at this point that the momentum wheels have no more capacity to control the attitude of the spacecraft by spinning up any faster. Outputs F (N): This output is the force required for the thruster to impart on the spacecraft in order to dump the required momentum. p_mass (kg): This output is the total mass of the propulsion system used to dump momentum when the momentum wheels become saturated throughout the lifetime of the spacecraft. Figure 3 Assumed thruster locations on spacecraft A Hydrazine monopropellant propulsion system is assumed to be used for momentum dumping for this spacecraft. Although the specific impulse, Isp, for Hydrazine can be as high as 245 seconds,18 a conservative value for the Isp for Hydrazine is assumed in this module. This conservative Isp value is 200 seconds. Theory & Equations The force required to dump momentum is obtained by using the following equation. F= Although the thrusters and propulsion system plumbing are not directly sized, a conservative value for propulsion system mass is estimated. It is assumed that 5 H Lt (9) In the above equations, L is the moment arm of the thruster to the required spin axis (meters) and t is the thruster firing time. Table 3 FireSat-like test case parameters The propellant mass, mP, and total propulsion system mass, p_mass, are then calculated using the following equations. 3 * lifetime * 365.25 (10) total _ pulses = sat _ rate F (total _ pulses )t (11) mP = I sp g m p _ mass = P 0.85 (12) Total_pulses is the total number of thruster firings required throughout the lifetime of the spacecraft to counteract the environmental disturbance torques. The numerator of the total_pulses equation has a factor of 3 in it because it is assumed that all three wheels will need to be desaturated at each time momentum is required to be dumped. The Isp in the equation above is for Hydrazine. The 0.85 factor is explained in the “Description of Code” section. Parameter Value Description veh.dim [ 1.7 1 1.7] [m], vehicle dimensions (x,y,z) veh.CG [0.2 0 0] [m], vehicle CG offset from geometric center veh.mass 200 [kg], vehicle mass veh.mat 9 vehicle surface material code (corresponds to 0.63 reflectance) veh.life 4 [years], vehicle design life OE.a 7,078,000 [m], orbit semi-major axis OE.e 0.0 orbit eccentricity OE.i 45 [deg], orbit inclination OE.Om 0 [deg], argument of periapsis OE.om 0 [deg], longitude of the ascending node planet.mu 3.986e14 [m3/s2], earth gravity constant planet.r_pol 6,357,000 [m], earth polar radius planet.r_equ 6,378,000 [m], earth equatorial radius Based on the input parameters in the above table, the MATLAB module calculated the cyclical and secular angular momentum required to counteract disturbance torques on the spacecraft for one orbit, the required momentum wheel mass for the ACS system, the thrust required for each instance angular momentum of the spacecraft needs to be dumped, and the ACS propulsion system mass required for dumping angular momentum for the life of the spacecraft. Results Module test case In order to verify that the MATLAB code is working properly, the module was used to determine the disturbance torques and calculate the required ACS mass for a test case. The test case used was the main spacecraft example in Space Mission Analysis and Design, by Wertz and Larson. This example is the FireSat satellite. The main parameters used to simulate the FireSat example to test this module with are shown below in Table 3. A plot of the orbit of the FireSat spacecraft example is shown below in Figure 4. Figure 4 FireSat orbit 6 are not included in SMAD and therefore can not be compared to the module solution. The disturbance torques imposed on the spacecraft are shown in the polar plot in following figure. Zero degrees on the plot corresponds to the orbital ascending node. Module output The following section will display the capabilities of the MATLAB module by presenting data collected by running the module for various orbits and spacecraft sizes. First, required angular momentum to counteract cyclic disturbance torques on the spacecraft is calculated for various orbit eccentricities and altitudes at perigee. The results are shown below in Figure 6. Figure 5 Disturbance Torques for FireSat Test Case The curves in the plot in Figure 5 show the disturbance torque due to gravity gradient is the largest contributor to the overall disturbance torque on the spacecraft. The torque due to magnetic fields is the second-largest contributor to the total disturbance torque. Torque due to solar radiation and aerodynamic drag are the two minor contributors to overall disturbance torque. Solar pressure on the spacecraft is a weak force, which results in the torque from solar radiation being small. The torque due to drag is low because the spacecraft is in an 700km orbit. The Earth’s atmosphere at an altitude of 700km is almost nonexistent. Figure 6 Cyclic angular momentum vs. eccentricity for various radii of perigee The curves in the figure show several interesting trends which correspond with reality. First, the highlyeccentric orbits generally have lower required angular momentum storage capability. This is due to the fact that these orbits are far away from the Earth for most of the orbit. This nearly eliminates disturbance torques from drag, magnetic fields, and gravity gradient. The remaining torque, caused from solar radiation, is a small contributor to disturbance torque. A comparison of the results from SMAD and results from the MATLAB code is shown below in Table 4. Table 4 FireSat-like test case solution vs. Module Value SMAD Solution Cyclical angular momentum (per orbit) 0.3845 [Nms] 0.4 [Nms] Secular angular momentum (per orbit) 0.0961 [Nms] N/A ACS Momentum wheel mass 1.41 [kg] N/A ACS Thruster system mass 2.46 [kg] 2.43 [kg] Parameter The reason the angular momentum increases slightly near the high end of eccentricity is most likely due to the fact that the satellite in those orbits will be passing by Earth at extremely high speeds compared to the smaller-eccentricity orbits. This high speed near Earth may contribute to much larger torque due to drag for that portion of the orbit. Another interesting observation from the figure above is that the disturbance torque for satellites around 400 or 500km altitude is much greater than that of satellites in slightly higher low-Earth orbits of 600km. This shows a potential significant cost savings for an ACS system if a LEO satellite were placed in a 600km orbit as opposed to a 400km orbit, for example. It can be seen in the above table that the values for cyclical angular momentum and ACS thruster system mass are nearly identical. The values for secular angular momentum and ACS momentum wheel mass 7 Figure 7 below graphs ACS mass versus orbit eccentricity for various radii of perigee. Figure 8 Angular momentum vs. altitude for various inclinations of circular orbits Figure 8 shows a clear trend that for all circular orbit inclinations, there is an altitude at which the required angular momentum storage of the ACS system has reduced significantly and levels-off. This altitude is around 525km. This means there could be a significant benefit to putting LEO spacecraft in circular orbits above 525km in altitude in order to minimize the angular momentum capacity required for the ACS system. Figure 7 ACS mass vs. eccentricity for various radii of perigee Figure 7 is interesting because it shows that the ACS mass for highly-eccentric orbits does not follow the exact same trend as that shown in the previous figure. The reason the ACS mass is reduced for the highlyeccentric orbits and does not continue with the same trend as the angular momentum is due to the fact that the highly-eccentric orbits have large orbit periods. Since the lifetime of the spacecraft is kept constant for these trend studies, a longer orbit period would result in fewer ACS angular momentum dumping situations during the lifetime of the satellite. It would take more time for a longer orbit for the attitude control system to become saturated with angular momentum and then require thrust to dump the momentum. This reduces the mass of the ACS propulsion system and produces the trends seen in the above figure. In addition, there is a difference in angular momentum between various inclinations. It is especially noticeable at lower altitudes. For example, at 400km altitude, the angular momentum required for the 0 degree inclination (equatorial) orbit is approximately 20 N*m*s greater than that of the 90 degree inclination (polar) orbit. This trend also reverses itself after the 525km altitude mark and the polar orbit becomes the maximum angular momentum case. The figure below shows the ACS mass versus orbit altitude for circular orbits of various inclinations. The figure below plots angular momentum versus altitude for various inclinations for circular orbits. 8 and may experience less torque due to the Earth’s magnetic fields. Figure 10 also clearly shows that as the spacecraft grows in size, the angular momentum requirement increases as well. This is due to increased aerodynamic drag. Figure 11 below shows the ACS mass versus vehicle volume for the same orbits as in the previous figure. Figure 9 ACS mass vs. altitude for various inclinations of circular orbits Figure 9 clearly exhibits the same trend of that shown in Figure 8. This shows that the ACS mass is directly related to the required angular momentum storage of the spacecraft for circular orbits. The figure below shows the angular momentum storage required versus vehicle size for three different types of orbits. The vehicle density was held constant while volume and the corresponding mass were varied. The three investigated orbits are the Molniya orbit, a circular LEO orbit of 400km, and a GEO orbit. Figure 11 ACS mass vs. vehicle volume for three orbit types Figure 11 above shows the trend that the ACS mass increases roughly at the same rate as the angular momentum requirement for orbits of the same type as vehicle volume is increased. Conclusion A MATLAB module was created which took inputs of orbital parameters, spacecraft dimensions, and the spacecraft environment and output cyclical angular momentum, secular angular momentum, ACS momentum wheel mass, and ACS thruster system mass. The module was checked by running the FireSat example. In addition, other cases were checked with the module and realistic data was output from the module. Figure 10 Angular momentum vs. vehicle volume for three orbit types The values of angular momentum and ACS mass are of use for preliminary design of an attitude control system. The designer can input preliminary information about a spacecraft design and get rough numbers for preliminary sizing of the ACS system and its impact on the spacecraft mass budget. Although the vertical axis in the above figure is in a logarithmic scale, it can be seen that the angular momentum required for a 400km LEO orbit is much greater than that for a Molniya or a GEO orbit. The main reason for this is the fact that satellites in Molniya or GEO orbits spend all or most of an orbit period far from Earth. This means the spacecraft in those orbits will not experience much aerodynamic drag, will experience a reduced torque due to gravity gradient, Future Work The module can be expanded to include reaction wheels, control-moment gyros, and magnetic torquers. 9 In addition, the database of momentum wheels can be expanded to include more wheels and include other ACS actuator specifications as well. Another area for future work is to validate the module by running a real life example of a satellite and comparing the results to experimental data. 10 Appendix A FireSat-like Test Case Input %############################### %## ACS Env %## %## Initialization file %## %# Vehicle Properties #% veh.dim = [1.7 1 1.7]; %[m], Length of edges on rectangular-prism shaped vehicle (length, width, depth) veh.CG = [0.2 0 0]; %[m], center of gravity offset from geometric center veh.mass = 200; %[kg] (minus ACS system) veh.mat = 9; %Surface material code veh.life = 4; %[yrs], Vehicle lifespan %# Orbital Elements #% OE.a = 7078000; %[m], semi-major axis OE.e = 0.0; % eccentricity OE.i = 45 *pi/180; %[rad], inclination OE.Om = 0 *pi/180; %[rad], argument of periapsis (angle from ascending node to periapsis) OE.om = 0; % [rad],longitude of the ascending node (angle between x and asc. node) %# Planet Properties #% planet.mu = 3.986e14; %[m^3/s^2], Earth gravity constant planet.r_pol = 6357000; %[m], Polar radius planet.r_equ = 6378000; %[m], Equitorial radius 11 ACS Env Main Program % ## ACS sizing for Environmental Disturbance torques ## % % Input: % OE - Orbital elements for vehicle orbit % veh - vehicle parameters % planet - planet parameters % % Output: % ang_mom_cyc - cyclical ang momentum storage requirement [Nms] % ang_mom_sec - secular ang momentum received per orbit [Nms] % w_mass - spec for momentum wheel mass [kg] % t_mass - spec for thruster mass [kg] % thrust - spec for thruster thrust [N] % function [ang_mom_cyc, ang_mom_sec, w_mass, t_mass, thrust] = ACS_env(OE,veh,planet); %## Calculate vehicle moment of inertia ##% Ixx Iyy Izz I = = veh.mass*(veh.dim(1)^2 + veh.dim(2)^2)/12; %[kg*m^2] = veh.mass*(veh.dim(1)^2 + veh.dim(3)^2)/12; %[kg*m^2] = veh.mass*(veh.dim(2)^2 + veh.dim(3)^2)/12; %[kg*m^2] diag([Ixx Iyy Izz]) %# Calculate time for one orbit t = 2*pi*sqrt(OE.a^3/planet.mu); %[sec] %# Calculate Solar radiation torques Ts = torque_solar(veh.dim, veh.CG, 0, veh.mat); ang_step = pi/50; %[rad] ang_range = [0:ang_step:2*pi]; ii = 1; %## Calculate max disturbance torque around one complete orbit ##% for ang = ang_range %# Calculate orbital postion and velocity [r,v] = oe2rv([OE.a OE.e OE.i OE.Om OE.om ang], planet.mu); v_mag = norm(v); %[m/s], Calculate velocity magnitude R(ii,:) = r; %[m] V(ii,:) = v; %[m/s] V_MAG(ii,:) = v_mag; E = acos((OE.e + cos(ang))/(1+OE.e*cos(ang))); %Eccentric anomoly time(ii,:) = sqrt(OE.a^3/planet.mu)*(E-OE.e*sin(E)); %time to 'ang' %# Post-process orbit elevation (latitude) aa = sqrt(r(1)^2 + r(2)^2); lat = atan2(r(3), aa); %[rad] LAT(ii,:) = lat; %# Post-process altitude r_planet = (planet.r_pol*planet.r_equ)/... sqrt(planet.r_pol^2*cos(lat)^2 + planet.r_equ^2*sin(lat)^2); %[m], Calculate planet radius assuming oblate sphereoid alt = norm(r)-r_planet; %[m], Subtract planet radius from vehicle position vector magnitude ALT(ii,:) = alt; %# Calculate Aerodynamic torques 12 Ta = torque_aero(veh.dim, v_mag, alt, veh.CG); TA(ii,:) = Ta; %# Calculate Magnetic torques Tm = torque_magnetic(lat, norm(r), r_planet); TM(ii,:) = Tm; %# Calculate Gravity torques Tg = torque_gravity(norm(r), planet.mu, I); TG(ii,:) = Tg; TS(ii,:) = Ts; %# Sum all disturbance torques T(ii,:) = Ts + Ta + Tm + Tg; %[Nm] ii=ii+1; %Increment counter end %## post-process time values ##% max_t = ceil(length(time)/2); for jj = [max_t+1: length(time)]; time(jj) = 2*time(max_t)-time(jj); end %## Integrate max torques around orbit to find total ang mom ##% ang_mom = trapz(time,T); %[Nms], total angular momentum around one complete orbit ang_mom_cyc = 0.8 * ang_mom; %[Nms], cyclical angular momentum per orbit ang_mom_sec = ang_mom - ang_mom_cyc; %[Nms], Secular angular momentum per orbit %## Size ACS actuators for cyclical momentum storage##% wheel_data = get_wheel_data; %Loads Reaction wheel data (mass vs Nms) w_mass = polyval(wheel_data, ang_mom_cyc); %## Size ACS thrusters for secular momentum dumping ##% orb_sat = 1; %[orbits/saturation] day_sat = orb_sat*t/86400; %[day/saturation] [thrust, t_mass] = prop_system(veh.dim, veh.CG, veh.life, ang_mom_sec, day_sat); %## Plot results ##% plot_planet_3D(R,planet); figure(23); polar(ang_range', T, 'k') hold on polar(ang_range', TS, 'b'); polar(ang_range', TA, 'm'); polar(ang_range', TM, 'r'); polar(ang_range', TG, 'g'); hold off legend('total', 'solar', 'aero', 'magnetic', 'gravity'); xlabel('[deg]') 13 torque_aero.m % Bill Nadir % 16.851 Satellite Engineering % 10/11/2003 % Module for calculating external spacecraft torque caused by Aerodynamic forces function Ta = torque_aero(size,V,h,CG) % Here the force on the S/C, F, is calculated % % % % % % INPUTS size = V = h = CG = edge lengths of the S/C: vector (x,y,z) side lengths (meters) S/C velocity (m/s) S/C altitude (m) location of the center of gravity (x,y,z) for the S/C (assumed offset from the geometric center of (0,0,0)) (m) % OUTPUT % Ta = Torque on S/C due to aerodynamic drag (Nm) % % % % rho is the atmospheric density at the location of the S/C An atmospheric model for the upper atmosphere (h>25000m) is used to approximate the density of the upper atmosphere T is the atmospheric temperature, p is atmospheric pressure T = -131.21 + .00299*h; % in deg C p = 2.488*(((T + 273.1)/216.6)^-11.388); % pressure in KPa rho = p / (.2869*(T + 273.1)); % in kg/m^3 % C_D is the drag coefficient of the cube-shaped S/C (assumed = 2.2) C_D = 2.2; % Cpa % % = location of the center of aerodynamic pressure (x,y,z) (assumed at the center of the face of one side of the cube which is facing directly into the atmosphere = max drag) % Here the surface areas of the sides of the S/C are determined % This is used to assume the worst-case drag on the vehicle % [x*z y*z x*y] => find max area = [size(1)*size(3) size(2)*size(3) size(1)*size(2)]; max_area = max(area); F = 0.5*rho*C_D*(max_area^2)*(V^2); % here the external aerodynamic torque on the S/C is calculated Ta = F*max(abs(CG)); 14 torque_gravity.m % % % % % % % % % % Disturbance torque from gravity gradient Input: r = vehicle radius [m] mu = planet gravity constant [m^3/s^2] I = vehicle moment of inertia [kg*m^2] Output: T_grav = gravity gradient torque [Nm] function T_grav = torque_gravity(r, mu, I) %# Max moment Imax = max(diag(I)); %[kg*m^2] %# Min moment Imin = min(diag(I)); %[kg*m^2] %# Angle deviation from vertical theta = 45*pi/180; %[rad], worst case angle chosen %# Calc gravity gradient torque T_grav = 3*mu*sin(2*theta)*(Imax - Imin)/(2*r^3); %[N*m] torque_magnetic.m % % % % % % % % % % Disturbance torque from magnetic field Input: lat = vehicle latitude [rad] r = vehicle position vector magnitude [m] re = earth radius [m] Output: T_mag = magnetic field torque [Nm] function T_mag = torque_magnetic(lat,r,re) %# Earth magnetic field (approx as dipole) B = (1 + sin(lat)^2)^(0.5) * 0.3/((r/re)^3); %[gauss] B_t = B*1e-4; %[tesla], [N/(A*m)] %# Vehicle residual dipole D = 1; %[A*m^2] %# Mag torque T_mag = B_t*D; %[Nm] 15 torque_solar.m function T_solar = torque_solar(A, CG, i, mat) % function to computer solar radiation pressure % % % % % INPUTS A: vector describing size of object CG: distance from center of solar pressure to center of mass (m) i: angle of incidence of the Sun (radians) mat: ID of material on outside of craft % OUTPUT % T_solar: solar radiation pressure, in N*m % some constants: % speed of light, m/s c = 3*10^8; % solar constant, W/m^2 F_s = 1367; % get reflectance, q, from file based on material used tmp = xlsread('material_prop.xls','abs'); q = tmp(mat,3); % find surface area of largest face of orbit A_s = A(1)*A(2); if(A(1)*A(3) > A_s) A_s = A(1)*A(2); end if(A(2)*A(3) > A_s) A_s = A(2)*A(3); end F = (F_s/c)*A_s*(1 + q)*cos(i); T_solar = F*(max(abs(CG))); 16 get_wheel_data.m function p = get_wheel_data wheel(1) = struct('name', 'Teldix RSI 01-5/15', 'ang_moment', 0.04, 'mass', 0.6); wheel(2) = struct('name', 'Teldix RSI 01-5/28', 'ang_moment', 0.12, 'mass', 0.7); wheel(3) = struct('name', 'LeoStar', 'ang_moment', 4.7, 'mass', 3.628); wheel(4) = struct('name', 'Dyncon MicroWheel 200', 'ang_moment', 0.18, 'mass', 0.93); wheel(5) = struct('name', 'Honeywell HR12', 'ang_moment', 50, 'mass', 9.5); wheel(6) = struct('name', 'Honeywell HR14', 'ang_moment', 75, 'mass', 10.6); wheel(7) = struct('name', 'Honeywell HR16', 'ang_moment', 100, 'mass', 12); wheel(8) = struct('name', 'Honeywell Miniature Reaction Wheel', 'ang_moment', 1.0, 'mass', 1.3); wheel(9) = struct('name', 'Honeywell HR0610', 'ang_moment', 12, 'mass', 5.0); wheel(10) = struct('name', 'Teldix DR23-0', 'ang_moment', 23, 'mass', 6.9); wheel(11) = struct('name', 'Teldix RDR68-6', 'ang_moment', 68, 'mass', 9.1); wheel(12) = struct('name', 'Teldix RSI 25-75/60', 'ang_moment', 25, 'mass', 6.3); wheel(13) = struct('name', 'Teldix RSI 68-75/60x', 'ang_moment', 68, 'mass', 8.5); wheel(14) = struct('name', 'Teldix RSI 4-75/60', 'ang_moment', 4, 'mass', 3.7); wheel(15) = struct('name', 'Teldix RSI 12-75/60x', 'ang_moment', 12, 'mass', 4.85); wheel(16) = struct('name', 'Teldix RSI 18-220/45', 'ang_moment', 18, 'mass', 6.3); wheel(17) = struct('name', 'Teldix RSI 30-280/30', 'ang_moment', 30, 'mass', 8.5); wheel(18) = struct('name', 'Teldix RSI 68-170/60', 'ang_moment', 68, 'mass', 8.9); wheel(19) = struct('name', 'Teldix RSI 02-25/30', 'ang_moment', 0.2, 'mass', 1.7); wheel(20) = struct('name', 'Teldix RSI 04-25/60', 'ang_moment', 0.4, 'mass', 1.7); wheel(21) = struct('name', 'Teldix RSI 1.6-25/60', 'ang_moment', 1.6, 'mass', 2.4); for(i=1:length(wheel)) %plot(wheel(i).ang_moment, wheel(i).mass, 'r*'); %hold on; ang(i) = wheel(i).ang_moment; mass(i) = wheel(i).mass; end [p,s] = polyfit(ang, mass, 4); %f = polyval(p, ang); %plot(ang, f, 'g*'); 17 oe2rv.m % % % % % % % % % % % % % % % CREDIT: Christopher D. Hall http://www.aoe.vt.edu/~cdhall/ oe2rv.m Orbital Elements to r,v [r,v] = oe2rv(oe,mu) oe = [a e i Om om nu] r,v expressed in IJK frame a = semi-major axis e = eccentricity i = inclination Om = argument of periapsis om = right ascension of the ascending node (longitude of ascending node) nu = true anomaly (at epoch). ***(location on orbit)*** function [ri,vi] = oe2rv(oe,mu) a=oe(1); e=oe(2); i=oe(3); Om=oe(4); om=oe(5); nu=oe(6); p = a*(1-e*e); r = p/(1+e*cos(nu)); rv = [r*cos(nu); r*sin(nu); 0]; % in PQW frame vv = sqrt(mu/p)*[-sin(nu); e+cos(nu); 0]; % % now rotate % cO = cos(Om); sO = sin(Om); co = cos(om); so = sin(om); ci = cos(i); si = sin(i); R = [cO*co-sO*so*ci -cO*so-sO*co*ci sO*si; sO*co+cO*so*ci -sO*so+cO*co*ci -cO*si; so*si co*si ci]; ri = (R*rv)'; vi = (R*vv)'; 18 prop_system.m % Bill Nadir % 16.851 Satellite Engineering % 10/11/2003 % Module for calculating spacecraft propulsion system mass for required % momentum dumping function [F, p_mass] = prop_system(size,CG,lifetime,H,sat_rate) % % % % % % % % % % % % % INPUTS size CG lifetime H sat_rate OUTPUTS p_mass F = edge lengths of the S/C: vector (x,y,z) side lengths (meters) = (x,y,z) coordinates of the location of the CG (offset from the geoometric center of the S/C) = required lifetime of the spacecraft [yrs] = maximum stored momentum in any one momentum wheel (saturation point of a momentum wheel) [N*m*s] = The rate of saturation of a momentum wheel (used to determine how often momentum needs to be dumped) [days/saturation] = total mass of the propulsion subsystem which will provide momentum dumping capability for the spacecraft [kg] = Thrust required to dump momentum % Hydrazine (monopropellant) is chosen as the fuel for this propulsion % system and a conservative specific impulse, Isp, is 200 seconds Isp = 200; % Here the earth's gravity constant is initialized (9.8 m/s^2) g = 9.8; % % % t Here the impulse time, t, of the thruster firing is set It is assumed that the thruster required for momentum dumping will fire for 1 second = 1; % Here the locations of the six required thrusters are initialized [x y z] % Each row is for a different thruster thruster = [0 size(2)/2 size(3)/2; 0 size(2)/2 -size(3)/2; size(1)/2 0 size(3)/2; -size(1)/2 0 size(3)/2; size(1)/2 -size(2)/2 0; -size(1)/2 -size(2)/2 0]; % Here the moment arms for the six thrusters from the CG are determined % For X-thrusters (spin about X-axis), moment arm is in Y-direction (cols 1,2) % For Y-thrusters (spin about Y-axis), moment arm is in Z-direction (cols 3,4) % For Z-thrusters (spin about Z-axis), moment arm is in X-direction (cols 5,6) moment_arms = [abs(CG(2) - thruster(1,2)) abs(CG(2)- thruster(2,2)) abs(CG(3) thruster(3,3)) abs(CG(3) - thruster(4,3)) abs(CG(1) - thruster(5,1)) abs(CG(1) - thruster(6,1))]; % Here we will assume the worst-case distance from the thruster to the CG % (shortest) which will require the largest thrust to impart the required % torque on the S/C for momentum dumping worst_moment_arm = min(moment_arms); % Here the thrust required to dump the momentum is calculated (per pulse) F = H / (worst_moment_arm * t); % Here the required propellant mass for this propulsion system is estimated total_pulses = (lifetime * 365.25) / sat_rate; % total thruster pulses required over lifetime m_prop = (F * total_pulses * t)/(Isp * g); % mass in kg % Here the total propulsion system mass is determined by assuming that 85% % of the propulsion system mass is propellant (SMAD, p. 660) - conservative p_mass = m_prop / 0.85; % mass in kg 19 material_prop.xls20 NAME OpticalSolarReflector QuartzOverSilver SilvercoatedFEP SilveredTeflon AluminizedTeflon WhiteEpoxy WhiteEnamel AluminizedFEP SilverPaint SolarCellFusedSilica BlackPaint Titanium6AL4V SteelAm350 Titanium6AL4V AluminiumTape Aluminum606T6 Gold Aluminum606T6 GoldizedKapton PolishedBeryllium Info SSE SMAD SSE SMAD SMAD Al.Substrate Al.Substrate SSE SSE SMAD Al.Substrate AsReceived AsReceived Polished SSE Polished AsRolled AsReceived SSE SSE MaterialNumber 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 20 Absorptivity 0.07 0.077 0.08 0.08 0.163 0.248 0.252 0.16 0.37 0.805 0.975 0.766 0.567 0.448 0.21 0.2 0.299 0.379 0.25 0.44 Reflectance 0.93 0.923 0.92 0.92 0.837 0.752 0.748 0.84 0.63 0.195 0.025 0.234 0.433 0.552 0.79 0.8 0.701 0.621 0.75 0.56 References 1 Meriam, J.L., and Kraige, L.G., Engineering Mechanic: Dynamic, 4th Ed, John Wiley, Inc. 1997. Hall, Christopher D., oe2rv.m software tool, http://www.aoe.vt.edu/~cdhall/ 3 Benson, Tom, http://www.grc.nasa.gov/WWW/K-12/airplane/atmosmet.html, Earth Atmosphere Model, Metric Units, NASA Glenn Research Center, 2002. 4 Chobotov, Vladimir A. (Editor), Orbital Mechanics, 2nd Ed., AIAA, 1996, p. 227 5 Benson, Tom, http://www.grc.nasa.gov/WWW/K-12/airplane/atmosmet.html, Earth Atmosphere Model, Metric Units, NASA Glenn Research Center, 2002. 6 Wertz, James, and Larson, Wiley, Space Mission Analysis and Design, 2nd Ed., Microcosm, Inc., 1997, p. 353. 7 E. Ahronovich, M. Balling, Reaction Wheel and Drive Electronics For LeoStar Class Space Vehicles, 12th Annual USU Conference on Small Satellites, 1998, www.sdl.usu.edu/conferences/smallsat/proceedings/12/ssc98/1/ssci5.pdf 8 Dynacon Enterprises Limited, Dynacon MicroWheel 200, www.dynacon.ca/pdf/files/productpdf_6.pdf 9 Honeywell Aerospace Electronic Systems, Constellation Series Reaction Wheels, http://content.honeywell.com/dses/assets/datasheets/constellation_series_reaction_wheels.pdf. 10 Honeywell Aerospace Electronic Systems, Miniature Reaction Wheels, http://content.honeywell.com/dses/assets/datasheets/mini-wheel_reaction_wheel.pdf 11 Honeywell Aerospace Electronic Systems, Honeywell Model HR 0610 Reaction Wheel, http://content.honeywell.com/dses/assets/datasheets/hr0610_reaction_wheel.pdf 12 Teldix Space Product Group, Momentum and Reaction Wheels 14-68 Nms with external Wheel Drive Electronics, http://www.teldix.de/P22/RDR23-68.pdf 13 Teldix Space Product Group, Momentum and Reaction Wheels 14-68 Nms with integrated Wheel Drive Electronics, http://www.teldix.de/P22/RSI25-68.pdf 14 Teldix Space Product Group, Momentum and Reaction Wheels 2-12 Nms with integrated Wheel Drive Electronics, http://www.teldix.de/P22/RSI4-12.pdf 15 Teldix Space Product Group, High motor torque Momentum and Reaction Wheels 14-68 Nms with integrated Wheel Drive Electronics, http://www.teldix.de/P22/RSI18-68.pdf 16 Teldix Space Product Group, Momentum and Reaction Wheels 0.04-0.12 Nms with integrated Wheel Drive Electronics, http://www.teldix.de/P22/RSI01.pdf 2 Teldix Space Product Group, Momentum and Reaction Wheels 0.2-1.6 Nms with integrated Wheel Drive Electronics, http://www.teldix.de/P22/RSI02.pdf 17 18 Sellers, Jerry Jon, Understanding Space: An Introduction to Astronautics, 2nd Ed., McGraw Hill, 2000, p. 570. Wertz, James, and Larson, Wiley, Space Mission Analysis and Design, 2nd Ed., Microcosm, Inc., 1997, p. 660. 20 Fortescue, Peter, Stark, John, and Swinerd, Graham, Spacecraft Systems Engineering. John Wiley and Sons Ltd., Third Edition. 19 21