Comet Dynamics II

So I finished the comet model for now (for a second time). I have reduced how long the model takes to run again from 33 seconds to 14 seconds by taking out the cross function and replacing it with the direct equation.

Changes from the last post:

• Variable inertia matrix

The inertia matrix now changes with the ejected mass. This makes the simulation more realistic. However, the shape of the comet does not change (ie. the change in mass does not change the comet thickness). 

• Random off-gassing jet placement and variable composition

Instead of off-gassing jets being able to be produced on every face, now the chance that a jet will be on a specific face is randomized (I currently have it at a 5% chance). The composition of those faces are then randomized at a maximum of 5% of the face being volatilizes. This also effects the initial inertia matrix.

One implication of this is with constant face jets, the comet usually will spin up along a specific rotation axis. This can be seen in the following video. The great thing is that the axis it settles on does not have the be a major or minor axis due to the off-gassing thrust.

Black line: Rotation axis.

Yellow line: Sun vector .

Comet Color: Proportional to temperature (Red is hot).

Dark blue line: Total force on body.

Rings: Torque direction (green is +, red is -) and normalized magnitude. 

Light blue lines: Direction and normalized magnitudes of off-gassing jets.

Off-Gassing Dynamics

The off-gassing of an asteroid also has implications in satellite dynamics. However, different size and velocity particles are created and can effect a spacecraft in different ways, ie. it can totally destroy it, puncture a solar panel, provide a torque, etc.. All of which would make precise data acquisition difficult. The second part of my thesis project will be analyzing this.

• Off-Gassing

On top of the trust off-gassing provides to the asteroid, it also provides a pressure on spacecraft, as well as accelerates dust particulate off the surface of the comet. To be able to analyze this, I have made a simple model to find the velocity that specific size dust would travel at. The velocity of each particle is a function of the radius of the comet, the coefficient of drag, the solar flux density, a diffusion coefficient, radius of the comet, the gas ejecta density and velocity, and the particle density and radius. 

In most literature I have seen, out-gas jets are taken to have a half-angle and not diverge from that. However, since the pressure will not be constant through the coma, I created the diffusion coefficient. There is no empirical data behind this except that a high pressure jet would not stay at a constant half-angle due to pressure differentials between the jet and a low density area / vacuum, ie. it would defuse. Also, literature says the majority of mass will be at 20% of the velocity of the jets and usually be between 20 and 100 m/s. My off-gassing model adheres to that for dust larger than 1 mm. I believe this makes sense, but that's because I made it this way. 

The following plot is the results from this model:

Dust Radius: [0.0001;0.0005;0.001;0.0025;0.005;0.01;0.02;0.03;0.04;0.05] (m)

Gas Velocity: 400 (m/s)     Gas Density: 0.004 (kg/m^3)     Dust Density: 300 (kg/m^3)     

Diffusion Ratio: 1.5           Coefficient of drag: 0.5              Comet Radius: 500 (m)

Left: Velocity of 10 different particles

Middle: Distance from initial location of particles

Right: Median acceleration of the particles

Conclusions:

