<< B-plane Cookbook Coordinates and frames >>

celestlab >> - Introduction - > Cookbook > Calculation of atmospheric drag

Calculation of atmospheric drag

Effect of atmospheric drag for a complex object

Introduction

This tutorial explains how to compute the drag area in a general case, and the effect of atmospheric drag on the orbital elements (semi major axis).

In particular, it shows the difference between considering a constant (average) value for the area exposed to the atmosphere and the actual area at each time.

1 - Create an object and an orbit

In order to compute the drag area, a model of an object (spacecraft) is needed if the computation cannot be done easily. The model may be created using STELA tool (an XML description can then be loaded in CelestLab), or it may be created directly. A simple satellite is created here: the body is a cube centered on the frame origin (size = l_body). It has 2 solar arrays of size l_array x w_array. The rotation axis is +X (X axis of body frame). The side with the solar cells initially looks towards +Z (Z axis of body frame) if the rotation angle is 0.

An orbit is also defined, and is propagated over one orbit period. The spacecraft frame is the defined such that the axes are parallel to the "qsw" local orbital frame, so that the solar arrays are correctly oriented even if they don't rotate.

// -----------------------------------------------------
// Load utility functions 
// (make_simple_spacecraft and plot_object)
// Look in specified directory for details
// -----------------------------------------------------
D = CL_path("SHAPES", fullfile(CL_home(), "data")); 
exec(fullfile(D, "shapes_util.sce"));

// -----------------------------
// Define a simple spacecraft 
// -----------------------------
// lengths in meters 
// margin_array: space between the body and the array
dim_body = [1; 1; 1]; // x, y, z (m)
dim_array = [1.5; 1]; // x, y (m)
gap = 0.2; // (m)
rot_angle = CL_deg2rad(0); // around x (rad)

obj = make_spacecraft(dim_body, dim_array, gap, rot_angle); 

// plot the spacecraft
scf(); 
plot_spacecraft(obj); 
xtitle("Spacecraft in Spacecraft''s frame"); 

// -----------------------------
// Create an orbit
// -----------------------------

// Define the orbit and propagate it 
// inclination = 90 deg, local time of asc. node = 12 h
cjd0 = CL_dat_cal2cjd(2020, 3, 21, 0, 0, 0); 

kep0 = [7.e6; 1.e-3; 98 * %pi / 180; 0; 0; 0]; 
kep0(5) = CL_op_locTime(cjd0, "mlh", 12, "ra"); 
T = CL_kp_params("per", kep0(1)); // orbit period 
t = linspace(0, T, 101); 
cjd = cjd0 + t / 86400; 

kep = CL_ex_kepler(cjd0, kep0, cjd); 
[pos, vel] = CL_oe_kep2car(kep); 

// Define the spacecraft frame (inertial frame -> spacecraft frame)
// Xsat: -angular momentum
// Zsat: radial direction (zenith)
Msat = CL_fr_locOrbMat("Wsq", pos, vel); 

// Velocity (relative to rotating Earth) in spacecraft frame 
// vel_ecf_eci: velocity / rotating Earth, coordinates in ECI   
[pos_ecf, vel_ecf] = CL_fr_convert("ECI", "ECF", cjd, pos, vel); 

// velocity of spacecraft / atmosphere (coord. in inertial frame) 
vel_atm = CL_fr_convert("ECF", "ECI", cjd, vel_ecf);

