main.tex

\documentclass[11pt,a4paper]{article}

% Essential packages for aerospace engineering and CFD
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{textcomp}  % For additional text symbols
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{siunitx}
\usepackage{float}
\usepackage[margin=1in]{geometry}
\usepackage{hyperref}

% Python code execution with PythonTeX
\usepackage{pythontex}
\usepackage{fancyvrb}

% Enhanced listings for code display
\usepackage{listings}
\usepackage{xcolor}

% Configure listings for Python
\lstset{
    language=Python,
    basicstyle=\ttfamily\small,
    keywordstyle=\color{blue},
    stringstyle=\color{red},
    commentstyle=\color{gray},
    numberstyle=\tiny,
    numbers=left,
    frame=single,
    breaklines=true
}

% Title and author information
\title{Aerospace Engineering CFD Analysis Template}
\author{Your Name\\
        Department of Aerospace Engineering\\
        Your Institution}
\date{\today}

\begin{document}

\maketitle

\begin{abstract}
This template demonstrates computational fluid dynamics (CFD) analysis techniques in aerospace engineering, including NACA airfoil analysis, aerodynamic performance calculations, propulsion system modeling, and flight dynamics simulations.

The template integrates Python computations with LaTeX for reproducible aerospace engineering reports.
\end{abstract}

\tableofcontents
\newpage

\section{Introduction}

Computational Fluid Dynamics (CFD) has become an essential tool in aerospace engineering for analyzing complex flow phenomena around aircraft, spacecraft, and propulsion systems.

This template provides a foundation for aerospace engineering reports that combine theoretical analysis with computational results.

\section{Aerodynamics Analysis}

\subsection{NACA Airfoil Geometry}

The NACA 4-digit series airfoils are widely used in aerospace applications. For a NACA MPXX airfoil, where M is the maximum camber percentage, P is the position of maximum camber, and XX is the thickness percentage.

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt

def naca_4digit(m, p, t, c=1.0, n=100):
    """
    Generate NACA 4-digit airfoil coordinates
    m: maximum camber (as fraction of chord)
    p: position of maximum camber (as fraction of chord)
    t: thickness (as fraction of chord)
    c: chord length
    n: number of points
    """
    # Cosine spacing for better resolution near leading edge
    beta = np.linspace(0, np.pi, n)
    x = c * (1 - np.cos(beta)) / 2

    # Thickness distribution
    yt = 5 * t * c * (0.2969 * np.sqrt(x/c) -
                      0.1260 * (x/c) -
                      0.3516 * (x/c)**2 +
                      0.2843 * (x/c)**3 -
                      0.1015 * (x/c)**4)

    # Camber line
    yc = np.zeros_like(x)
    dyc_dx = np.zeros_like(x)

    if p > 0:  # Cambered airfoil
        # Forward of maximum camber
        mask1 = x <= p * c
        yc[mask1] = (m * c * x[mask1] / (p * c)**2 *
                     (2 * p * c - x[mask1]) / c**2)
        dyc_dx[mask1] = 2 * m / p**2 * (p - x[mask1]/c)

        # Aft of maximum camber
        mask2 = x > p * c
        yc[mask2] = (m * c * (c - x[mask2]) / (1 - p)**2 / c**2 *
                     (1 + x[mask2]/c - 2*p))
        dyc_dx[mask2] = 2 * m / (1 - p)**2 * (p - x[mask2]/c)

    # Surface coordinates
    theta = np.arctan(dyc_dx)
    xu = x - yt * np.sin(theta)
    yu = yc + yt * np.cos(theta)
    xl = x + yt * np.sin(theta)
    yl = yc - yt * np.cos(theta)

    return xu, yu, xl, yl, x, yc

# Generate NACA 2412 airfoil
xu, yu, xl, yl, x_camber, y_camber = naca_4digit(0.02, 0.4, 0.12)

# Plot airfoil
plt.figure(figsize=(10, 6))
plt.plot(xu, yu, 'b-', linewidth=2, label='Upper surface')
plt.plot(xl, yl, 'r-', linewidth=2, label='Lower surface')
plt.plot(x_camber, y_camber, 'k--', linewidth=1, label='Camber line')
plt.xlabel('x/c')
plt.ylabel('y/c')
plt.title('NACA 2412 Airfoil Geometry')
plt.grid(True, alpha=0.3)
plt.legend()
plt.axis('equal')
plt.tight_layout()
plt.savefig('naca2412_geometry.pdf', dpi=300, bbox_inches='tight')
plt.show()

