Basic 2D simulation of lifting body spacecraft reentry

crew_dragon_reentry.py

#!/usr/bin/env python3

import numpy as np
from matplotlib import pyplot as plt
import matplotlib

### 130m/s deorbit burn from a 410km*410km orbit ###
#
# periapsis:      -31.87 km
# speed at EI:      7870.71 m/s
# Flight Path Angle at EI:  1.84 degrees

# Dragon data sources:
# https://smartech.gatech.edu/bitstream/handle/1853/26437/107-277-1-PB.pdf

# matplotlib.rcParams['font.size'] = 14
matplotlib.rcParams['mathtext.fontset'] = 'stix'
matplotlib.rcParams['font.family'] = 'STIXGeneral'
matplotlib.rcParams['lines.linewidth'] = 1.2

OUT_DIR       = './plots'
TRAJ_FPA      = -1.8
TRAJ_VEL      = 7870

CRAFT_CD      = 1.57
CRAFT_LD      = 0.18

CRAFT_RADIUS  = 3.66/2.0
CRAFT_MASS    = 8500
CRAFT_NAME    = "Crew Dragon"

CHUTE_CD = 1.75
DROGUE_CD = 1.75

DROGUE_AREA = 25
CHUTE_AREA = 1000
CRAFT_AREA = np.pi * CRAFT_RADIUS**2.0

EARTH_RADIUS = 6371e3
EARTH_GRAV_K = 3.986004418e14

TEMPERATURE = 273.15 + 15
START_ALT = 121e3
SIM_DT = 0.1

# Functions we'll use to get the atmospheric density during simulations. We could use the
# NASA Standard Atmosphere model, but scale height is probably a good enough approximation
# for this and it's much quicker to implement

SCALE_HEIGHT = 1e3 * TEMPERATURE * 1.38/(4.808*9.81)
    
def density(alt):
  return 1.221 * np.exp(- alt / SCALE_HEIGHT)

# Since we work in 2D, we need to be able to convert from x,y coordinates and lat/altitude
def cartesian(lat, alt):
  r = alt + EARTH_RADIUS
  return (r * np.cos(lat), r * np.sin(lat))

def polar(x, y):
  alt = np.sqrt(x**2 + y**2) - EARTH_RADIUS
  return (np.arctan2(x, y), alt)

def altitude(pos):
  return np.linalg.norm(pos) - EARTH_RADIUS

# Forces that will act on the body as it travels
#  - Gravity. Honestly, we could probably get similar results with 9.8
def gravity(pos, mass):
  r = np.linalg.norm(pos)
  return -((EARTH_GRAV_K * mass) / (r**2.0)) * (pos/r)

def mass(t):
  return CRAFT_MASS

def drag_area(alt):
  # TODO: we could try and change the returned value here to simulate chutes and all.
  #       Honestly, this is probably not worth it given how basic these sims are.
  return CRAFT_AREA * CRAFT_CD

#  - Aerodynamic drag + lift
def aero(pos, vel, drag_a):
  rho = density(altitude(pos))
  v = np.linalg.norm(vel)
  drag_mag = 0.5 * rho * v * drag_a
  lift_mag = drag_mag * CRAFT_LD
  
  drag = -drag_mag * vel
  lift = lift_mag * np.array([vel[1], vel[0]])
  #return -(0.5 * rho * v * drag_a) * vel
  return drag + lift


#===---------------------------------------------------------------------------------------------===
# And now we can actually simulate. 

def sim_run(alt, fpa, velocity, dt=0.1, max_it=100000):
  # Initialise our position, velocity, and acceleration components
  p = np.array([0, alt + EARTH_RADIUS])
  x = np.zeros(max_it)
  y = np.zeros(max_it)
  
  v = velocity * np.array([np.cos(fpa), np.sin(fpa)])
  vx = np.zeros(max_it)
  vy = np.zeros(max_it)
  
  a = np.array([0, 0])
  ax = np.zeros(max_it)
  ay = np.zeros(max_it)
  
  t = np.arange(0, max_it * dt, dt)
  
  # Very basic Euler integrator, running until impact with the ground
  # (doesn't take parachute into acount -- decent would slow down then)
  k = 0
  for _ in range(0, max_it):
    p = p + v * dt
    x[k], y[k] = p
    
    a =  (aero(p, v, drag_area(altitude(p))) + gravity(p, mass(t[k]))) / mass(t[k])
    ax[k], ay[k] = a
    
    v = v + a * dt
    vx[k], vy[k] = v
    
    k += 1
    if altitude(p) <= 10e3:
      print('done in %d iterations' % k)
      break
      
  return (
    np.resize(x, k),
    np.resize(y, k),
    np.resize(vx, k),
    np.resize(vy, k),
    np.resize(ax, k),
    np.resize(ay, k),
    np.resize(t, k)
  )


# Run the simulation
x, y, vx, vy, ax, ay, t = sim_run(
  START_ALT,
  np.radians(TRAJ_FPA),
  TRAJ_VEL,
  dt=SIM_DT,
  max_it=10000
)

title = r'Crew Dragon, $C_D=%.2f, L/D=%.2f$' % (CRAFT_CD, CRAFT_LD)
label = r'$fpa=%.2f, v_\mathrm{EI}=%.0f\mathrm{ms}^{-1}$' % (TRAJ_FPA, TRAJ_VEL)

# convert to useful cooridnates
lat, alt = polar(x, y)
downrange = (lat) * EARTH_RADIUS
vn = np.linalg.norm([vx, vy])

# Compute the velocity magnitude
v = np.sqrt(vx**2.0 + vy**2.0)

# Get the axial load ((ax,ay) projected onto (vx,vy))
a = np.abs((ax*vx + ay*vy)/v)

# Time and distance to go
tti = np.max(t) - t
dtg = np.max(downrange) - downrange


f1 = plt.figure(figsize=(10, 3))
plt.plot(downrange/1e3, alt/1e3, label=label)
plt.title(title)
plt.xlabel('downrange (km)')
plt.ylabel('altitude (km)')
plt.legend()
plt.tick_params(which='both', direction='in', length=5)
f1.tight_layout()
plt.savefig('%s/traj.pdf' % OUT_DIR)
plt.savefig('%s/traj.png' % OUT_DIR, dpi=300)

f2 = plt.figure()
plt.plot(a/9.81, alt/1e3, label=label)
plt.title(title)
plt.xlabel('axial loads (g)')
plt.ylabel('altitude (km)')
plt.legend()
plt.tick_params(which='both', direction='in', length=5)
f2.tight_layout()
plt.savefig('%s/load.pdf' % OUT_DIR)
plt.savefig('%s/load.png' % OUT_DIR, dpi=300)

f3 = plt.figure()
plt.plot(dtg/1e3, tti/60.0, label=label)
plt.title(title)
plt.xlabel('distance to splashdown (km)')
plt.ylabel('time to splashdown (min)')
plt.legend()
plt.tick_params(which='both', direction='in', length=5)
f3.tight_layout()
plt.savefig('%s/dtg.pdf' % OUT_DIR)
plt.savefig('%s/dtg.png' % OUT_DIR, dpi=300)