rocket_propulsion.tex

\documentclass[a4paper, 11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{booktabs}
\usepackage{algorithm2e}
\usepackage{subcaption}
\usepackage[makestderr]{pythontex}

% Theorem environments for lab report style
\newtheorem{definition}{Definition}
\newtheorem{theorem}{Theorem}
\newtheorem{remark}{Remark}

\title{Rocket Propulsion Analysis: Thrust Curves, Specific Impulse, and Staging Optimization\\
\large A Comprehensive Study of Chemical Rocket Performance}
\author{Propulsion Systems Division\\Computational Science Templates}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
This laboratory report presents a comprehensive analysis of rocket propulsion systems. We examine thrust curves for different propellant combinations, compare specific impulse values, and optimize multi-stage rocket configurations using the Tsiolkovsky equation. The analysis includes propellant mass flow rates, chamber pressure effects, and payload fraction optimization for orbital insertion missions.
\end{abstract}

\section{Objectives}
\begin{enumerate}
    \item Analyze thrust and specific impulse for various propellant combinations
    \item Compare single-stage and multi-stage rocket performance
    \item Optimize staging ratios for maximum payload fraction
    \item Evaluate thrust-to-weight ratios for different mission profiles
\end{enumerate}

\section{Theoretical Background}

\begin{definition}[Specific Impulse]
Specific impulse is the total impulse per unit weight of propellant:
\begin{equation}
I_{sp} = \frac{F}{\dot{m} g_0} = \frac{v_e}{g_0}
\end{equation}
where $F$ is thrust, $\dot{m}$ is mass flow rate, and $v_e$ is effective exhaust velocity.
\end{definition}

\subsection{Tsiolkovsky Rocket Equation}
The ideal velocity change achievable by a rocket:
\begin{equation}
\Delta v = v_e \ln\left(\frac{m_0}{m_f}\right) = I_{sp} g_0 \ln(MR)
\end{equation}
where $MR = m_0/m_f$ is the mass ratio.

\subsection{Thrust Equation}
\begin{theorem}[Rocket Thrust]
The thrust generated by a rocket engine:
\begin{equation}
F = \dot{m} v_e + (p_e - p_a) A_e
\end{equation}
where $p_e$ is exit pressure, $p_a$ is ambient pressure, and $A_e$ is exit area.
\end{theorem}

\subsection{Staging Analysis}
For an $n$-stage rocket with equal structural coefficients:
\begin{equation}
\lambda_{payload} = \left[\frac{1 - \epsilon \cdot MR_{stage}}{MR_{stage}}\right]^n
\end{equation}
where $\epsilon$ is the structural coefficient (typically 0.05-0.15).

\section{Computational Analysis}

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar
plt.rc('text', usetex=True)
plt.rc('font', family='serif')

np.random.seed(42)

g0 = 9.81  # m/s^2

# Propellant data: Isp (s), density (kg/m^3), cost index
propellants = {
    'Solid (APCP)': {'Isp': 265, 'density': 1800, 'T_chamber': 3400},
    'RP-1/LOX': {'Isp': 353, 'density': 1030, 'T_chamber': 3670},
    'LH2/LOX': {'Isp': 455, 'density': 360, 'T_chamber': 3250},
    'MMH/N2O4': {'Isp': 340, 'density': 1200, 'T_chamber': 3000},
    'CH4/LOX': {'Isp': 380, 'density': 830, 'T_chamber': 3550}
}

# Mission delta-v requirements (m/s)
missions = {
    'LEO': 9400,
    'GTO': 13500,
    'TLI': 15500,
    'Mars Transfer': 18000
}

# Single-stage payload fraction
def single_stage_payload(dv, Isp, epsilon=0.10):
    ve = Isp * g0
    MR = np.exp(dv / ve)
    if MR <= 1:
        return 0
    payload_frac = (1 - epsilon * MR) / MR
    return max(0, payload_frac)

# Multi-stage payload fraction
def multi_stage_payload(dv, Isp, n_stages, epsilon=0.08):
    ve = Isp * g0
    dv_per_stage = dv / n_stages
    MR_stage = np.exp(dv_per_stage / ve)

    if MR_stage <= 1:
        return 0

    stage_payload = (1 - epsilon * MR_stage) / MR_stage
    if stage_payload <= 0:
        return 0

    return stage_payload ** n_stages

# Thrust curve profiles
def thrust_curve(t, F_max, t_burn, profile='constant'):
    thrust = np.zeros_like(t)
    burning = t <= t_burn

    if profile == 'constant':
        thrust[burning] = F_max
    elif profile == 'regressive':
        thrust[burning] = F_max * (1 - 0.3 * t[burning]/t_burn)
    elif profile == 'progressive':
        thrust[burning] = F_max * (0.7 + 0.3 * t[burning]/t_burn)
    elif profile == 'boost-sustain':
        boost = t <= t_burn * 0.3
        sustain = (t > t_burn * 0.3) & (t <= t_burn)
        thrust[boost] = F_max
        thrust[sustain] = F_max * 0.5

    return thrust

# Time array
t = np.linspace(0, 150, 1000)
F_max = 2e6  # 2 MN
t_burn = 120  # s

# Generate thrust curves
profiles = ['constant', 'regressive', 'progressive', 'boost-sustain']
thrust_curves = {p: thrust_curve(t, F_max, t_burn, p) for p in profiles}

# Calculate total impulse for each profile
total_impulse = {p: np.trapz(thrust_curves[p], t) for p in profiles}

# Staging analysis
n_stages_range = range(1, 6)
staging_results = {}
for mission, dv in missions.items():
    staging_results[mission] = []
    for n in n_stages_range:
        pf = multi_stage_payload(dv, 455, n, 0.08)  # LH2/LOX
        staging_results[mission].append(pf * 100)

# Optimal staging ratio
def optimal_stage_ratio(dv_total, Isp, n_stages, epsilon):
    """Find optimal mass ratio per stage."""
    ve = Isp * g0

    def neg_payload(MR):
        if MR <= 1 or MR > 20:
            return 1e10
        dv_achieved = n_stages * ve * np.log(MR)
        if dv_achieved < dv_total:
            return 1e10
        stage_pf = (1 - epsilon * MR) / MR
        if stage_pf <= 0:
            return 1e10
        return -stage_pf ** n_stages

    result = minimize_scalar(neg_payload, bounds=(1.5, 15), method='bounded')
    return result.x, -result.fun

# Create comprehensive visualization
fig = plt.figure(figsize=(14, 12))

# Plot 1: Thrust curves
ax1 = fig.add_subplot(2, 3, 1)
colors = ['blue', 'red', 'green', 'purple']
for prof, color in zip(profiles, colors):
    ax1.plot(t, thrust_curves[prof]/1e6, color=color, linewidth=2,
             label=prof.replace('-', ' ').title())
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Thrust (MN)')
ax1.set_title('Thrust Curve Profiles')
ax1.legend(fontsize=8)
ax1.grid(True, alpha=0.3)

# Plot 2: Propellant comparison
ax2 = fig.add_subplot(2, 3, 2)
names = list(propellants.keys())
isps = [propellants[n]['Isp'] for n in names]
densities = [propellants[n]['density'] for n in names]

x = np.arange(len(names))
width = 0.35
bars1 = ax2.bar(x - width/2, isps, width, label='$I_{sp}$ (s)', color='steelblue')
ax2_twin = ax2.twinx()
bars2 = ax2_twin.bar(x + width/2, densities, width, label='$\\rho$ (kg/m$^3$)', color='coral')
ax2.set_xlabel('Propellant')
ax2.set_ylabel('$I_{sp}$ (s)', color='steelblue')
ax2_twin.set_ylabel('Density (kg/m$^3$)', color='coral')
ax2.set_xticks(x)
ax2.set_xticklabels([n.split()[0] for n in names], rotation=45, ha='right')
ax2.set_title('Propellant Performance')
ax2.legend(loc='upper left', fontsize=7)
ax2_twin.legend(loc='upper right', fontsize=7)

# Plot 3: Single-stage performance
ax3 = fig.add_subplot(2, 3, 3)
dv_range = np.linspace(1000, 12000, 100)
for name in ['RP-1/LOX', 'LH2/LOX', 'CH4/LOX']:
    pf = [single_stage_payload(dv, propellants[name]['Isp']) * 100 for dv in dv_range]
    ax3.plot(dv_range/1000, pf, linewidth=2, label=name)
ax3.axvline(x=9.4, color='gray', linestyle='--', alpha=0.7)
ax3.text(9.5, 12, 'LEO', fontsize=8)
ax3.set_xlabel(r'$\Delta v$ (km/s)')
ax3.set_ylabel('Payload Fraction (\\%)')
ax3.set_title('Single-Stage Performance')
ax3.legend(fontsize=8)
ax3.grid(True, alpha=0.3)
ax3.set_ylim([0, 25])

# Plot 4: Staging benefits
ax4 = fig.add_subplot(2, 3, 4)
markers = ['o', 's', '^', 'd']
for (mission, results), marker in zip(staging_results.items(), markers):
    if mission in ['LEO', 'GTO', 'Mars Transfer']:
        ax4.plot(list(n_stages_range), results, marker=marker,
                 linewidth=2, markersize=6, label=mission)
ax4.set_xlabel('Number of Stages')
ax4.set_ylabel('Payload Fraction (\\%)')
ax4.set_title('Staging Optimization (LH2/LOX)')
ax4.legend(fontsize=8)
ax4.grid(True, alpha=0.3)
ax4.set_xticks(list(n_stages_range))

# Plot 5: Mass ratio requirements
ax5 = fig.add_subplot(2, 3, 5)
Isp_range = np.linspace(250, 500, 100)
for mission in ['LEO', 'GTO', 'TLI']:
    MR = np.exp(missions[mission] / (Isp_range * g0))
    ax5.plot(Isp_range, MR, linewidth=2, label=mission)
ax5.set_xlabel('$I_{sp}$ (s)')
ax5.set_ylabel('Required Mass Ratio')
ax5.set_title('Mass Ratio vs Specific Impulse')
ax5.legend(fontsize=8)
ax5.grid(True, alpha=0.3)
ax5.set_ylim([0, 100])

# Plot 6: Delta-v budget breakdown
ax6 = fig.add_subplot(2, 3, 6)
mission_names = list(missions.keys())
dvs = [missions[m]/1000 for m in mission_names]
colors = plt.cm.viridis(np.linspace(0.2, 0.8, len(missions)))
bars = ax6.bar(mission_names, dvs, color=colors)
ax6.set_xlabel('Mission')
ax6.set_ylabel(r'$\Delta v$ (km/s)')
ax6.set_title('Mission Requirements')
for bar, dv in zip(bars, dvs):
    ax6.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.3,
             f'{dv:.1f}', ha='center', va='bottom', fontsize=9)
