Finding the Force -- Consistent Particle Seeding for Satellite Aerodynamics
FFinding the Force—Consistent Particle Seeding forSatellite Aerodynamics
J. Brent Parham ∗ MIT Lincoln Laboratory, Lexington, MA 02421, U.S.A.
L. A. Barba † Mechanical Engineering, Boston University, Boston, MA 02215, U.S.A.
When calculating satellite trajectories in low-earth orbit, engineers need to adequatelyestimate aerodynamic forces. But to this day, obtaining the drag acting on the complicatedshapes of modern spacecraft suffers from many sources of error. While part of the problemis the uncertain density in the upper atmosphere, this works focuses on improving the mod-eling of interacting rarified gases and satellite surfaces. The only numerical approach thatcurrently captures effects in this flow regime—like self-shadowing and multiple molecularreflections—is known as test-particle Monte Carlo. This method executes a ray-tracingalgorithm to follow particles that pass through a control volume containing the spacecraftand accumulates the momentum transfer to the body surfaces. Statistical fluctuations in-herent in the approach demand particle numbers in the order of millions, often making thisscheme too costly to be practical. This work presents a parallel test-particle Monte Carlomethod that takes advantage of both GPUs and multi-core CPUs. The speed at which thismodel can run with millions of particles allowed exploring a regime where a flaw in themodels initial particle seeding was revealed. Our new model introduces an analytical fixbased on seeding the calculation with an initial distribution of particles at the boundary ofa spherical control volume and computing the integral for the correct number flux. Thiswork includes verification of the proposed model using analytical solutions for several sim-ple geometries and demonstrates uses for studying aero-stabilization of the Phobos-GruntMartian probe and pose-estimation for the ICESat mission.
Nomenclature k B = Boltzmann constant, J/K T = Temperature, K m = Molecular mass of gas particles, kg A ref = Reference Area, m n = Number density, ρ = mass density ( nm ), kg/m C D = Drag coefficient, (Drag Force)/ (cid:0) ρV A ref (cid:1) Γ = Number flux through surface c mp = Most probable thermal speed, (cid:112) k B T ∞ /mS = Speed ratio, V /c mp Subscripts ∞ = Free stream equilibrium conditions W = Properties evaluated at the satellite surface ∗ Lincoln Scholar at Boston University, Department of Mechanical Engineering; AIAA Member. † New address: Mechanical and Aerospace Engineering, George Washington University; AIAA Member.
Distribution Statement A. Approved for public release; distribution is unlimited. a r X i v : . [ phy s i c s . c o m p - ph ] D ec . Introduction N early five years ago , the active communications satellite Iridium 33 crashed into Cosmos 2251—a defunct Russian satellite—scattering thousands of satellite fragments into low-earth orbit. a Thismajor space collision, among other similar events, alerted satellite operators to the dangers of space-debrisproliferation and motivated world-wide response [1]. With a sharp increase in new objects to track in space,accurately estimating the orbits and trajectories of satellites while reducing uncertainties is now a criticalaspect of preventing future collisions and reducing the cost of maneuvers to evade debris.Just a year before the Iridium 33 crash, Vallado [2] had pointed out that a lack of research in the area ofphysically consistent drag estimation was an obstacle for predicting satellite orbits accurately. The situationremains the same today. Aerodynamic forces induce the most uncertain accelerations encountered by asatellite in Low-Earth Orbit (LEO), but they remain elusive to calculate. The drag coefficient itself is ageometry-dependent factor multiplying the dynamic pressure to give the drag force, as shown in equation(1), with atmospheric neutral density ρ , and orbital velocity V . F Drag = C D A ref ρV (1)The main uncertainties in this equation arise with specifying the neutral density and calculating the dragarea ( C D A ref ) for a given satellite, which can have an arbitrarily complicated geometry. While recent workhas improved the accuracy of atmospheric neutral density models [3, 4] by inferring densities from satellitetracking data—e.g., the HASDM model [5]—, calculating satellite coefficients numerically in the context oforbit determination needs to be revisited and standardized. Uncertainty reduction in all terms of the dragforce is the only way to increase the precision of its estimate.The earliest attempts to model satellite aerodynamics were analytical [6, 7], with aerodynamic forcesdirectly integrated from momentum flux over a surface. But closed-form solutions were only obtained forsimple geometries: flat plate, cone, cylinder and sphere. These solutions omitted effects such as shadowingand multiple reflections, which may change aerodynamic properties. As computational power grew throughthe years, the Test-Particle Monte Carlo (TMPC) method [8, 9, 10] appeared as an alternative to analyticalmodels. Test-particle methods trace simulated particles through a control volume that contains the space-craft. Each particle can then interact with the surface to transfer momentum and energy through a surfaceinteraction model (e.g., diffuse, specular, etc.). The procedure reduces to a Monte Carlo initialization and ageometric ray-tracing algorithm. TPMC methods remain the only way to adequately account for shadowingand multiple reflections with arbitrary surface interaction. But, like all Monte Carlo methods, TPMC needsa large number of particles to reduce the statistical scatter of the solutions. The computational cost thereforeseverely limits its application to orbital analysis.A TPMC simulation is set up by placing particles randomly with a uniform distribution on the faces ofa rectangular control volume. This initialization, however, is inconsistent with kinetic theory for arbitraryfree-stream incidence. Here, we construct a new method to initialize a TPMC simulation that solves this in-consistency. By choosing a spherical control volume, we reveal an exact solution for initial-position likelihoodof incident particles and the total number flux into the volume. We developed a code implementing TPMC,then extended it to include the new method and to run in parallel using modern computational hardware—both multi-core CPUs and many-core GPUs. In this paper, we describe the mathematical derivation for thephysically consistent test-particle seeding along with the validation of the method and verification of thecode using analytical tests. We also demonstrate the application of the new method to orbital analysis withactual satellite data. a A software reconstruction of the collision is available on video at http://youtu.be/_o7EKlqCE20 , courtesy of AGI.2 of 14American Institute of Aeronautics and Astronautics
I. Methods
II.A. Test-particle Monte Carlo (TPMC) method
A rarefied-gas flow is well modeled by a drifting Maxwellian gas, i.e., a gas described by the following velocitydistribution, f ( v i ) = (cid:114) m πk B T exp (cid:20) − m ( v i − u i ) k B T (cid:21) , (2)where the bulk average velocity components, u i , are superimposed on the individual distribution components, v i and an isotropic particle density distribution in space is assumed, f ( x i ) = constant.In free molecular flow—where the mean free path of molecules is large in relation to the volume considered—the gas particles do not interact with each other in the free stream nor with the particle distribution createdby surface reflections. The TPMC method takes advantage of this limit by ignoring particle-particle in-teractions and treating particles independently. It then obtains the average momentum transfer to a bodyfrom the sum of discrete particle reflections. “Particles” here are meant to represent a large number ofgas molecules, for which a weighting factor is introduced. Using test-particle methods requires accuratesampling of the free-stream flux into the specified control volume and a model of particle interactions withbody surfaces.Fritsche and Klinkrad [11] explain a flux-sampling scheme of the TPMC method implemented in thesoftware that is used by the European Space Agency, called ANGARA. The sampling consists of placingparticles with a random uniform distribution on the six faces of a rectangular control volume. The numberof particles to cross each face derives from the analytical molecular flux across a flat plate with given normalvelocity component of the free stream and the total number of particles in the TPMC simulation (a parameterdictating accuracy). Once the particles are initialized on the control surface, they are followed via a ray-tracing algorithm until they intersect with the surface of the spacecraft. Satellite surfaces are represented bya triangulation, allowing the use of well-known ray-tracing methods. To account for shadowing, the closestintersection point to the initial position becomes the first collision with the satellite and the particle is movedto this point. At the new location, a surface-interaction model determines the new velocity for the particle.The surface-interaction model most used in satellite aerodynamics is diffuse reflection, which is bothsimple and has been proven to adequately model satellite surface interactions in LEO [12]. In this model,the incident particles settle in the surface material and are reemitted with a velocity distribution determinedby the satellite surface temperature. In TPMC, a random departure direction from the surface is specified andthe reemission speed is chosen from the Maxwellian speed distribution (the product of the three componentsin equation (2) integrated over all angles with spherical coordinates in velocity space). The differencebetween the original and final velocities are then scaled by the particle’s effective mass and added to thetotal momentum exchange with the spacecraft. To calculate the particle’s effective mass, the total number ofincident molecules, derived from the molecular flux, is divided by the number of simulated particles (specifiedby the user) and multiplied by the molecular mass. The final resultant aerodynamic force is then the totalmomentum exchange that includes contributions from all particles. II.B. Analytical solutions for simple geometries
The following three analytical solutions provide a way to scrutinize the current model. The solutions assumea diffuse reemission and the resultant force derives from a surface integral of the momentum flux on thesatellite surface. The influx of momentum is integrated using the free-stream distribution and the outfluxis calculated with the reemission distribution described by the satellite temperature. To simplify, we definethe ratio of the free-stream speed magnitude, U ∞ , to the most-probable thermal speed, c mp = (cid:112) k B T ∞ /m ,as the “speed ratio,” S = U ∞ /c mp . II.B.1. Flat plate at angle of incidence
This solution apears in Sentman’s work.[6] The first three terms represent the influx of momentum integratedover a planar surface and the last term is the reemission contribution to the total momentum flux. The drag oefficient is then given by C D = 2 √ πS (cid:18) exp[ − ( S cos α ) ] + √ πS cos α (cid:18) S (cid:19) erf[ S cos α ] + πSS W cos α (cid:19) , (3)where S W is the speed ratio using the wall temperature of the body in c mp and α is the angle of attack. II.B.2. Convex spherical segment
As given by Bird [9], this result is recreated by integrating equation (3) over a spherical surface with theappropriate change in angle of attack and choice of differential area element. The drag coefficient is C D = exp[ − S ] √ πS (1 + 2 S ) + 4 S + 4 S − S erf[ S ] + 2 √ π S W . (4) II.B.3. Concave spherical segment in hyperthermal flow
Pratts’s thesis on concave bodies in hyperthermal flow [13] generalizes the drag coefficient of a hemisphericalbody with full thermal accommodation to the following expression (in his notation): C D = 2 + (cid:15) D ( θ ) √ πS (cid:114) T W T ∞ . (5)Here, (cid:15) D ( θ = 90 ◦ ) = 1 . II.C. Satellite-tracking derived coefficients
The accelerations caused by drag forces can manifest themselves in tracking data of satellites. By trackingsatellites for a length of time, the orbital state and ballistic coefficient are derived with an orbit-determinationroutine, most commonly a least-squares optimization. The North American Aerospace Defense Command(NORAD) releases the resulting states publicly through the web service Space-Track. b Vallado standardizedthe details of the Simplified General Perturbations (SGP4) dynamical model used to compute these two-lineelement sets (TLEs) with his interpretation of the original specifications [14]. We use his interpretation forthe analysis of TLEs.For an LEO object, the TLE contains information about the drag force during the period that the elementset was created. This ballistic coefficient B ∗ of the TLE maps directly to the the drag area ( R ⊕ = 6378 . ρ = 2 . × − kg/m ): C D A ref = 2 B ∗ M satellite R ⊕ ρ (6)The International Laser Ranging Service (ILRS) c also releases tracking data for participating satellites.An orbit determination routine—along with an initial orbit estimate from a TLE—creates a precise orbitalstate with this range data. The precision estimate serves as a check to the TLE data for the ICESat missionand is processed with an orbit determination routine based on the work of Montenbruck and Gill [15]. III. Results
III.A. Physically consistent particle seeding for TPMC
TPMC simulations begin by placing a uniform random distribution of particles on the faces of a rectangulardomain. The increasing demands—i.e., accuracy for orbital analysis—expose the biases introduced by seedingon a rectangular domain: the method constructs spurious high-density regions at the edges of the domain b c http://ilrs.gsfc.nasa.gov/data_and_products igure 1: Percent error of a 10-million-particle TPMC calculation with respect to exact sphere drag. Bothmodels assume fully diffuse reflection. The TPMC result is not a constant function of incidence angle to therectangular control volume.faces. As a result, the drag coefficient (in this test, of the sphere) shows an artificial increase as the freestream is swept around the rectangular control volume (see Figure 1).To avoid these density artifacts, we modified the TPMC method by developing a particle-seeding schemeon a spherical control volume. We obtained closed-form solutions for this control volume by solving theassociated probability integrals for a drifting Maxwell gas. The solution gives the correct number flux andparticle distribution on the spherical shell, which in turn gives the correct particle weights, sampled positionsand sampled velocities. III.A.1. Position sampling
The number flux and particle distribution are set up as an integral of the velocity distribution along thespherical surface. We start by defining a flux, Γ, of particles into a surface, with a given velocity distribution f ( v ), and write the following expression (where n is the number density):Γ = n (cid:90) Velocity (cid:20)(cid:90)
Area v f ( v ) · d a (cid:21) d v . (7)We consider a set of orthogonal coordinates that are aligned with the normal and tangent vectors, asshown in Figure 2 (from here onwards ˆ t = ˆ t and ˆ t = ˆ t × ˆ n ). Aligning the inward velocity with the sphericalnormal, we rewrite the integral as shown in equation (8). The bulk velocity components u, v, and w areexpressed in the local coordinate system.Γ = n ( c mp π ) / (cid:90) Area (cid:20)(cid:90) ∞ (cid:90) ∞−∞ (cid:90) ∞−∞ v n exp (cid:18) − ( v n − u ) + ( v t − v ) + ( v t − w ) c mp (cid:19) dv t dv t dv n (cid:21) da = n ( c mp π ) / (cid:90) π/ φ = − π/ (cid:90) πθ =0 c mp π (cid:18) c mp e − u c mp + √ πuc mp (cid:20) (cid:18) uc mp (cid:19)(cid:21)(cid:19) r cos φdθdφ (8)The free stream is aligned with the { φ, θ } = { π/ , } direction so that only one component remains inspecified spherical coordinates ( u = v sin φ ), and we can integrate to find the result: = nr ( c mp π ) / π c mp e − v c mp + π / c mp (cid:0) c mp + 2 v (cid:1) erf (cid:16) vc mp (cid:17) v . (9)Simplifying further, by assuming that both v and c mp are positive and defining the speed ratio as S = v/c mp , we obtain Γ = (cid:34) √ πc mp e − v c mp + (cid:32) πc mp v + πv (cid:33) erf (cid:18) vc mp (cid:19)(cid:35) nr = (cid:104) √ πe − S + (cid:16) π S + πS (cid:17) erf ( S ) (cid:105) c mp nr . (10)This number flux is an exact expression for the flux into a sphere for all Maxwellian gases in equilibrium.With this result, a likelihood g ( φ | S ) of incidence angle φ ∈ [ − π/ , π/
2] into a sphere can be constructed andcharacterized only by the speed ratio, as follows: g ( φ | S ) = cos φ (cid:104) e − S sin φ + S √ π sin φ (1 + erf( S sin φ )) (cid:105) e − S + √ π (1+2 S )erf( S )2 S . (11)Figure 2: Spherical geometry used to solvemolecular influx to the control volume. The evaluation of this distribution for several speed ratiosis plotted in Figure 3a. For low speed ratio the distributionapproaches the limit lim S → g ( φ | S ) = 12 cos φ. (12)When the density in equation (12) is combined with auniform distribution sampling θ ∈ [0 , π ], we obtain a uni-form sampling on a spherical surface. This limit exhibits thespatially isotropic thermal influxes from a resting gas that isin equilibrium. With large speed ratio, the distribution ap-proaches lim S →∞ g ( φ | S ∞ ) = sin(2 φ ) φ ≥ φ < φ and θ ) exist, we have succeeded and can write the position of the particle in Cartesianspace. III.A.2. Velocity sampling
With a position sampled, we choose the velocity using the canonical velocity component distributions for fluxpast a flat plate, i.e., a surface element with the normal aligned with the position vector. The distributionsin velocity space transform into the local coordinate frame with respect to the sphere’s normal vector, and
80 −60 −40 −20 0 20 40 60 8000.20.40.60.81 φ [deg] P r o b a b ili t y D e n s i t y S = 0.1S = 1S = 10 (a) Likelihood for distribution of particles based on angle of incidence andvarious speed ratios.
S=0.1 S=1 S=10
Flow Direction (b) Initial position sampling of 10,000 particles for various speed ratios
Figure 3: Physically consistent particle sampling, applied to a spherical volume.are represented in Equations (14a), (14b) and (14c). p ( v n ) = 2 v i exp[ − (cid:16) v n c mp − S n (cid:17) ]exp[ − S n ] + √ πS n (1 + erf[ S n ]) (14a) p ( v t ) = 1 √ π exp (cid:34) − (cid:18) v t c mp − S t (cid:19) (cid:35) (14b) p ( v t ) = 1 √ π exp (cid:34) − (cid:18) v t c mp − S t (cid:19) (cid:35) (14c)The free stream velocity component-wise speed-ratio in the local coordinate system can be written as, S t = | S | cos φ , and S n = | S | sin φ (inwardly pointing here) leaving the third component, S t , always equalto zero. III.A.3. Rotation to body coordinates
The position vectors are rotated from the free-stream centered frame, with axis ˜ z , back into the satellitereference frame, with axis ˆ z , with a rotation about the vector axis ˆ u = ˜ z × ˆ z by the angle β = cos − (˜ z · ˆ z ).This rotation is implemented using Rodrigues’ rotation formula, given as: x (cid:48) = x cos β + (ˆ u × x ) sin β + ˆ u (ˆ u · x )(1 − cos β ) (15)The velocity components transform back to the satellite reference frame via multiplication by the transposeof the local coordinate matrix, consisting of the vectors ˆ n, ˆ t and ˆ t . N particles S p ee d u p [ C P U ] / [ G P U ] (a) NVIDIA K20 GPU compared to Intel Core i7 CPU −3 −2 −1 N particles σ C D / h C D i Simulation1 / p N particles (b) Decreasing statistical scatter with more particles Figure 4: Simulation workload and convergence characterizations
III.B. Software implementation and characterization
We implemented the revised TPMC ray-tracing method in CUDA—parallelized over particles—and testedon a Kepler-architecture NVIDIA Tesla K20 GPU. We also wrote an OpenMP version to run efficiently onan Intel Core i7 CPU. Chip-to-chip comparison shows up to 35-times speedup on the GPU, with appropriateproblem size (Figure 4a). For all calculations in the following results, we ran the CUDA version with aparticle count of 10 million (providing a statistical scatter of approximately 0.07% for the comparison withthe flat-plate analytical solution, as shown in Figure 4b).
III.C. Verification with analytical solutions
The comparison with analytical solutions allows a strict verification of the new boundary conditions and codeimplementations against free-molecular theory. For the cases shown, we set the free-stream temperature to T ∞ = 922 K and the body-surface temperature to T W = 300 K. These conditions are typical of thoseexperienced at an altitude of about 400 km. The triangular surface mesh of the simple geometric shapeswere created using Gmsh [16].Figure 5 shows that when the code is compared to the analytical solution of a double-sided flat plate(i.e., a plate that can absorb free-stream molecules from the front and back), our model exactly capturesthe variation with the speed ratio from near-thermal flows speeds to the extremely fast flows encounteredby spacecraft. As angle of attack changes from normal to the plate ( α = 0 ◦ ) to parallel, C D decreasesas expected. The thermal fluctuations near the conditions at altitude ( S ∞ = 7) create a non-zero drageven on surfaces that are parallel to the flow. This effect often causes underestimates in drag forces whenthe hyper-thermal approximation is used for analysis. The full TPMC method accurately reproduces thisphysical phenomenon and can readily be extended to more complicated bodies.To verify the method against a slightly more complex body, we turn to the analytical solution for asphere, as seen in Figure 6. Again, the physical variation with the speed ratio is captured by the TPMCmethod and we can confidently say that the code reproduces forces on convex surfaces. As a check to see ifthe code can accurately model the effects of multiple reflections, we compare Pratt’s hyper-thermal modelto the TPMC calculation of a concave shell over the whole domain of speed ratio (also in Figure 6). TheTPMC solution asymptotically approaches the hyper-thermal solution as S ∞ increases. Although there isno direct solution of the reflection integrals for the whole domain, we can reasonably draw a conclusion thatthis verifies the reflection algorithm implemented in the method. Free Stream Speed Ratio, S ∞ D r a g C o e ff i c i e n t , C D SimulationTheory = 0= 30= 60= 90
Figure 5: The analytical solution of a double-sided flat plate at angle of attack compares favorably to theoutput of the TPMC code developed. In this figure, free stream speed ratio and angle of attack are variedto show the dependance of the drag coefficient on those variables. S ∞ C D Analytical SphereHyperthermal Concave Shell − Pratt 1963Sphere SimulationConcave Shell Simulation
Figure 6: The simulation for a sphere with varying speed ratio compared to its analytical solution. Thesolution of a concave spherical shell is also shown, and it asymptotes to the hyper-thermal approximationgiven by Pratt.
II.D. Demonstration with satellite-tracking data
III.D.1. Phobos-Grunt aero-stabilization
Nov. 15 Dec. 1 Dec. 15 Jan. 1 Jan. 1510152025303540 050100150200250
Mean Altitude [km] C D A ref [m ] Figure 7: The Phobos-Grunt probe aero-stabilized by themiddle of December 2011, and Two Line Element (TLE) setderived force coefficients allow for comparison of model es-timates to actual dynamical characteristics during this sta-bilization.Phobos-Grunt, a multinational scientific mis-sion that failed to escape its initial parkingorbit, reentered in early January 2012. Thisfailure provided a unique opportunity to ob-serve a large spacecraft’s uncontrolled interac-tions with the near-Earth environment. Fromits launch in November 2011 to its reentry,the craft was observed by many amateur as-tronomers d and appeared to aero-stabilize inmid-December 2011. This behavior also ap-pears in TLE-derived force coefficients (Figure7) and can be compared to an estimated modelof the craft’s geometry. Using a model basedon Phobos-Grunt, we can simulate the torquesencountered at altitude—with conditions spec-ified in Table 1.All moments are shifted to a “quarter-chord” position that is measured from the fueltank end. The total length of the craft, c ≈ . T ∞ T W
300 K V ρ ∞ . × − kg/m q ∞ × − N/m To understand the stability of the spacecraft, we vary the angle of incidence about two axes and calculatethe moment about the axis that would induce pitch; see Figure 8. This calculation reveals that for increasingangle α there is an increasing moment that serves to decrease the magnitude of α . As shown in the figure,this characteristic moment is roughly symmetric and produces a stability point very close to α = 0 ◦ . Thereis however asymmetry as we rotate about the roll angle β .To better understand the asymmetry we look at the spacecraft geometry. Choosing two β -angles sothat the line connecting the center of each solar panel is either “Horizontal” (parallel to the pitch axis) or“Vertical” (perpendicular to the pitch axis), we recast the moments in terms of a pitching moment coefficient.To do this, we divide the moment by the dynamic pressure q ∞ and the chord length c , leaving the coefficientscaled by an unknown reference area. Defining positive moments as anything that pitches the spacecraftup, we obtain the results in Figure 9. The asymmetry of the moments appears at high angle of attack, andmanifests itself as a decreased moment for the “Vertical” configuration. We hypothesize that this decreasein moment coefficient is due to a shadowing of a solar panel by the tank end of the spacecraft as it increasesangle of attack. If the panel was not shadowed, the dominant role of the drag force (lacking a substantial lift d http://legault.perso.sfr.fr/phobos-grunt.html
10 of 14American Institute of Aeronautics and Astronautics (cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:2)(cid:8)(cid:9)(cid:4)(cid:10)(cid:11)(cid:6)(cid:5)(cid:12)(cid:13)(cid:2)(cid:8)(cid:14)(cid:15)(cid:8)(cid:1)(cid:2)(cid:16)(cid:5)(cid:14)(cid:17)(cid:6)(cid:11)(cid:10)(cid:8)(cid:9)(cid:14)(cid:18)(cid:2)(cid:11)(cid:5)
Figure 8: The restoring moment—aerodynamic moment about c/ c = 5 . α —for the spacecraft model. The moment is not axisymmetric about the rollaxis, but it is monotonically increasing with increasing α .force) in creating the free-molecular moment would mean the moment coefficient would be more axisymmetricabout the roll axis. This phenomenon therefore can only appear when complex bodies are simulated withthe appropriate considerations for flow shadowing. III.D.2. ICESat pose estimation
Figure 10: Surface discretization used forICESat simulations. Solar panels are ar-ticulated. For the Airplane mode calcu-lations, we use the average of calculationswith the panels perpendicular and parallelto the body axis.The ICESat science mission operated in two attitude modesduring its mission: “airplane” and “sailboat.” These twomodes, chosen to best present the solar panels to the sun,can clearly appear in derived force coefficients from the In-ternational Laser Ranging Service (ILRS) tracking data. Toshow how this discrete attitude change is filtered through thedrag-coefficient estimation with a typical orbit-determinationroutine, we create a model of ICESat (Figure 10) and imposedthe attitude constraints based on Webb’s work [17].The orbit, along with the drag area C D A ref , is estimatedwith a least-squares filter that uses a three-day sliding windowof laser ranges for its calculation. The atmospheric conditionssampled from NRLMSISE-00 [18] and orbital velocity at themiddle of each span are fed into the TPMC method to modela free-molecular drag coefficient with the TPMC method. InFigure 11, we compare the estimated drag area to the modeleddrag area. We notice that the free-molecular drag area does notappreciably vary over the course of the year,despite changes intemperature of a few hundred Kelvin and density fluctuations by an order of magnitude. The tracking dataestimates, however, have a slow secular drift during “sailboat” mode and a wide scatter during the “airplane”portion of flight. The simulation does allow us to precisely understand the step response of drag area filteredthrough tracking data, but the ILRS data only provides enough precision to roughly understand the physicsof the aerodynamic forces on such a complicated body.
11 of 14American Institute of Aeronautics and Astronautics
40 −20 0 20 40−15−10−5051015 Angle of Attack, α [deg] C M A r e f VerticalHorizontal
Figure 9: Calculation of a pitching moment coefficient about the quarter-chord ( c = 5 . C D A r e f Day of year, 2009 − ρ m a ss [ k g / m ] Day of year, 2009
ILRS Data DerivedSimulated "Sailboat" "Airplane"
Figure 11: ICESat C D as a function of the time in 2009. The drag coefficient is simulated using the on orbitvelocities and atmospheric properties sampled from the NRLMSISE-00 model, and the second panel showsthe atmospheric density from that model.
12 of 14American Institute of Aeronautics and Astronautics
V. Conclusion
By demonstrating an improved approach to free-molecular aerodynamics, we hope to propel a stagnantfield forward. The new physically consistent boundary condition ensures confidence in the calculation of forcesand moments for a satellite model and with the power of modern workstation hardware, calculations becomeusable in a wider range of orbital analysis. Beyond that, we provide verification tests (even accounting forconcave bodies) for the TPMC method—often overlooked previously—and explore two intriguing applicationsin satellite dynamics.With this tool, we show that the uncontrolled dynamics of a complicated surface in the free-molecularregime may be dependent on accurate modeling of self-shadowing geometry, and prove aerodynamic stabilityfor a spacecraft based on the Phobos-Grunt probe. We can also extend the analysis drag forces with trackingdata to complicated bodies with known geometry and attitude, such as ICESat. Although attitude andgeometry are known for this spacecraft, the data provided by the ILRS service—some of the most precisetracking data available—when mixed with the NRLMSISE-00 ICESat atmospheric model is only accurateenough to qualitatively compare drag-force estimates to the physical model. This is enough, however, tounderstand large attitude changes with orbit determination methods based on metric tracking data.To summarize, we have shown the following: • The TPMC method with its application to rarefied gas dynamics is ideal for implementation on newmassively parallel high-bandwidth coprocessors such as GPUs. • We can derive set of boundary conditions that describe the correct distribution of states for influxparticles that enter a spherical control volume. • TPMC codes can be verified appropriately using available analytical solutions. • TPMC codes can be used to calculate aerodynamic stability of satellites in the free-molecular regime. • We obtain qualitative agreement between free-molecular aerodynamic methods and real-world satellitedynamics with the ICESat mission.The analysis and tool resulting from this work promises that future work may include physically correctdynamical simulation for attitude and orbital-state propagation. As a result, finding the force for satellitedrag becomes a little less uncertain.
Acknowledgments
This work is sponsored by the U.S. Department of the Air Force under Air Force Contract FA8721-05-C-0002. Opinions, interpretations, conclusions and recommendations are those of the authors and are notnecessarily endorsed by the United States Government. The first author wishes to express special thanks toDavid Gonzales for many years of support and insightful discussions about satellite dynamics and rarefiedgases. He also extends his thanks to MIT Lincoln Laboratory for supporting him through the LincolnScholars Program.
References [1] United Nations Office for Outer Space Affairs, “Space Debris Mitigation Guidlines of the Committee on the PeacefulUses of Outer Space,” Vienna, 2010.[2] Vallado, D. A. and Finkleman, D., “A Critical Assessment of Satellite Drag and Atmospheric Density Modeling,”
AIAA/AAS Astrodynamics Specialist Conference , Honolulu, Hawaii, August 2008. doi:10.2514/6.2008-6442.[3] Sutton, E. K., “Normalized Force Coefficients for Satellites with Elongated Shapes,”
Journal of Spacecraft and Rockets ,Vol. 46, January 2009, pp. 112–116. doi:10.2514/1.40940.[4] Marcos, F. A., Lai, S. T., Huang, C. Y., Lin, C. S., Retterer, J. M., Delay, S. H., and Sutton, E., “Towards Next LevelSatellite Drag Modeling,”
AIAA Atmospheric and Space Environments Conference , Toronto, Ontario Canada, August2010. doi:10.2514/6.2010-7840.[5] Storz, M. F., Bowman, B. R., Branson, M. J. I., Casali, S. J., and Tobiska, W. K., “High Accuracy Satellite Drag Model(HASDM),”
Advances in Space Research , Vol. 36, 2005, pp. 2497–2505. doi:10.1016/j.asr.2004.02.020.[6] Sentman, L. H., “Free Molecule Flow Theory and Its Application to the Determination of Aerodynamic Forces,” Tech.Rep. 448514, Lockheed Missiles and Space Company, Sunnyvale, California, October 1961.13 of 14American Institute of Aeronautics and Astronautics7] Schaaf, S. A. and Chambre, P. A.,
Flow of Rarefied Gases , Princeton Aeronautical Paperbacks, Princeton UniversityPress, 1961.[8] Bogacheva, A. A., Perepukhov, V. A., and Rukhman, E. Y., “Monte Carlo Computation of Aerodynamic Characteristicsof Bodies of Complex Shape in Free Molecule Flow,” Technical Translation TT F-12-360, NASA, July 1969.[9] Bird, G. A.,
Molecular Gas Dynamics and the Direct Simulation of Gas Flows , Oxford University Press, 1994.[10] Klinkrad, H., KoppenwaIlner, G., Johannsmeier, D., Ivanov, M., and Kashkovsky, A., “Free-Molecular and TransitionalAerodynamics of Spacecraft,”
Advanced Space Research , Vol. 16, No. 12, 1995, pp. (12)33–(12)36.[11] Fritsche, B. and Klinkrad, H., “Accurate Prediction of Non-Gravitational Forces for Precise Orbit Determination - PartI: Principles of the Computation of Coefficients of Force and Torque,”
AIAA/AAS Astrodynamics Specialist Conferenceand Exhibit , American Institute of Aeronautics and Astronautics, 2013/05/13 2004. doi:10.2514/6.2004-5461.[12] Moe, K. and Moe, M. M., “Gas-Surface Interactions in Low-Earth Orbit,”
AIP Conference Proceedings , Vol. 1333, 2011,p. 1313.[13] Pratt, M.,
The Free-Molecule Characteristics of Concave Surfaces , Ph.D. thesis, Cranfield University, 1963.[14] Vallado, D., Crawford, P., Hujsak, R., and Kelso, T. S., “Revisiting Spacetrack Report
AIAA/AAS Astro-dynamics Specialist Conference and Exhibit , American Institute of Aeronautics and Astronautics, 2013/04/08 2006.doi:10.2514/6.2006-6753.[15] Montenbruck, O. and Gill, E.,
Satellite Orbits: Models, Methods, and Applications , Physics and astronomy online library,Springer, 2000.[16] Geuzaine, C. and Remacle, J.-F., “Gmsh: a Three-Dimensional Finite Element Mesh Generator with Built-in Pre-and Post-processing Facilities,”
International Journal for Numerical Methods in Engineering , Vol. 709, No. 11, 2009,pp. 1309–1331.[17] Webb, C. E.,
Radiation Force Modeling for ICESat Precision Orbit Determination , Ph.D. thesis, University of Texas atAustin, 2007.[18] Picone, J., Hedin, A., Drob, D., and Aikin, A., “NRL-MSISE-00 Empirical Model of the Atmosphere: Statistical Com-parisons and Scientific Issues,”