Section 4-7 - Venus SmallSat Aerocapture - Performance Analysis
[1]:
from AMAT.planet import Planet
from AMAT.vehicle import Vehicle
[2]:
import numpy as np
from scipy import interpolate
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
from matplotlib.patches import Polygon
Small Run N = 10
[3]:
# Create a planet object
planet=Planet("VENUS")
planet.h_skip = 150000.0
planet.h_low=10.0E3
# Load an nominal atmospheric profile with height, temp, pressure, density data
planet.loadAtmosphereModel('../../../atmdata/Venus/venus-gram-avg.dat', 0 , 1 ,2, 3)
# Set up a vehicle object
vehicle=Vehicle('SmallSat1', 37, 20, 0.0, 1.767, 0.0, 0.35, planet)
vehicle.setInitialState(150.0,-6.22,23.16,10.8191,89.9982,-5.20,0.0,0.0)
vehicle.setSolverParams(1E-6)
vehicle.setDragModulationVehicleParams(20, 7.5)
# Set up the drag modulation entry phase guidance parameters.
vehicle.setDragEntryPhaseParams(2.0, 20.0, 101, -150.0)
# Set the target orbit parameters.
vehicle.setTargetOrbitParams(200.0, 2000.0, 20.0)
# Define the path to atmospheric files to be used for the Monte Carlo simulations.
atmfiles = ['../../../atmdata/Venus/LAT20N.txt']
# Set up the Monte Carlo simulation for drag modulation.
# NPOS = 151, NMONTE = 200
# Target EFPA = -5.30 deg
# EFPA 1-sigma error = +/- 0.033 deg
# Nominal beta_1 = 20 kg/m2
# beta_1 1-sigma = 0.0
# guidance time step for entry = 1.0s (Freq. = 1 Hz)
# guidance time step after jettison = 1.0 s
# max. solver time step = 0.1 s
# max. time used by solver = 2400 s
vehicle.setupMonteCarloSimulationD(151, 200, atmfiles, 0 , 1, 2, 3, 4, True,
-5.20, 0.0333, 20.0, 0.0,
1.0, 1.0, 0.1, 2400.0)
# Run 10 trajectories
vehicle.runMonteCarloD2(10, '../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/MCBX1')
RUN #: 1, PROF: ../../../atmdata/Venus/LAT20N.txt, SAMPLE #: 197, EFPA: -5.22, SIGMA: -0.07, APO : 1997.74
RUN #: 2, PROF: ../../../atmdata/Venus/LAT20N.txt, SAMPLE #: 2, EFPA: -5.26, SIGMA: -1.89, APO : 1998.12
RUN #: 3, PROF: ../../../atmdata/Venus/LAT20N.txt, SAMPLE #: 170, EFPA: -5.16, SIGMA: 1.32, APO : 2028.43
RUN #: 4, PROF: ../../../atmdata/Venus/LAT20N.txt, SAMPLE #: 154, EFPA: -5.20, SIGMA: 1.40, APO : 2090.60
RUN #: 5, PROF: ../../../atmdata/Venus/LAT20N.txt, SAMPLE #: 62, EFPA: -5.19, SIGMA: 1.01, APO : 2036.95
RUN #: 6, PROF: ../../../atmdata/Venus/LAT20N.txt, SAMPLE #: 71, EFPA: -5.29, SIGMA: 1.85, APO : 1196.56
RUN #: 7, PROF: ../../../atmdata/Venus/LAT20N.txt, SAMPLE #: 130, EFPA: -5.20, SIGMA: -1.04, APO : 2059.06
RUN #: 8, PROF: ../../../atmdata/Venus/LAT20N.txt, SAMPLE #: 174, EFPA: -5.20, SIGMA: -0.33, APO : 1873.30
RUN #: 9, PROF: ../../../atmdata/Venus/LAT20N.txt, SAMPLE #: 1, EFPA: -5.21, SIGMA: -0.45, APO : 1945.85
RUN #: 10, PROF: ../../../atmdata/Venus/LAT20N.txt, SAMPLE #: 132, EFPA: -5.12, SIGMA: -0.10, APO : 2029.68
Production Run N = 1000
Run the file section-4-7-venus-smallsat-performance-analysis-prd.py
[3]:
p1 = np.loadtxt('../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/MCB-Venus-1/terminal_periapsis_arr.txt')
x1 = np.loadtxt('../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/MCB-Venus-1/terminal_apoapsis_arr.txt')
y1 = np.loadtxt('../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/MCB-Venus-1/acc_net_g_max_arr.txt')
z1 = np.loadtxt('../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/MCB-Venus-1/q_stag_max_arr.txt')
h1 = np.loadtxt('../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/MCB-Venus-1/heatload_max_arr.txt')
u1 = np.loadtxt('../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/MCB-Venus-1/periapsis_raise_DV_arr.txt')
v1 = np.loadtxt('../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/MCB-Venus-1/apoapsis_raise_DV_arr.txt')
[10]:
fig, axes = plt.subplots(1, 2, figsize = (6.5, 3.25))
plt.rc('font',family='Times New Roman')
params = {'mathtext.default': 'regular' }
plt.rcParams.update(params)
fig.tight_layout()
plt.subplots_adjust(wspace=0.30)
axes[0].axvline(x=2000, linewidth=1.5, color='k', linestyle='dashed')
axes[0].hist(x1, bins=20, color='xkcd:periwinkle')
axes[0].set_xlabel('Apoapsis altitude, km', fontsize=12)
axes[0].set_ylabel('Number of cases', fontsize=12)
axes[0].tick_params(direction='in')
axes[0].yaxis.set_ticks_position('both')
axes[0].xaxis.set_ticks_position('both')
axes[0].tick_params(axis='x',labelsize=12)
axes[0].tick_params(axis='y',labelsize=12)
axes[1].axhline(y=2000, linewidth=1.5, color='k', linestyle='dashed')
axes[1].plot(p1, x1, 'bo', markersize=5, mfc='None')
axes[1].set_xlabel('Periapsis altitude, km',fontsize=12)
axes[1].set_ylabel('Apoapsis altitude, km', fontsize=12)
axes[1].tick_params(direction='in')
axes[1].yaxis.set_ticks_position('both')
axes[1].xaxis.set_ticks_position('both')
axes[1].tick_params(axis='x',labelsize=12)
axes[1].tick_params(axis='y',labelsize=12)
for ax in axes:
ax.xaxis.set_tick_params(direction='in', which='both')
ax.yaxis.set_tick_params(direction='in', which='both')
ax.xaxis.set_tick_params(width=1, length=8)
ax.yaxis.set_tick_params(width=1, length=8)
ax.xaxis.set_tick_params(width=1, length=6, which='minor')
ax.yaxis.set_tick_params(width=1, length=6, which='minor')
ax.xaxis.grid(which='major', color='k', linestyle='dotted', linewidth=0.5)
ax.xaxis.grid(which='minor', color='k', linestyle='dotted', linewidth=0.0)
ax.yaxis.grid(which='major', color='k', linestyle='dotted', linewidth=0.5)
ax.yaxis.grid(which='minor', color='k', linestyle='dotted', linewidth=0.0)
for axis in ['top', 'bottom', 'left', 'right']:
ax.spines[axis].set_linewidth(2)
plt.savefig('../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/venus-smallsat-apoapsis-vs-periapsis.png', dpi= 300,bbox_inches='tight')
plt.savefig('../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/venus-smallsat-apoapsis-vs-periapsis.pdf', dpi=300,bbox_inches='tight')
plt.savefig('../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/venus-smallsat-apoapsis-vs-periapsis.eps', dpi=300,bbox_inches='tight')
plt.show()
[7]:
fig, axes = plt.subplots(1, 2, figsize = (6.5, 3.25))
plt.rc('font',family='Times New Roman')
params = {'mathtext.default': 'regular' }
plt.rcParams.update(params)
fig.tight_layout()
plt.subplots_adjust(wspace=0.30)
axes[0].plot(y1, z1, 'ro', markersize=5, mfc='None')
axes[0].set_xlabel('Peak deceleration, Earth g', fontsize=12)
axes[0].set_ylabel('Peak stagnation-point heat rate, '+r'$W/cm^2$', fontsize=12)
axes[0].tick_params(direction='in')
axes[0].yaxis.set_ticks_position('both')
axes[0].xaxis.set_ticks_position('both')
axes[0].tick_params(axis='x',labelsize=12)
axes[0].tick_params(axis='y',labelsize=12)
axes[0].axhline(y=np.percentile(z1, 99.87), linewidth=1.5, color='k', linestyle='dotted')
axes[0].axvline(x=np.percentile(y1, 99.87), linewidth=1.5, color='k', linestyle='dotted')
axes[1].plot(z1, h1, 'go', markersize=5, mfc='None')
axes[1].set_xlabel('Peak stagnation-point heat rate, '+r'$W/cm^2$', fontsize=12)
axes[1].set_ylabel('Peak stagnation-point heat load, '+r'$kJ/cm^2$', fontsize=12)
axes[1].tick_params(direction='in')
axes[1].yaxis.set_ticks_position('both')
axes[1].xaxis.set_ticks_position('both')
axes[1].tick_params(axis='x',labelsize=12)
axes[1].tick_params(axis='y',labelsize=12)
axes[1].axhline(y=np.percentile(h1, 99.87), linewidth=1.5, color='k', linestyle='dotted')
axes[1].axvline(x=np.percentile(z1, 99.87), linewidth=1.5, color='k', linestyle='dotted')
for ax in axes:
ax.xaxis.set_tick_params(direction='in', which='both')
ax.yaxis.set_tick_params(direction='in', which='both')
ax.xaxis.set_tick_params(width=1, length=8)
ax.yaxis.set_tick_params(width=1, length=8)
ax.xaxis.set_tick_params(width=1, length=6, which='minor')
ax.yaxis.set_tick_params(width=1, length=6, which='minor')
ax.xaxis.grid(which='major', color='k', linestyle='dotted', linewidth=0.5)
ax.xaxis.grid(which='minor', color='k', linestyle='dotted', linewidth=0.0)
ax.yaxis.grid(which='major', color='k', linestyle='dotted', linewidth=0.5)
ax.yaxis.grid(which='minor', color='k', linestyle='dotted', linewidth=0.0)
for axis in ['top', 'bottom', 'left', 'right']:
ax.spines[axis].set_linewidth(2)
plt.savefig('../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/venus-smallsat-g-load-heat-rate.png', dpi= 300,bbox_inches='tight')
plt.savefig('../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/venus-smallsat-g-load-heat-rate.pdf', dpi=300,bbox_inches='tight')
plt.savefig('../../../data/mdpi-aerospace/smallsat-mission-concepts/venus/venus-smallsat-g-load-heat-rate.eps', dpi=300,bbox_inches='tight')
plt.show()
[8]:
fig, axes = plt.subplots(1, 2, figsize = (6.5, 3.25))
plt.rc('font',family='Times New Roman')
params = {'mathtext.default': 'regular' }
plt.rcParams.update(params)
fig.tight_layout()
plt.subplots_adjust(wspace=0.30)
axes[0].hist(u1, bins=10, color='xkcd:periwinkle')
axes[0].set_xlabel('Periapse raise '+r'$\Delta V$'+', m/s', fontsize=10)
axes[0].set_ylabel('Number of cases', fontsize=12)
axes[0].tick_params(direction='in')
axes[0].yaxis.set_ticks_position('both')
axes[0].xaxis.set_ticks_position('both')
axes[0].tick_params(axis='x',labelsize=12)
axes[0].tick_params(axis='y',labelsize=12)
axes[1].hist(np.abs(v1), bins=10, color='xkcd:orange')
axes[1].set_xlabel('Apoapsis correction '+r'$\Delta V$'+', m/s', fontsize=12)
axes[1].set_ylabel('Number of cases', fontsize=12)
axes[1].tick_params(direction='in')
axes[1].yaxis.set_ticks_position('both')
axes[1].xaxis.set_ticks_position('both')
axes[1].tick_params(axis='x',labelsize=12)
axes[1].tick_params(axis='y',labelsize=12)
for ax in axes:
ax.xaxis.set_tick_params(direction='in', which='both')
ax.yaxis.set_tick_params(direction='in', which='both')
ax.xaxis.set_tick_params(width=1, length=8)
ax.yaxis.set_tick_params(width=1, length=8)
ax.xaxis.set_tick_params(width=1, length=6, which='minor')
ax.yaxis.set_tick_params(width=1, length=6, which='minor')
ax.xaxis.grid(which='major', color='k', linestyle='dotted', linewidth=0.5)
ax.xaxis.grid(which='minor', color='k', linestyle='dotted', linewidth=0.0)
ax.yaxis.grid(which='major', color='k', linestyle='dotted', linewidth=0.5)
ax.yaxis.grid(which='minor', color='k', linestyle='dotted', linewidth=0.0)
for axis in ['top', 'bottom', 'left', 'right']:
ax.spines[axis].set_linewidth(2)
plt.show()
[4]:
def percent_between(a, lower, upper, closed_left=True, closed_right=True):
"""
Finds the percentage of values between a range for a numpy array
Parameters
----------
a: np.ndarray
numpy array to calculate percentage
lower: int, float
lower bound
upper: int, float
upper bound
closed_left:
closed left bound ( > vs >= )
closed_right:
closed right bound ( < vs <= )
"""
l = np.greater if not closed_left else np.greater_equal
r = np.less if not closed_right else np.less_equal
c = l(a, lower) & r(a, upper)
return np.count_nonzero(c) / float(a.size)
[12]:
percent_between(x1, 1600, 2400)
[12]:
0.914
[13]:
percent_between(x1, 1400, 2600)
[13]:
0.962
[14]:
percent_between(x1, 1200, 2800)
[14]:
0.984
[15]:
percent_between(x1, 1000, 3000)
[15]:
0.99
[5]:
print("--------------------------------------------------------------------------------------------------------")
print("Parameter | Min. 5-percentile Mean 95-percentile Max")
print("--------------------------------------------------------------------------------------------------------")
print("Apoapsis altitude, km | " + str(round(np.min(x1))) + " " + str(round(np.percentile(x1, 5))) +\
" " + str(round(np.average(x1))) +
" " + str(round(np.percentile(x1, 95)))+\
" " + str(round(np.max(x1))))
print("Periapsis altitude, km | " + str(round(np.min(p1), 1)) + " " + str(round(np.percentile(p1, 5), 1)) +\
" " + str(round(np.average(p1), 1)) +
" " + str(round(np.percentile(p1, 95), 1))+\
" " + str(round(np.max(p1), 1)))
print("Peak deceleration, g | " + str(round(np.min(y1), 2)) + " " + str(round(np.percentile(y1, 5), 2)) +\
" " + str(round(np.average(y1), 2)) +
" " + str(round(np.percentile(y1, 95), 2))+\
" " + str(round(np.max(y1), 2)))
print("Peak heat rate, W/cm2 | " + str(round(np.min(z1))) + " " + str(round(np.percentile(z1, 5))) +\
" " + str(round(np.average(z1))) +
" " + str(round(np.percentile(z1, 95)))+\
" " + str(round(np.max(z1))))
print("Peak heat load, kJ/cm2 | " + str(round(np.min(h1),1)) + " " + str(round(np.percentile(h1, 5), 1)) +\
" " + str(round(np.average(h1), 1)) +
" " + str(round(np.percentile(h1, 95), 1))+\
" " + str(round(np.max(h1), 1)))
print("Periapsis raise dV, m/s| " + str(round(np.min(u1), 1)) + " " + str(round(np.percentile(u1, 5), 1)) +\
" " + str(round(np.average(u1), 1)) +
" " + str(round(np.percentile(u1, 95), 1))+\
" " + str(round(np.max(u1), 1)))
print("Apoapsis corr. dV, m/s | " + str(round(np.min(v1), 2)) + " " + str(round(np.percentile(v1, 5), 2)) +\
" " + str(round(np.average(v1), 2)) +
" " + str(round(np.percentile(v1, 95), 2))+\
" " + str(round(np.max(v1), 2)))
--------------------------------------------------------------------------------------------------------
Parameter | Min. 5-percentile Mean 95-percentile Max
--------------------------------------------------------------------------------------------------------
Apoapsis altitude, km | 1399 1865 2107 2510 5418
Periapsis altitude, km | 97.5 99.3 101.2 102.6 103.4
Peak deceleration, g | 4.13 5.47 7.11 8.81 9.51
Peak heat rate, W/cm2 | 129 144 157 168 191
Peak heat load, kJ/cm2 | 10.4 12.3 16.1 19.7 22.2
Periapsis raise dV, m/s| 22.9 26.1 26.5 27.1 28.2
Apoapsis corr. dV, m/s | -552.78 -101.53 -21.05 28.29 131.42