# Particulate Contamination Transport on Orbiting Satellites

Lately I’ve been working on simulating transport of particulate contaminants on an orbiting satellite. I am not going to get into the results but wanted to at least summarize the approach as I thought this was an interesting problem. And since my background in orbital mechanics is quite limited, this way I also hope to bring in a wider audience. **Don’t hesitate to chime in if the method is totally bogus!**

So what’s the issue? Spacecraft thermal surfaces are often painted with black or white paint, depending on the desired thermal properties. Paints can “particulate”, or shed small, dust-like particles. This process may arise from acoustic loads during launch, or from deployment of appendages on orbit. The question is: if particulates are released, what fraction of them comes back to contaminate the satellite?

Now, we know that in absence of forces, particles will move in a straight line and only impact surfaces within their line of sight. But, once forces are considered, the answer is less clear. But what forces act on particulates in space? First, **aerodynamic drag** can be neglected due to low pressure, especially for geostationary satellites. There are also **electrostatic** forces as particles can charge up. But if we assume that charge over mass remains small, these forces should also be negligible. There is however one more force that needs to be considered: **gravity**. You may be thinking, these dust flakes are very tiny, and the force of gravity will also be small. That may be the case, but what we are specifically talking about is the force of gravity affecting **orbital bodies**.

Take a look at the sketch in Figure 1. The bottom face of the spacecraft bus is a source of particulates. Even though there are no surfaces within a direct line of sight, we may end up with particulate contamination on the spacecraft once orbital motion is considered. Why is this? Well, in Keplerian orbital mechanics, there is a one-to-one relationship between position along an orbit and the orbital velocity. This means that if a particle is launched at some point \(\vec{x}_0\) with velocity offset from the parent body \(\Delta \vec{v}_0\ne \vec{0}\), then this particle will necessarily fall onto a different orbit from the parent object. Furthermore, in the absence of drag or perturbations (more on this below), the two orbits will again intersect at this point. Now, the period of the new orbit will be different from the parent body. But, as long as the initial delta velocity is small enough, the temporal difference between the two crossings will be small, and the dust flake will impact the satellite. This is basically the gist of what I attempted to model.

## Gravitational Force

The force acting on all orbiting bodies is given by Newton’s Law of Universal Gravitation,

$$\vec{F}_G = -\frac{GMm}{r^2} \frac{\vec{r}}{r}$$

In this formulation, \(\vec{F}_G\) is the force acting on a body with mass \(m\) by the body with mass \(M\). Vector \(\vec{r}\) is \(\vec{r}_2-\vec{r}_1\) and \(r=|\vec{r}|\). The gravitational constant \(G=6.674\times 10^{-11}\) Nm^{2} kg^{-2}. Once the force is computed, the incremental velocity is computed from

$$\Delta\vec{v}_G=\frac{\vec{F}}{m}\Delta t$$

and the position can be integrated using the standard Leapfrog integrator (care must be taken to use small enough time steps).

## Python Implementation

I used Python to test out the algorith before implemeting it in CTSP. The code (newton.py) is listed below.