While I showed one run of the model here, it can be seen that the particles settle on a certain velocity relatively quick (withing 15 seconds for any "reasonable" parameter combination). I need to do further research or develop a formula to find the velocity from the above parameters without having to run a quick simulation. On top of that, at close quarters to the comet (within 1000 meters, the acceleration due to the gas dominates, whereas the acceleration due to gravity and solar radiation pressure gain prominence after a while. 

So whats next:

• Satellite Dynamics

Here is a quick model of the satellite I will be using (it's the same as past satellite pictures, but now it has an "engine" in the back and a "camera" in the front). Behold the JimSat in all it's glory:

JimSat

I will be using this to analyze a satellites ability to acquire data in highly dynamic offgassing/asteroid environment.

• Coma and Tail Model

I'm not too sure if I will be able to do this, but with the given trajectory and composition I want to predict the tail and brightness of the coma. Maybe even predict if the comet will give the Earth a meteor shower, I don't know.

Comet Dynamics

So I finished the comet model for now. It works well for me to move on and work on other parts of the model. I have verified it with literature and it gets similar results (which was a nice thing to see), but it still needs some work. It should be noted that the only "volatile" taken into account so far is ice/water vapor. 

Before I get into the bulk of this post, I found a few edits which make the model run magnitudes faster (~33 seconds opposed to 40 minutes). 

   - Due to the differential equation solver (ode113), if the axes lengths of the ellipse are equal or if the rotation axis is along a ellipse axis, the model will take forever (about 40 minutes, depending on simulation time) to run. This is because the differential equation solver wants to keep the angular velocity along the same axis, or having an angular acceleration of [0;0;0]. By making the rotation axis slightly offset from an ellipse axis, this is avoided since it neither wants to keep a constant rotation axis as well as constant angular velocity. I currently have the rotation axis at [pi/1000;pi/1000], or 0.18 degrees, and it works amazingly.

• Comet Model (Asteroid model with off-gassing)

The difference between the asteroid model and comet model is the thrust and energy use from off-gassing. The temperature of the face is used to predict the sublimation rate. Since the comet is in a vacuum, the velocity of the gas should be mach 1, which is found easily using values for water vapor for molar mass and heat capacity ratio. The force is also found and is shown to greatly effect the rotation of the comet.

So without further explanation, here are the same plots and videos from the last two posts, but with off-gassing. The first plot is a comet with axes a=5000, b=3000, and c=2000;, the second is scaled down by 100. 

Large Comet:

Small Comet: 

Conclusions: 

The time scale for both the first and second simulations was the same (just under 1 day's time). Notice the relatively constant temperature of the fast rotating comet, as shown in the video, this makes the force tend to be in the direction of travel, instead of along the anti-sun vector. This has some large implications for small comets because this means that they would (rather quickly) increase their semi-major axis until the off-gassing reduced to a negligible amount. Also, it can be concluded that small comets either will break apart due to their spin rate, become tumblers, or break apart from internal pressure build up. 

YORP Update

This will be quick post. I have found a few errors in the YORP model and corrected them. I was wondering why the force plot looked as it did, the reason is because it was only recording one face. As you can see below, this has been fixed.

Following are a few things to note and the new output plot as well as a video I made to validate some aspects of the model.

Things to note:

• The torques and forces are in the body frame. To make it useful it needs to be put into the orbital frame (check the video).
• The average temperature reaches steady state as the number of faces increases and the asteroid becomes circular about the rotation axis (if the rotational and orbital angular momentum axes also align).
• The warmest time on the asteroids is in the afternoon (around 2-4).

For the video below, the yellow line is the direction of solar rays, the red line is the direction of the force due to solar radiation pressure and thermal emittance. The color of the asteroid is relative to it's temperature with the colormap [1 0 0] (red) being the total hottest temperature and [0 0 0] being absolute zero. 

Asteroid Rotation Video viewed along the Z-axis

(Sorry about the resolution, uploading it requires compressing)

New output plot

Asteroid Dynamics

Finals are done and I am on winter break, so I finally have time to work on my thesis model.

Before getting into this, here are a few definitions:
    -Solar Radiation Pressure: Momentum transfer from photons which provide a force onto any body.
    -Yarkovsky Effect: This effect takes place due to the radiation pressure that a rotating body emits. Due to it's rotation, the force applied to the body can, over large periods of time, alter the orbital parameters of the asteroid.
    -YORP Effect: Second order Yarkovsky effect where, if the force is off axis, it provides a torque to the body and affects the rotational state of the asteroid.

Over the last few days, I have been putting a YORP model together, which I will then put into a Yarkovsky / bulk model, which I have yet to create (I had to make one first, I decided I would make the YORP). On top of that, I also need to add a Cometary function into the YORP model, which will be done after the bulk model. This model will include off-gassing. The total model will be broken up into bulk and body parameters. The bulk parameters being orbital based (radius and velocity) while the body parameters being rotation based (angle and angular velocity).

The way I currently see it going, the YORP model will be nested inside the Yarkovsky. I hope to find a way to simplify the model in an individual time dependent equation for force and torque (change in body and bulk parameters), but I don't think this will be possible for a more accurate asteroid model (not an ellipsoid). On that note, right now I use an ellipsoid for ease of computation. The lower the number of faces, the quicker it works. Once I have the model validated for the simplified model, I will run it on a more complex shape model.

Above is current output. Going counter-clockwise from the top left, the individual subplots are as follows:

-Picture of the asteroid and the rotation axis (red)
-Average temperature
-Temperature of every face
-Unit quaternion describing the angle
-Angular velocity in x-y-z directions
-Total force on surface (x,y,z)
-Torque on body

Asteroid Gravity Models

I did this a while back but decided only now to post it. Right now I have 3 gravity models (2 working). The problem with gravity models is that for concave structures, a tetrahedron from a face to the origin may have elements external to the body or overlap other tetrahedra. This "creates" mass/volume outside of the actual object volume which can cause huge errors for simulating an orbiting satellite.

A few ways around this are:
   1). to create two "origins" and set the origin for each face as the one where the inner product (dot product) of the CCW face normal and radius vector are positive (if both are positive take the one with the smaller radius).
   2). Define the mass to exclusively be on the surface.
   3). Do some math and define the gravity as a field affected by the body volume.