print(f"NACA 2412 Airfoil Parameters:")
print(f"Maximum camber: 2% at 40% chord")
print(f"Maximum thickness: 12% chord")
print(f"Leading edge radius: {5 * 0.12**2:.4f}c")
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.8\textwidth]{naca2412_geometry.pdf}
\caption{NACA 2412 airfoil geometry showing upper surface, lower surface, and camber line.}
\label{fig:naca2412}
\end{figure}

\subsection{Aerodynamic Performance Analysis}

\begin{pycode}
# Simplified thin airfoil theory calculations
def thin_airfoil_theory(alpha_deg, camber_max, camber_pos):
    """
    Calculate lift coefficient using thin airfoil theory
    alpha_deg: angle of attack in degrees
    camber_max: maximum camber as fraction of chord
    camber_pos: position of maximum camber as fraction of chord
    """
    alpha_rad = np.radians(alpha_deg)

    # Lift curve slope (per radian)
    a0 = 2 * np.pi

    # Zero-lift angle of attack for cambered airfoil
    alpha_L0 = -2 * camber_max * (1 - 2 * camber_pos) / np.pi * np.pi/180

    # Lift coefficient
    cl = a0 * (alpha_rad - alpha_L0)

    return cl, alpha_L0 * 180/np.pi

# Performance analysis for NACA 2412
alpha_range = np.linspace(-5, 15, 21)
cl_values = []
alpha_L0_theory = None

for alpha in alpha_range:
    cl, alpha_L0 = thin_airfoil_theory(alpha, 0.02, 0.4)
    cl_values.append(cl)
    if alpha_L0_theory is None:
        alpha_L0_theory = alpha_L0

cl_values = np.array(cl_values)

# Plot lift curve
plt.figure(figsize=(10, 6))
plt.plot(alpha_range, cl_values, 'bo-', linewidth=2, markersize=6,
         label='Thin Airfoil Theory')
plt.axhline(y=0, color='k', linestyle='--', alpha=0.5)
plt.axvline(x=alpha_L0_theory, color='r', linestyle='--', alpha=0.7,
           label=f'Zero-lift AoA = {alpha_L0_theory:.2f}°')
plt.xlabel('Angle of Attack (degrees)')
plt.ylabel('Lift Coefficient, $C_L$')
plt.title('NACA 2412 Lift Curve - Thin Airfoil Theory')
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig('naca2412_lift_curve.pdf', dpi=300, bbox_inches='tight')
plt.show()

print(f"Aerodynamic Analysis Results:")
print(f"Zero-lift angle of attack: {alpha_L0_theory:.2f} degrees")
print(f"Lift curve slope: {2*np.pi:.3f} per radian "
      f"({2*np.pi*180/np.pi:.1f} per degree)")
print(f"At alpha = 5 degrees: CL = {cl_values[np.argmin(np.abs(alpha_range - 5))]:.3f}")
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.8\textwidth]{naca2412_lift_curve.pdf}
\caption{Lift coefficient variation with angle of attack for NACA 2412 airfoil using thin airfoil theory.}
\label{fig:lift_curve}
\end{figure}

\section{CFD Analysis Fundamentals}

\subsection{Governing Equations}

The Navier-Stokes equations govern fluid flow in aerospace applications:

\begin{align}
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) &= 0 \label{eq:continuity}\\
\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) &= -\nabla p + \nabla \cdot \boldsymbol{\tau} + \rho \mathbf{f} \label{eq:momentum}\\
\frac{\partial (\rho E)}{\partial t} + \nabla \cdot ((\rho E + p)\mathbf{u}) &= \nabla \cdot (\mathbf{k} \nabla T) + \nabla \cdot (\boldsymbol{\tau} \cdot \mathbf{u}) + \rho \mathbf{f} \cdot \mathbf{u} \label{eq:energy}
\end{align}

where $\rho$ is density, $\mathbf{u}$ is velocity vector, $p$ is pressure, $\boldsymbol{\tau}$ is viscous stress tensor, $E$ is total energy per unit mass, and $\mathbf{f}$ represents body forces.

\subsection{Dimensionless Parameters}

Key dimensionless parameters in aerospace CFD:

\begin{pycode}
import pandas as pd