ax6.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('rocket_propulsion_plot.pdf', bbox_inches='tight', dpi=150)
print(r'\begin{center}')
print(r'\includegraphics[width=\textwidth]{rocket_propulsion_plot.pdf}')
print(r'\end{center}')
plt.close()

# Key results
best_2stage = staging_results['LEO'][1]
best_3stage = staging_results['LEO'][2]
opt_MR, opt_pf = optimal_stage_ratio(9400, 455, 2, 0.08)
\end{pycode}

\section{Algorithm}

\begin{algorithm}[H]
\SetAlgoLined
\KwIn{Mission $\Delta v$, propellant $I_{sp}$, number of stages $n$, structural coefficient $\epsilon$}
\KwOut{Payload fraction $\lambda$}
$v_e \leftarrow I_{sp} \cdot g_0$\;
$\Delta v_{stage} \leftarrow \Delta v / n$\;
$MR_{stage} \leftarrow \exp(\Delta v_{stage} / v_e)$\;
$\lambda_{stage} \leftarrow (1 - \epsilon \cdot MR_{stage}) / MR_{stage}$\;
$\lambda \leftarrow \lambda_{stage}^n$\;
\Return{$\lambda$}
\caption{Multi-Stage Payload Fraction Calculation}
\end{algorithm}