Originally, I was only using a point mass tetrahedra model. This model took the faces, connected them to the origin (creating a tetrahedron), found the volume, and put a point mass at the center of the tetrahedron proportional to the volume and density. This model does not work for very irregularly shaped objects, however it works best for the majority of asteroids.

The other working model is a shell model. It takes all of the mass at the surface of the asteroid where each unit of surface area contains the same amount of mass. While this model is better for more irregularly shaped objects, by putting all the mass at the surface it creates a gravity singularity at the surface and can mess up data if trying to analyze an orbit close to the asteroid. This model is best for analyzing the gravity field for very irregularly shaped asteroids at a reasonable distance from the surface.

The following plots are for validating these two models. I will discuss the differences in the plots and what denotes a good "fit" at the end.

It should be noted that the scale at the bottom is the number of points in the analysis direction and does not correlate to the size of the objects.

Eros
    

        

Left: Point Gravity Tetrahedra Model

Right: Surface Area Model

Churyumov-Gerasimenko (67P)

Left: Point Gravity Tetrahedra Model

Right: Surface Area Model

What to take from this:

Point Mass Tetrahedra

• Underestimates radius since it takes the center of the tetrahedra.
• This conserves mass while avoiding the gravity singularity at the surface
• For very irregularly shaped objects (67P), the volume is disconfigured and does not give accurate information about the shape/gravity of the object.

Surface Area

• Preserves shape of irregularly shaped objects (67P)

Conclusions

• The Point Mass Tetrahedra model should be used for regular shaped objects.
• The Surface Area model should be used for irregularly shaped objects. But may create large errors is simulating a close approach.
• Another model needs to be created to better represent this.
• A mixture of the two may work. I am also looking into Spherical Harmonics as well as a triangulation method I read about a while ago. But for now these will do.

Not SpaceSim

I will discuss 2 different programs in this post: a 2D finite element analysis program and a Robot arm. In my next post I will discuss my progress for my thesis.

I will not go too in-depth, so if there are any questions feel free to comment. 

Here is a quick picture of Dawn. I found the .stl file online and wrote a quick parser to make my program able to plot it.

2D FEA:

Without getting to deep into it, I made an FEA program for my Intro to Finite Element's class (my partner is doing the 3D analysis in Abaqus). The object being analyzed is the Philae Lander's leg. We are doing static analyses for the stress, strain, internal forces, and deformation for a resting Philae as well as when it was compacted before the first bounce. Assuming no energy dissipation (pretty bad assumption since there was lots of kinetic dissipation), the maximum internal energy is proportional to the square of the velocity after the bounce, which is estimated to be around 38 cm/s.

Here is what it looks like:

2D FEA Program

Robotics:

While we have not yet completed the project, my main contribution so far has been this robot which is able to get out of a randomly generated enclosure. It does this by always turning right, after the first left. It saves where it has been as an artificial "memory" (very basic AI), and if it "sees" the same spot over and over it randomizes it's movements until a new wall has been found.

I created the simulated robot since we were having trouble interfacting with the real robot (PhantomOmni).

Oh yeah, on that note I figured out how to take complex videos on matlab, the buffer needs to be set as "zbuffer" otherwise it does not do something that it has to do.

Anyways, here is a video. For best results try to play both videos at the same time:

Robot Arm

Robot Artificial Memory