# Create table of important dimensionless parameters
parameters = {
    'Parameter': ['Reynolds Number', 'Mach Number', 'Prandtl Number',
                  'Lift Coefficient', 'Drag Coefficient', 'Pressure Coefficient'],
    'Symbol': ['$Re$', '$Ma$', '$Pr$', '$C_L$', '$C_D$', '$C_p$'],
    'Definition': ['$\\frac{\\rho U L}{\\mu}$', '$\\frac{U}{a}$', '$\\frac{\\mu c_p}{k}$',
                   '$\\frac{L}{\\frac{1}{2}\\rho U^2 S}$', '$\\frac{D}{\\frac{1}{2}\\rho U^2 S}$',
                   '$\\frac{p - p_\\infty}{\\frac{1}{2}\\rho U^2}$'],
    'Physical Meaning': ['Inertial/Viscous forces', 'Flow/Sound speed', 'Momentum/Thermal diffusivity',
                        'Lift/Dynamic pressure', 'Drag/Dynamic pressure', 'Local/Dynamic pressure']
}

df = pd.DataFrame(parameters)
print("Key Dimensionless Parameters in Aerospace CFD:")
print("=" * 40)
for i, row in df.iterrows():
    print(f"{row['Parameter']:<20} {row['Symbol']:<8} {row['Physical Meaning']}")
\end{pycode}

\section{Propulsion Analysis}

\subsection{Turbojet Engine Performance}

\begin{pycode}
def turbojet_performance(Ma0, gamma=1.4, cp=1005, T0=288.15, pi_c=8, pi_b=0.95,
                        T04=1400, eta_c=0.85, eta_t=0.88, eta_n=0.98):
    """
    Simplified turbojet engine performance analysis
    Ma0: flight Mach number
    gamma: specific heat ratio
    cp: specific heat at constant pressure (J/kg/K)
    T0: ambient temperature (K)
    pi_c: compressor pressure ratio
    pi_b: burner pressure ratio
    T04: turbine inlet temperature (K)
    eta_c: compressor efficiency
    eta_t: turbine efficiency
    eta_n: nozzle efficiency
    """
    R = cp * (gamma - 1) / gamma

    # Station 0: Ambient conditions
    T00 = T0 * (1 + (gamma - 1) / 2 * Ma0**2)
    p00 = 101325 * (T00 / T0)**(gamma / (gamma - 1))

    # Station 2: Compressor exit (ideal)
    T02_ideal = T00 * pi_c**((gamma - 1) / gamma)
    T02 = T00 + (T02_ideal - T00) / eta_c
    p02 = p00 * pi_c

    # Station 3: Burner exit
    T03 = T04  # Assuming complete combustion
    p03 = p02 * pi_b

    # Station 4: Turbine exit
    # Work balance: turbine work = compressor work
    T04_ideal = T03 - (T02 - T00) / eta_t
    T04_actual = T03 - (T02 - T00) / eta_t
    pi_t = (T04_actual / T03)**(gamma / (gamma - 1))
    p04 = p03 * pi_t

    # Station 9: Nozzle exit
    p9 = 101325  # Assume perfectly expanded
    T09_ideal = T04_actual * (p9 / p04)**((gamma - 1) / gamma)
    T9 = T04_actual - eta_n * (T04_actual - T09_ideal)

    # Exhaust velocity
    V9 = np.sqrt(2 * cp * (T04_actual - T9))

    # Flight velocity
    V0 = Ma0 * np.sqrt(gamma * R * T0)

    # Specific thrust and specific fuel consumption (simplified)
    F_specific = V9 - V0  # N per kg/s of air

    return {
        'V0': V0,
        'V9': V9,
        'F_specific': F_specific,
        'T02': T02,
        'T04': T04_actual,
        'pi_c': pi_c,
        'pi_t': pi_t
    }

# Performance analysis for different flight conditions
mach_numbers = np.linspace(0, 2.0, 21)
results = []

for Ma in mach_numbers:
    perf = turbojet_performance(Ma)
    results.append([Ma, perf['F_specific'], perf['V9'], perf['V0']])

results = np.array(results)

# Plot performance
plt.figure(figsize=(12, 8))

plt.subplot(2, 2, 1)
plt.plot(results[:, 0], results[:, 1], 'b-o', linewidth=2, markersize=4)
plt.xlabel('Flight Mach Number')
plt.ylabel('Specific Thrust (N·s/kg)')
plt.title('Turbojet Specific Thrust vs Flight Mach')
plt.grid(True, alpha=0.3)

plt.subplot(2, 2, 2)
plt.plot(results[:, 0], results[:, 2], 'r-', linewidth=2, label='Exhaust Velocity')
plt.plot(results[:, 0], results[:, 3], 'g--', linewidth=2, label='Flight Velocity')
plt.xlabel('Flight Mach Number')
plt.ylabel('Velocity (m/s)')
plt.title('Velocities vs Flight Mach')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(2, 2, 3)
efficiency = results[:, 1] * results[:, 3] / (0.5 * (results[:, 2]**2 - results[:, 3]**2))
plt.plot(results[:, 0], efficiency, 'purple', linewidth=2)
plt.xlabel('Flight Mach Number')
plt.ylabel('Propulsive Efficiency')
plt.title('Propulsive Efficiency vs Flight Mach')
plt.grid(True, alpha=0.3)