// plot the spacecraft in a specific frame ("AN" frame): 
// X -> ascending node, Z -> north pole
// Note : 2nd argument of plot_object: transformation matrix 
// from this frame to the spacecraft's frame
scf(); 
plot_spacecraft(obj, Msat(:,:,1) * CL_rot_angles2matrix(3, kep0(5))'); 
xtitle("Spacecraft in ''AN'' frame at initial epoch");

2 - Area over one orbit (fixed solar array)

The average drag area over one orbit is computed. The solar array is fixed (in the body frame). The orientation of the spacecraft is such that: Z axis points towards the zenith, and +Y points in the same direction as the velocity (relative to the atmosphere).

// -----------------------------------------------------
// Area over one orbit - fixed solar array.  
// -------------------------------------------------------
// First, run section 1 
// Then ...  
// -------------------------------------------------------

// Compute area 
opts = struct("resol", 200000); 
area = CL_stela_area(obj, "fix", Msat * vel_atm, opts = opts);

// Compute average area 
area_mean = CL_mean(area, meth = "trap"); 

scf(); 
plot(t, area); 
plot(t, area_mean * ones(t), "thickness", 2); 
xtitle("Area and mean area (m^2)"); 
CL_g_stdaxes(); 

// Note : 
// If the velocity / atmosphère was the inertial velocity, 
// the expected area would be : 
// l_body*l_body = 9 m^2

3 - Area over one orbit - rotating solar array (Sun pointing).

// -----------------------------------------------------
// Area over one orbit - rotating solar array.  
// -------------------------------------------------------
// First, run sections 1  
// Then ...  
// -------------------------------------------------------

// Compute the rotation angle
// =>> angle between Z axxis of solar array and Sun direction (in spacecraft 
// frame) as small as possible. 
dirsun = CL_eph_sun(cjd); 
dirsun_sat = Msat * dirsun; 
k = find(obj.id == "array1");
cell_axis = obj.z_axis(:,k); 
rot_axis = obj.x_axis(:,k)
[qa, rot_angle, Inok] = CL_rot_defRot1Ax(cell_axis, rot_axis, dirsun_sat, 0); 

// compute the rotation structure for the computation of the mean area
rotdesc_list = list(); 
rotdesc = struct(); 
rotdesc.id = "array1"; 
rotdesc.axis = rot_axis; 
rotdesc.ang = rot_angle; 
rotdesc_list($+1) = rotdesc; 

rotdesc = struct(); 
rotdesc.id = "array2"; 
rotdesc.axis = rot_axis; 
rotdesc.ang = rot_angle; 
rotdesc_list($+1) = rotdesc; 

opts = struct("resol", 200000); 
area = CL_stela_area(obj, "fixr", Msat * vel_atm, rotdesc_list, opts = opts); 

// compute average area 
area_mean = CL_mean(area, meth = "trap"); 

scf(); 
plot(t, area); 
plot(t, area_mean * ones(t), "thickness", 2); 
xtitle("Area and mean area (m^2)"); 
CL_g_stdaxes(); 

// Note : 
// If the velocity / atmosphere was the inertial velocity, the expected area 
// would be: l_body*l_body + 2 * (l_array * w_array) * (2/%pi) = 28.1 m^2

4 - Area x density over one orbit

// -----------------------------------------------------
// Area x density over one orbit  
// -----------------------------------------------------
// First, run section 1 and (2 or 3) 
// Then ...  
// -------------------------------------------------------

// compute atmosphere density (average solar activity) 
Flux = 150; 
Ap = 15; 

pos_ecf_ell = CL_co_car2ell(pos_ecf); 

rho = CL_mod_atmMSIS00(cjd, pos_ecf_ell(1,:), pos_ecf_ell(2,:), ...
                   pos_ecf_ell(3,:), Flux, Flux, Ap, res="dens"); 

// compares mean(density x area) with mean(density) x mean(area)
rho_mean = CL_mean(rho, meth = "trap"); 
rho_area_mean = CL_mean(area .* rho, meth = "trap"); 

UNIT = 1.e-13; 
scf(); 
plot(t, area .* rho / UNIT, "r"); 
plot(t, rho_area_mean * ones(t) / UNIT, "r", "thickness", 2); 
plot(t, area_mean * rho / UNIT, "b"); 
plot(t, area_mean * rho_mean * ones(t) / UNIT, "b", "thickness", 2); 
xtitle("Area x density (1.e-13 kg/m)"); 
CL_g_legend(gca(), ["var", "var (mean)", "const", "const (mean)"]); 
CL_g_stdaxes();

5 - Semi-major axis decrease rate over one orbit

// -----------------------------------------------------
// Effect of drag on semi major axis  
// -------------------------------------------------------
// First, run section 1, (2 or 3) and 4
// Then ...  
// -------------------------------------------------------

mass = 300;
Cx = 2.7;  
M_gauss = CL_op_gauss("kep", kep, "inertial"); 

// acceleration due to drag (variable area)
acc = CL_dMult(-vel_atm, CL_norm(vel_atm) .* rho .* area) * Cx / mass; 

dkepdt = M_gauss * acc; 
smadot = dkepdt(1,:); 
smadot_mean = CL_mean(dkepdt(1,:), meth = "trap");

// acceleration due to drag (mean area)
acc2 = CL_dMult(-vel_atm, CL_norm(vel_atm) .* rho .* area_mean) * Cx / mass; 

dkepdt2 = M_gauss * acc2; 
smadot2 = dkepdt2(1,:); 
smadot2_mean = CL_mean(smadot2, meth = "trap");

UNIT = 1 / 86400; 
scf(); 
plot(t, smadot / UNIT, "r"); 
plot(t, smadot_mean * ones(t) / UNIT, "r", "thickness", 2); 
plot(t, smadot2 / UNIT, "b"); 
plot(t, smadot2_mean * ones(t) / UNIT, "b", "thickness", 2); 
xtitle("Decrease rate of semi major axis (m/day)"); 
CL_g_legend(gca(), ["var", "var (mean)", "const", "const (mean)"]); 
CL_g_stdaxes();

6 - Area x Cx over one orbit *** TO DO ***

// -----------------------------------------------------
// mean value of Area x Cx  
// -----------------------------------------------------


Report an issue
<< B-plane Cookbook Coordinates and frames >>