Eth: 0x00cce8E2e56a543abc084920eee3f88eFD0921ea

Thursday, March 1, 2012

How to Calculate Drag in LEO using Matlab

I'm going to highlight a simple model that can be used to calculate drag in Low Earth Orbit (LEO). Objects in orbit around earth experience atmospheric drag at altitudes of up to 1200km off the surface of the Earth. Spacecraft that orbit in LEO must compensate for drag or else their orbit will decay and they will burn up in the atmosphere. Anything that travels below 100km altitude is now so low that it will eventually either burn up or hit the surface of the planet. Drag is a complex force and there are many variables that one should take into account. However, we have simplified it down to the bare bones so you can understand the essentials in a simplified model. 

Things you will need: 
Data on your spacecraft you want to calculate drag on including: orbit altitude (assume circular orbit), surface area (assume spherical), mass of spacecraft, and the velocity of the spacecraft. Optional: A copy of Matlab to run example code.

Example Problem: 
For our example we will assume an altitude of 455 (km). From the velocity equation of a circular orbit Vcirular=sqrt(mu/radius). Where mu for the Earth is 398,600 (km^3/s^2). Radius = (altitude + radius of earth) = 455 (km) + 6378 (km) = 6833 (km).  Our velocity is therefore 7.64 (km/s) = 7640 (m/s). We will assume a spacecraft mass of 150 (kg). Furthermore we will assume our spacecraft is a sphere with a diameter of 5 (meters). 

Drag Equation:
Where p is the density of the atmosphere at that altitude, Cd is the coefficient of drag, V is our velocity and Area is the frontal area exposed to drag which is 1/2 the area of a sphere for our example. A good place to start is assuming a coefficient of drag of 2 to 2.2. This is typical of hypervelocity trajectories at the altitudes we are dealing with. Our velocity is 7640 (m/s) and our area for the front half of a sphere with a diameter of 5m is 39.26 (m^2). That leaves the density which is the most complicated part of the whole equation.

Estimate Density:
We will estimate the density of the environment using a model developed by James Wertz. This model assumes standard solar conditions and is broken up into altitude blocks. The reason it is broken up into blocks is because this creates a better estimate. The equation we will be applying is:
This equation says that the density a a point is equal to the nominal density (p0) multiplied by e to the power of the altitude over the scale height (H). Use the table below taken from Wertz' work to find the nominal density and scale height for the altitude of 455 (km).
Source: David A. Vallado, Wayne D. McClain - Fundamentals of astrodynamics and applications

For a 450 (km) altitude we will find our scale height to be 60.828 km and our nominal density to be 1.585*10^(-12) kg/(m^3). Plugging these values into the density equation we find a p value of 8.9425*10^(-16) kg/(m^3). 

Finalizing Drag:
We've found our density for the altitude we are orbiting (8.9425*10^(-16) kg/(m^3)), now we need to plug that along with everything else into our drag equation. The the area we found before is 39.26 (m^2). The drag coefficient (Cd) we are assuming is 2. Finally we found our velocity 7640 (m/s). Applying all these elements to the drag equation we find Fdrag =  -2.05*10^(-6) (kg*m/s^2). This force we can translate into an acceleration using the F=m*a relationship, where m is the mass of the spacecraft 150 (kg) and a is the acceleration (change in velocity) the spacecraft experiences from drag. 
Matlab Example File:
I've included a Matlab example file which you can call by putting it into your matlab directory and running: 
Note: This program is case sensitive. Altitude should be in km and Diameter should be in meters. This program calculates the velocity just like we did assuming a circular orbit and calculates the frontal surface area.
Happy modeling!

1 comment:

Mike Luby said...

Thanks for this! I just wrote it out in python.

I'm interested in how many orbits a thing will manage, which is what the orbital_decay function does a very crude job of estimating.

Do you have any suggestions for how to improve that estimation?