plt.subplot(2, 2, 4)
thrust_ratio = results[:, 1] / results[0, 1]  # Normalized to Ma=0
plt.plot(results[:, 0], thrust_ratio, 'orange', linewidth=2)
plt.xlabel('Flight Mach Number')
plt.ylabel('Normalized Specific Thrust')
plt.title('Thrust Variation with Altitude')
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('turbojet_performance.pdf', dpi=300, bbox_inches='tight')
plt.show()

print("Turbojet Performance Summary:")
print(f"Sea level static specific thrust: {results[0, 1]:.1f} N·s/kg")
print(f"Cruise (Ma=0.8) specific thrust: {results[16, 1]:.1f} N·s/kg")
print(f"Maximum analyzed Mach: {results[-1, 0]:.1f}")
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{turbojet_performance.pdf}
\caption{Turbojet engine performance characteristics showing specific thrust, velocities, propulsive efficiency, and thrust variation with flight Mach number.}
\label{fig:turbojet}
\end{figure}

\section{Flight Dynamics}

\subsection{Aircraft Stability Analysis}

\begin{pycode}
def static_stability_analysis(cg_position, ac_position, cl_alpha, cm_alpha_wing,
                             cm_alpha_tail, tail_volume):
    """
    Analyze longitudinal static stability
    cg_position: center of gravity position (fraction of MAC)
    ac_position: aerodynamic center position (fraction of MAC)
    cl_alpha: lift curve slope (per radian)
    cm_alpha_wing: wing pitching moment curve slope
    cm_alpha_tail: tail pitching moment curve slope
    tail_volume: tail volume coefficient
    """

    # Static margin
    static_margin = ac_position - cg_position

    # Total pitching moment curve slope
    cm_alpha_total = cm_alpha_wing + cm_alpha_tail * tail_volume

    # Neutral point
    neutral_point = ac_position - cm_alpha_total / cl_alpha

    # Stability criterion
    is_stable = static_margin > 0

    return {
        'static_margin': static_margin,
        'neutral_point': neutral_point,
        'cm_alpha_total': cm_alpha_total,
        'is_stable': is_stable
    }

# Example aircraft stability analysis
aircraft_data = {
    'cg_range': np.linspace(0.20, 0.35, 100),  # CG positions from 20% to 35% MAC
    'ac_position': 0.25,  # Aerodynamic center at 25% MAC
    'cl_alpha': 2 * np.pi,  # Lift curve slope
    'cm_alpha_wing': -0.1,  # Wing contribution
    'cm_alpha_tail': -0.8,  # Tail contribution
    'tail_volume': 0.5  # Tail volume coefficient
}

stability_results = []
for cg in aircraft_data['cg_range']:
    result = static_stability_analysis(
        cg, aircraft_data['ac_position'], aircraft_data['cl_alpha'],
        aircraft_data['cm_alpha_wing'], aircraft_data['cm_alpha_tail'],
        aircraft_data['tail_volume']
    )
    stability_results.append([cg, result['static_margin'], result['is_stable']])

stability_results = np.array(stability_results)

# Plot stability analysis
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(stability_results[:, 0] * 100, stability_results[:, 1] * 100,
         'b-', linewidth=2)
plt.axhline(y=0, color='r', linestyle='--', linewidth=2, label='Neutral Stability')
stable_mask = stability_results[:, 1] > 0
plt.fill_between(stability_results[stable_mask, 0] * 100,
                 stability_results[stable_mask, 1] * 100,
                 alpha=0.3, color='green', label='Stable Region')
plt.xlabel('CG Position (% MAC)')
plt.ylabel('Static Margin (% MAC)')
plt.title('Longitudinal Static Stability')
plt.grid(True, alpha=0.3)
plt.legend()

plt.subplot(1, 2, 2)
# Control authority analysis (simplified)
elevator_deflection = np.linspace(-20, 20, 41)
cm_delta_e = -0.02  # Elevator effectiveness (per degree)
trim_moments = []

for delta_e in elevator_deflection:
    cm_trim = cm_delta_e * delta_e
    trim_moments.append(cm_trim)

plt.plot(elevator_deflection, trim_moments, 'g-', linewidth=2,
         label='Elevator Moment')
plt.axhline(y=0, color='k', linestyle='-', alpha=0.5)
plt.xlabel('Elevator Deflection (degrees)')
plt.ylabel('Pitching Moment Coefficient')
plt.title('Control Authority')
plt.grid(True, alpha=0.3)
plt.legend()