# -*- coding: utf-8 -*- """ Demo of particulate contaminant transport on an orbiting satellite, solar drag not considered Created on Mon Jan 11 10:13:01 2016 https://www.particleincell.com/2016/particulate-orbit/ @author: lbrieda """ import numpy as np import pylab as pl def mag(a): return np.sqrt(np.dot(a,a)) Ra=42164e3 #Ra is the initial apoapsis radius Rp=1*Ra #Rp is the initial periapsis radius G= 6.67408e-11 #gravitational constant Me=5.972e24 GM = G*Me #GM=3.98e14 #G*M for earth v0 = 3070 dt=60 #dt is the time-step for numerical solution num_it=20000 #init pos=np.zeros((num_it,3)) pos[0]=np.array([Rp,0,0]); vel=np.zeros((num_it,3)) vel[0]=np.array([0,v0,0]) vel_dust = np.zeros_like(vel); pos_dust = np.zeros_like(pos); pos_dust[0,:] = pos[0,:] vel_dust[0,:] = vel[0,:] dp = [0,0,0] dv = [-0.05,0,0] pos_dust[0] += dp; vel_dust[0] += dv; #vel_dust[0] += np.append(200*(1-2*np.random.random(2)),0) #rewind r = pos[0] r_mag = mag(r) a=-GM*r/(r_mag**3) vel[0] -= 0.5*a*dt r = pos_dust[0] r_mag = mag(r) a=-GM*r/(r_mag**3) vel_dust[0] -= 0.5*a*dt for it in range(1,num_it): #r-vector, earth is at (0,0,0) r = pos[it-1] r_mag = mag(r) a=-GM*r/(r_mag**3) vel[it] = vel[it-1]+a*dt pos[it] = pos[it-1]+vel[it]*dt #dust r = pos_dust[it-1] r_mag = mag(r) a=-GM*r/(r_mag**3) vel_dust[it] = vel_dust[it-1]+a*dt pos_dust[it] = pos_dust[it-1]+vel_dust[it]*dt #load CTSP results #ctsp_sat = np.loadtxt("cg_pos.txt") #ctsp_part = np.loadtxt("part_pos.txt") fig = pl.figure(figsize=(18,8)) ax = fig.add_subplot(121) ax.plot(pos[:,0],pos[:,1],'k')#marker='o',markersize=4); # plot initial orbit ax.hold(True); ax.plot(pos_dust[:,0],pos_dust[:,1],'b'); # plot initial orbit #add CTSP orbit #ax.plot(ctsp_sat[:,0],ctsp_sat[:,1],'k:',linewidth="2",marker='o',markersize=4,markevery=100); #ax.plot(ctsp_part[:,3],ctsp_part[:,4],'b:',linewidth="2",marker='o',markersize=4,markevery=100); pl.legend(['Python sat','Python dust','CTSP sat','CTSP dust']) pl.title('Orbits for dv=(%.2f,%.2f,%.2f)'%(dv[0],dv[1],dv[2])) pl.axis("equal") ax = fig.add_subplot(122) d=pos_dust-pos ax.plot(d[:,0],d[:,1],'r-'); # plot initial orbit ax.hold(True); #ax.plot(ctsp_part[:,0],ctsp_part[:,1],'r:',linewidth="2",marker='o',markersize=4,markevery=10); pl.legend(['Python','CTSP']) pl.title('Relative position of dust') pl.savefig("plot.png",dpi=180,format="PNG",pad_inches=0.0,bbox_inches="tight")

This code includes several commented out lines that were used to verify the implementation in CTSP (when it comes to V&V of simulation codes, verification checks that equations are implemented correctly, while validation checks that the correct equations are implemented). Getting the code working in CTSP required performing transformations between local and earth inertial coordinate frames hence the need for this verification. Results with large injection delta velocity are shown below in Figure 2. Each of the four sets of images shows the orbit of the satellite (black) and dust (blue) on left. In this figure, Earth is at (0,0). The second set shows the relative position between the satellite and the particulate. Here (0,0) indicates collision. In both sets of plots, solid lines are results from the Python simulation and dotted lines with markers correspond to data loaded from CTSP. You won’t have this second set of data when running the demo code. The two lie on top of each which is what I wanted to see.

As you can see, the orbit of the dust relative to the satellite can take on some quite funky shapes! However, we don’t really see collisions, but this is due to the large injection velocity. Now, I don’t know at what velocity the particulates leave, but it’s probably closer to centimeters per second, instead of the 100 m/s considered above. When we use 5 cm/s delta velocity, we get the plot as shown below in Figure 3. Now the relative position routinely crosses through the (0,0) point, indicating the particle impacting the satellite.

## Solar Pressure

In reality, the relative orbit will not have the shape as shown above due to various perturbations, such as non-uniformities in the gravitational field. I neglected this one. Another effect is drag due to solar pressure. The scalar magnitude at 1 AU is obtained from the following equation

$$P_S=a \frac{W}{c} \cos^2\alpha$$