\section{Results and Discussion}

\subsection{Propellant Performance}

\begin{pycode}
print(r'\begin{table}[h]')
print(r'\centering')
print(r'\caption{Propellant Performance Comparison}')
print(r'\begin{tabular}{lccc}')
print(r'\toprule')
print(r'Propellant & $I_{sp}$ (s) & Density (kg/m$^3$) & $\rho \cdot I_{sp}$ \\')
print(r'\midrule')
for name, data in propellants.items():
    density_isp = data['Isp'] * data['density'] / 1000
    print(f"{name} & {data['Isp']} & {data['density']} & {density_isp:.0f} \\\\")
print(r'\bottomrule')
print(r'\end{tabular}')
print(r'\end{table}')
\end{pycode}

\begin{remark}[Propellant Selection Trade-offs]
LH2/LOX provides the highest $I_{sp}$ (455 s) but lowest density, requiring larger tanks. RP-1/LOX offers good performance with higher density, making it preferred for first stages. CH4/LOX (Methalox) is gaining popularity for its balance of performance, storability, and potential for in-situ resource utilization on Mars.
\end{remark}

\subsection{Staging Benefits}

For LEO insertion ($\Delta v = 9.4$ km/s) with LH2/LOX propellant:
\begin{itemize}
    \item Single stage: \py{f"{staging_results['LEO'][0]:.2f}"}\% payload fraction
    \item Two stages: \py{f"{best_2stage:.2f}"}\% payload fraction
    \item Three stages: \py{f"{best_3stage:.2f}"}\% payload fraction