plt.tight_layout()
plt.savefig('flight_dynamics_stability.pdf', dpi=300, bbox_inches='tight')
plt.show()

print("Flight Dynamics Analysis:")
print(f"Aerodynamic center: {aircraft_data['ac_position']*100:.1f}% MAC")
print(f"Stable CG range: {stability_results[stable_mask, 0].min()*100:.1f}% to {stability_results[stable_mask, 0].max()*100:.1f}% MAC")
print(f"Maximum static margin: {stability_results[:, 1].max()*100:.1f}% MAC")
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{flight_dynamics_stability.pdf}
\caption{Longitudinal static stability analysis showing static margin variation with CG position and elevator control authority.}
\label{fig:stability}
\end{figure}

\section{Computational Methods}

\subsection{Finite Difference Discretization}

For CFD applications, the governing equations are discretized using finite difference, finite volume, or finite element methods. A simple example of finite difference discretization for the 1D heat equation:

\begin{equation}
\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}
\end{equation}

\begin{pycode}
def heat_equation_1d(nx=50, nt=100, dx=0.02, dt=0.001, alpha=0.01):
    """
    Solve 1D heat equation using finite differences
    Demonstrates numerical methods used in CFD
    """
    x = np.linspace(0, 1, nx)
    T = np.zeros((nt, nx))

    # Initial condition: Gaussian temperature distribution
    T[0, :] = np.exp(-50 * (x - 0.5)**2)

    # Boundary conditions (Dirichlet)
    T[:, 0] = 0
    T[:, -1] = 0

    # Time stepping using explicit finite differences
    r = alpha * dt / dx**2  # Stability parameter

    for n in range(1, nt):
        for i in range(1, nx-1):
            T[n, i] = T[n-1, i] + r * (T[n-1, i+1] - 2*T[n-1, i] + T[n-1, i-1])

    return x, T

# Solve heat equation
x, T_solution = heat_equation_1d()

# Plot temperature evolution
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
time_steps = [0, 20, 50, 99]
for t in time_steps:
    plt.plot(x, T_solution[t, :], label=f't = {t*0.001:.3f}s')
plt.xlabel('Position (x)')
plt.ylabel('Temperature')
plt.title('1D Heat Equation Solution')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(1, 2, 2)
X, T_time = np.meshgrid(x, np.linspace(0, 0.099, 100))
plt.contourf(X, T_time, T_solution, levels=20, cmap='hot')
plt.colorbar(label='Temperature')
plt.xlabel('Position (x)')
plt.ylabel('Time (s)')
plt.title('Temperature Evolution')

plt.tight_layout()
plt.savefig('heat_equation_solution.pdf', dpi=300, bbox_inches='tight')
plt.show()

print("Numerical Method Demonstration:")
print(f"Grid points: {len(x)}")
print(f"Time steps: {T_solution.shape[0]}")
print(f"Stability parameter r = {0.01 * 0.001 / 0.02**2:.3f} (should be <= 0.5)")
print(f"Maximum temperature at t=0: {T_solution[0, :].max():.3f}")
print(f"Maximum temperature at final time: {T_solution[-1, :].max():.3f}")
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{heat_equation_solution.pdf}
\caption{Solution of the 1D heat equation demonstrating finite difference methods commonly used in CFD applications.}
\label{fig:heat_equation}
\end{figure}

\section{Conclusion}

This template demonstrates the integration of theoretical aerospace engineering concepts with computational analysis using Python and LaTeX. The combination of:

\begin{itemize}
\item NACA airfoil geometry generation and analysis
\item Aerodynamic performance calculations using thin airfoil theory
\item Turbojet engine performance modeling
\item Flight dynamics and stability analysis
\item Numerical methods for CFD applications
\end{itemize}

provides a comprehensive foundation for aerospace engineering reports that require both analytical and computational approaches. The PythonTeX integration allows for reproducible results and seamless combination of code, calculations, and professional documentation.

\section{References}

\begin{enumerate}
\item Anderson, J.D., \textit{Fundamentals of Aerodynamics}, 6th Edition, McGraw-Hill, 2017.
\item Mattingly, J.D., \textit{Elements of Propulsion: Gas Turbines and Rockets}, AIAA, 2006.
\item Etkin, B. and Reid, L.D., \textit{Dynamics of Flight: Stability and Control}, 3rd Edition, Wiley, 1996.
\item Blazek, J., \textit{Computational Fluid Dynamics: Principles and Applications}, 3rd Edition, Butterworth-Heinemann, 2015.
\end{enumerate}

\end{document}