reentry.py

#!/usr/bin/env python3

import numpy as np
from matplotlib import pyplot as plt

# Parachute data:
# http://www.laboratoridenvol.com/space/gnusoyuz/gnusoyuz.en.html
#
# Timings
# http://russianspaceweb.com/soyuz-ms-10.html#scenario

SOYUZ_MASS = 3000
FAIRING_MASS = 2500

SOYUZ_CD = 1.5
FAIRING_CD = 0.5
CHUTE_CD = 1.75
DROGUE_CD = 1.75

SOYUZ_RADIUS = 2.72
FAIRING_RADIUS = 3.0

DROGUE_AREA = 25
CHUTE_AREA = 1000
SOYUZ_AREA = np.pi * SOYUZ_RADIUS**2.0
FAIRING_AREA = np.pi * FAIRING_RADIUS**2.0

EARTH_RADIUS = 6371e3
EARTH_GRAV_K = 3.986004418e14

TEMPERATURE = 273.15 + 15
START_ALT = 50e3
SIM_DT = 0.5

# 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 dragcoef(t):
  return SOYUZ_CD if t > 80 else 0.25

def mass(t):
  return SOYUZ_MASS if t > 80 else SOYUZ_MASS + 3500

#returns Cd * A for all phases of flight
def drag_area(t):
  if t > 600:
    return SOYUZ_CD * SOYUZ_AREA + np.min([(t-600.0) / 10.0, 1.0]) * CHUTE_CD * CHUTE_AREA
  if t > 550:
    return SOYUZ_CD * SOYUZ_AREA + np.min([(t-550.0) / 3.0, 1.0]) * DROGUE_CD * DROGUE_AREA
  if t > 80:
    return SOYUZ_CD * SOYUZ_AREA
  # add parachute case
  return FAIRING_CD * FAIRING_AREA

#  - Aerodynamic drag
def drag(pos, vel, drag_a):
  rho = density(altitude(pos))
  v = np.linalg.norm(vel)
  return -(0.5 * rho * v * drag_a) * vel


#===---------------------------------------------------------------------------------------------===
# 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 =  (drag(p, v, drag_area(t[k])) + 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) <= 0:
      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)
  )

# Get the maximum acceleration from a run
@np.vectorize
def max_g(fpa, vel, dt=0.1, max_it=100000):
  x, y, vx, vy, ax, ay, t = sim_run(START_ALT, np.radians(fpa), vel, dt=dt, max_it=max_it)
  return np.max(np.sqrt(ax**2.0 + ay**2.0)) / 9.81

@np.vectorize
def max_t(fpa, vel, dt=0.1, max_it=100000):
  x, y, vx, vy, ax, ay, t = sim_run(START_ALT, np.radians(fpa), vel, dt=dt, max_it=max_it)
  return t[-1]

@np.vectorize
def plot_accel_t(fpa, vel, dt=0.1, max_it=100000):
  x, y, vx, vy, ax, ay, t = sim_run(START_ALT, np.radians(fpa), vel, dt=dt, max_it=max_it)
  a = np.sqrt(ax**2 + ay**2)
  plt.plot(t+80, a/9.81, label='%.0f°, %.0fm/s' % (fpa, vel))

@np.vectorize
def plot_traj(fpa, vel, dt=0.1, max_it=100000):
  x, y, vx, vy, ax, ay, t = sim_run(START_ALT, np.radians(fpa), vel, dt=dt, max_it=max_it)
  lat, alt = polar(x, y)
  downrange = (lat) * EARTH_RADIUS
  plt.plot(downrange/1e3, alt/1e3, label='%.0f°, %.0fm/s' % (fpa, vel))

@np.vectorize
def plot_vy(fpa, vel, dt=0.1, max_it=100000):
  x, y, vx, vy, ax, ay, t = sim_run(START_ALT, np.radians(fpa), vel, dt=dt, max_it=max_it)
  lat, alt = polar(x, y)
  #downrange = (lat) * EARTH_RADIUS
  print('vy=%f' % vy[-1])
  plt.plot(t, vy, label='%.0f°, %.0fm/s' % (fpa, vel))

#Do a parameter space analysis -- vary flight path angle and velocity

angles1 = np.linspace(40.0, 70.0, 20)
vels1 = np.linspace(1600.0, 2000.0, 20)

V, A = np.meshgrid(vels1, angles1)
G = max_g(A, V, dt=SIM_DT)

plt.figure()
CS = plt.contour(V, A, G)
plt.clabel(CS, inline=True)
plt.title('Maximum Soyuz Ballistic Acceleration')
plt.ylabel(r'$\alpha$ at sep.')
plt.xlabel(r'|v| at sep')
plt.tick_params(direction='in', which='both')
plt.savefig('soyuz-accel.pdf')
plt.savefig('soyuz-accel.png')
plt.show()

angles2 = np.linspace(50.0, 60.0, 2)
vels2 = np.linspace(1700, 1900, 3)
A2,V2 = np.meshgrid(angles2, vels2)

plt.figure()
plot_traj(A2, V2, dt=SIM_DT)
plt.title('Soyuz Abort Trajectory')
plt.xlabel('downrange (km)')
plt.ylabel('altitude (km)')
plt.legend()
plt.savefig('soyuz-traj.pdf')
plt.savefig('soyuz-traj.png')

plt.figure()
plot_accel_t(A2, V2, dt=SIM_DT)
plt.title('Soyuz Abort Acceleration Profile')
plt.xlabel('elapsed time since abort (s)')
plt.ylabel('acceleration (g)')
plt.legend()
plt.savefig('soyuz-at.pdf')
plt.savefig('soyuz-at.png')
plt.show()