\end{itemize}

Optimal mass ratio for 2-stage LEO: $MR = $ \py{f"{opt_MR:.2f}"} yielding \py{f"{opt_pf*100:.2f}"}\% payload.

\subsection{Thrust Profile Analysis}

\begin{pycode}
print(r'\begin{table}[h]')
print(r'\centering')
print(r'\caption{Thrust Profile Total Impulse Comparison}')
print(r'\begin{tabular}{lc}')
print(r'\toprule')
print(r'Profile & Total Impulse (MN$\cdot$s) \\')
print(r'\midrule')
for prof in profiles:
    I_total = total_impulse[prof] / 1e6
    print(f"{prof.replace('-', ' ').title()} & {I_total:.1f} \\\\")
print(r'\bottomrule')
print(r'\end{tabular}')
print(r'\end{table}')
\end{pycode}

\begin{remark}[Thrust Profile Selection]
Regressive profiles (common in solid rockets) provide higher initial thrust for liftoff, while boost-sustain profiles optimize gravity losses. Progressive profiles are rare but can reduce initial structural loads.
\end{remark}

\section{Limitations and Extensions}

\subsection{Model Limitations}
\begin{enumerate}
    \item \textbf{Ideal rocket}: Neglects nozzle losses, incomplete combustion
    \item \textbf{Constant $I_{sp}$}: Real engines vary with altitude
    \item \textbf{Gravity/drag losses}: Not included in $\Delta v$ budget
    \item \textbf{Fixed structural coefficient}: Varies with stage size
\end{enumerate}

\subsection{Possible Extensions}
\begin{itemize}
    \item Trajectory optimization with gravity and drag
    \item Parallel staging (boosters) analysis
    \item Reusability impact on payload fraction
    \item Electric propulsion for upper stages
\end{itemize}

\section{Conclusions}
\begin{itemize}
    \item LH2/LOX provides best $I_{sp}$ (455 s) for upper stages
    \item Staging dramatically improves payload fraction (factor of 2-3)
    \item Optimal number of stages is 2-3 for most Earth-orbit missions
    \item Propellant density matters for first stages (tank mass)
    \item Modern trends favor CH4/LOX for reusability and ISRU
\end{itemize}

\section*{References}
\begin{itemize}
    \item Sutton, G. P., \& Biblarz, O. (2016). \textit{Rocket Propulsion Elements}. Wiley.
    \item Turner, M. J. L. (2008). \textit{Rocket and Spacecraft Propulsion}. Springer.
    \item Humble, R. W., Henry, G. N., \& Larson, W. J. (1995). \textit{Space Propulsion Analysis and Design}. McGraw-Hill.
\end{itemize}

\end{document}