Here, \(W\) is the solar constant valued at 1361 W/m^{2} and \(c\) is the speed of light. The factor \(a\) has value of 1 for perfectly absorbing surfaces and 2 for reflecting surfaces. A colleague suggested to use 1.4 for diffusely reflecting surface so that’s what I used. The angle \(\alpha\) is the angle between the surface normal and the light incidence angle. In the code, I assumed that the flake, approximated as a flat plate, is randomly tumbling. I hence used an average value over \([-\pi/2:\pi/2]\) range (don’t quote me on this as not sure if this approximation is correct). This can be computed from

$$\overline{\cos^2\alpha}=\frac{1}{\frac{\pi}{2}-\frac{-\pi}{2}}\int_{-\pi/2}^{\pi/2}\cos^2\alpha \text{d}\alpha$$

The “cosine squared” integral is computed using the half angle identity,

$$\begin{align}

\overline{\cos^2\alpha} &= \frac{1}{\pi}\int_{-\pi/2}^{\pi/2}\left(\sqrt{\frac{1+\cos 2\alpha}{2}}\right)^2\alpha \text{d}\alpha\\

&=\frac{1}{2\pi}\int_{-\pi/2}^{\pi/2}1+\cos 2\alpha \text{d}\alpha\\

&=\frac{1}{2\pi}\left[\alpha+\frac{1}{2}\sin 2\alpha \right]_{-\pi/2}^{\pi/2}\\

&=\frac{1}{2}

\end{align}

$$

What’s interesting about solar pressure is that for half of the orbit, it pushes on the particle, while for the half, it pushes against it. This also means we need to know the position of the sun relative to the particle. I used a function from Cowel’s Method for Earth Satellites code by David Eagle to compute the position of the sun in the ECI coordinates. We can then compute the solar radiation pressure force acting on the particle as

$$\vec{F}_S = P_SA \frac{\vec{x}_{p}-\vec{x}_{S}}{|\vec{x}_{p}-\vec{x}_{S}|}$$

where the last term is simply the unit vector denoting the direction from the sun to the particle. Note that in case of the gravitational force, we didn’t actually need the mass of the particle, since that term dropped out when computing \(d\vec{v}/dt\equiv\vec{F}/m\). But we need it to compute the acceleration due to solar pressure. In the code below, I assumed that the dust flake is a brick 1mm x 2mm x 0.5mm with density of 2000 kg/m^{3}. These was just an assumption, so again, don’t quote me on this.

In order to keep this post shorter, I skipped including the text of the second code here, but can find it here: sun.py. When you run this code, you will get a plot like shown in Figure 4. The first figure is result after 2 days, while the second one is after a year. The black line is the unperturbed relative position, while the magenta line corresponds to the result with solar drag. As you can see, they are quite different!

And that’s it. Again, feel free to comment especially if you find some bugs in the code.

### Source Code

The source code for the two Python codes is here:

newton.py

sun.py

Why do you assume charge effects would remain small? Could sunlight cause particle charging by photoemission? Then there might be a small force between the neutral satellite and charged particle.

Hi Steve, I agree that this assumption is bit hairy and needs work. The real answer is that this effect was out of scope of the analysis. But regardless, I still would expect the q/m ratio to be quite small for the reason that I imagine only the “surface layer” atoms possibly getting an electron knocked off by photoemission, while all the internal atoms, making up bulk of the weight, remain neutral. But perhaps this is completely incorrect. I would be glad to hear your insight if this is something you specialize in.

It’s not my specialty, unfortunately, but a paper I read on satellite electrostatic modeling mentioned the UV photoemission effect kept the satellite positively charged.

Hmm, some googling shows that solar plasma negatively charges satellites up to the kilovolt range when the satellite is in shadow, but in sunlight, photoemission neutralizes positively charges it (to lower voltage, though). http://holbert.faculty.asu.edu/eee460/spacecharge.html

I bet you’re right that the net charge compared to the mass is small.. but even a tiny force effect will have an impact on the long term dynamics.

We’d have to learn more about such charging intensities and put it in your simulation to see if it had any effect.