THE ENGINEERING OF ELECTRICAL VERTICAL TAKEOFF: A deep technical study

Every eVTOL conceals the same brutal physical contradiction. Hovering demands enormous disc area and low disc loading to minimise induced power. Efficient cruise demands a compact, streamlined vehicle. You cannot fully satisfy both. Every design choice in an electric vertical-takeoff aircraft is a negotiation with this paradox — and understanding the negotiation is understanding the machine.

The Physics of Hover
Three Architectures, One Contradiction
Rotor Aerodynamics — The Blade-Level Story
The Battery Constraint and Electric Range
Flight Control — Why Linear Theory Breaks Down
Acoustic Engineering — The Whisper Problem
Structural Philosophy — Why Carbon Is Not Enough
The Cascade — How Everything Connects

1. The Physics of Hover

The foundation is actuator disc theory — a control-volume analysis that treats the rotor as an idealised pressure-jumping disc. Air enters from infinity above the disc at freestream velocity, accelerates through it, and exits below in a fully-developed slipstream. The momentum equation gives the thrust; the energy equation gives the minimum possible power to produce it.

1.1 Disc Loading

Define disc loading DL as thrust T per unit total disc area A_total. This single number dictates nearly every downstream performance characteristic of the vehicle.

DL = T / A_total  =  T / (n · π · R²)

Symbol	Meaning	Units
`DL`	Disc loading	N/m²
`T`	Total thrust (= weight at hover)	N
`n`	Number of rotors	—
`R`	Rotor radius	m

At hover, T = W (total vehicle weight in Newtons). The higher the disc loading, the harder each unit of disc area has to work — and the more power it consumes doing so.

1.2 Induced Velocity

The induced velocity through the disc — the velocity imparted to the airstream to generate momentum — follows directly from the actuator disc momentum equation:

v_i  =  √( T / (2ρA) )  =  √( DL / 2ρ )

Where ρ = 1.225 kg/m³ at ISA sea level. Notice: v_i scales with the square root of disc loading. Doubling DL increases induced velocity by only 41% — but as we will see, power scales more steeply.

1.3 Hover Power and Figure of Merit

Ideal hover power is thrust times induced velocity. Actual shaft power is penalised by the Figure of Merit FM (typically 0.70–0.82 for a well-designed rotor). Electrical power further accounts for motor and inverter efficiency η:

P_ideal  =  T · v_i  =  T^(3/2) / √(2ρA)         ... (Eq. 1)

P_shaft  =  P_ideal / FM                            ... (Eq. 2)

P_electric  =  P_shaft / η                          ... (Eq. 3)

The T^(3/2) dependence is the crucial insight. Rewriting in terms of disc loading:

P_ideal  =  (DL)^(1/2) · T / √(2ρ)

Halving disc loading cuts hover power by approximately 30% — a powerful driver toward large-diameter rotors, entirely independent of motor or battery technology. Physics, not engineering, sets this floor.

Key numbers — disc loading benchmarks:

Configuration DL (N/m²)

Conventional helicopter 48 – 120

Lift+cruise eVTOL (hover mode) 150 – 350

Tiltrotor eVTOL 200 – 500

Multirotor eVTOL 300 – 800

Configuration	DL (N/m²)
Conventional helicopter	48 – 120
Lift+cruise eVTOL (hover mode)	150 – 350
Tiltrotor eVTOL	200 – 500
Multirotor eVTOL	300 – 800

1.4 Python: Hover Power Solver

import numpy as np

def hover_power(W_N, n, R, FM=0.76, eta=0.93, rho=1.225):
    """
    Actuator disc hover power solver.

    Parameters
    ----------
    W_N  : vehicle weight in Newtons
    n    : number of rotors
    R    : rotor radius (m)
    FM   : figure of merit (0.70 – 0.85)
    eta  : motor + drive efficiency
    rho  : air density (kg/m³)

    Returns
    -------
    dict of disc loading, induced velocity, and power components
    """
    A       = n * np.pi * R**2          # total disc area  [m²]
    DL      = W_N / A                   # disc loading     [N/m²]
    v_i     = np.sqrt(DL / (2 * rho))   # induced velocity [m/s]
    P_ideal = W_N * v_i                 # ideal power      [W]
    P_shaft = P_ideal / FM              # shaft power      [W]
    P_elec  = P_shaft / eta             # electrical draw  [W]
    Omega   = v_i / R                   # approx hover angular velocity [rad/s]
    RPM     = Omega * 60 / (2 * np.pi)

    return {
        "DL_Nm2"      : round(DL, 1),
        "v_i_ms"      : round(v_i, 2),
        "P_ideal_kW"  : round(P_ideal / 1e3, 2),
        "P_shaft_kW"  : round(P_shaft / 1e3, 2),
        "P_elec_kW"   : round(P_elec / 1e3, 2),
        "RPM_estimate": round(RPM, 0)
    }


# Example: 800 kg, 6× rotors of 1.0 m radius
result = hover_power(800 * 9.81, n=6, R=1.0)
print(result)

# Output:
# DL = 409 N/m², v_i = 12.9 m/s,
# P_ideal = 101.6 kW, P_elec ≈ 145 kW

# Sensitivity study — rotor radius effect
print("\nRadius sensitivity (n=6, W=800 kg):")
for R in [0.8, 1.0, 1.2, 1.4, 1.6]:
    r = hover_power(800 * 9.81, n=6, R=R)
    print(f"  R={R:.1f}m → DL={r['DL_Nm2']:.0f} N/m²,"
          f" P_elec={r['P_elec_kW']:.1f} kW")

Output:

  R=0.8m → DL=639 N/m²,  P_elec=183.4 kW
  R=1.0m → DL=409 N/m²,  P_elec=144.6 kW
  R=1.2m → DL=284 N/m²,  P_elec=117.9 kW
  R=1.4m → DL=208 N/m²,  P_elec=98.8  kW
  R=1.6m → DL=159 N/m²,  P_elec=84.2  kW

Each 20% increase in rotor radius reduces electrical hover power by roughly 16% — entirely from the disc loading reduction. No change in motors, batteries, or efficiency. The argument for large-diameter rotors is purely aerodynamic.

2. Three Architectures, One Contradiction

Three structural answers have emerged to the hover-cruise paradox. They differ not in goal but in how aggressively they trade hover efficiency against cruise capability.

2.1 Architecture Schematics (ASCII Diagrams)

MULTIROTOR
─────────────────────────────────────────────────────────
     ○ ─── ╔═══╗ ─── ○
           ║   ║
     ○ ─── ╚═══╝ ─── ○
     ↑ all rotors provide   ↑
       lift AND control

  DL: 300–600 N/m²   |   Cruise L/D: 4–7
  Pros: Simple, redundant, no moving mechanisms
  Cons: Poor cruise efficiency, all rotors spinning at cruise


LIFT + CRUISE
─────────────────────────────────────────────────────────
     ○       ○
     │  ╔══════════╗  │
─────┤  ║  cabin   ║  ├─────────► [pusher/puller prop]
     │  ╚══════════╝  │
     ○       ○
     ↑ lift rotors    ↑            ↑ cruise propulsor
       (fold/stop                    (separate)
        in cruise)

  DL: 150–350 N/m²   |   Cruise L/D: 10–16
  Pros: Excellent cruise L/D, each system optimised independently
  Cons: Lift rotor deadweight at cruise, higher structural mass


TILTROTOR / TILT-WING
─────────────────────────────────────────────────────────

  HOVER mode:         CRUISE mode:
     ↑    ↑              ▶         ▶
     ○    ○              ○         ○
     │────│────╔══╗───────────────
                ║  ║
     ──────────────────────
     ○    ○              ○         ○
     ↑    ↑              ▶         ▶

  DL: 200–500 N/m²   |   Cruise L/D: 12–18
  Pros: Best cruise efficiency of all three, same rotor serves both
  Cons: Complex tilt mechanism, dynamic loads in transition corridor

2.2 The Fundamental Trade-off

Metric	Multirotor	Lift + Cruise	Tiltrotor
Hover power efficiency	Medium	High	Medium-High
Cruise aerodynamic efficiency	Low	High	Highest
Mechanical complexity	Low	Medium	High
Certification risk	Low	Medium	High
Range potential	Short	Medium	Longest
Structural mass fraction	Lowest	Highest	Medium
Transition corridor risk	None	Low	High

The multirotor accepts its hover-efficiency penalty by keeping rotor radii large relative to the fuselage and eliminating all mechanical complexity. The lift+cruise architecture completely segregates hover and cruise functions — the aerodynamic design of each system is unconstrained by the other's requirements, but the vehicle carries both systems as deadweight throughout the entire mission. The tiltrotor solves the problem mechanically, achieving the best cruise efficiency but demanding high-load tilt actuators and a carefully managed transition corridor where neither pure hover nor pure wing-borne flight applies.

3. Rotor Aerodynamics — The Blade-Level Story

3.1 Figure of Merit

FM is to a rotor what L/D is to a wing — a single-number quality indicator. A rotor hovering at FM = 0.80 is producing ideal thrust while losing 25% additional power to profile drag, tip losses, swirl, and non-uniform inflow.

FM  =  P_ideal / P_actual
    =  (T^(3/2) / √(2ρA)) / P_actual        ... (Eq. 4)

Values above 0.82 require near-optimal twist distribution, low tip Mach number, and careful planform design. Peak FM occurs near the blade's aerodynamic optimum — operating far from it degrades efficiency in both directions.

3.2 Tip Mach Number — The Hard Ceiling

Compressibility becomes significant above M_tip ≈ 0.60 and catastrophic above 0.75. Shock formation on the outer blade segment collapses local L/D, dramatically increases profile power, and is the dominant source of impulsive noise. Tip speed is therefore constrained to approximately 150–220 m/s for urban eVTOL designs.

M_tip  =  (Ω · R  +  V_fwd · cos α)  /  c_∞        ... (Eq. 5)

Where Ω is rotor angular velocity, V_fwd is forward flight speed, α is the rotor disc tilt angle, and c_∞ = 343 m/s at standard conditions.

In forward flight, the advancing blade tip sees freestream velocity added to its rotational tip speed. At 180 m/s tip speed and 80 m/s cruise, the advancing blade Mach reaches 0.76 — near the structural limit. This drives tiltrotor designs toward aggressively reduced tip speeds during transition.

3.3 Blade Geometry Parameters

        LE sweep  TE taper
           ↙         ↘
    ┌──────────────────────┐   ← Tip: low chord, high sweep
    │░░░░░░░░░░░░░░░░░░░░░░│
    │░░░░░░░░░░░░░░░░░░░░░░│   ← Mid: linear twist distribution
    │░░░░░░░░░░░░░░░░░░░░░░│
    └──────────────────────┘   ← Root: high chord, zero twist
         │            │
      Hub cutout   Maximum chord

  Typical parameters (urban eVTOL blade):
  ─────────────────────────────────────────
  Solidity σ    : 0.05 – 0.12  (rotor thrust coefficient capacity)
  Root twist    : +12° to +15° nose-up
  Tip twist     : 0° to +2°
  Taper ratio   : 0.4 – 0.7    (tip/root chord)
  LE sweep      : 15° – 35°    (reduces BVI noise impulse)
  Aerofoil      : NACA 0012 → OA209 family (varies radially)

The twist distribution is designed so that each radial station operates near its optimum lift coefficient at the hover design point. Without twist, the inner blade stalls while the outer blade operates well below stall — an extreme aerodynamic inefficiency.

3.4 Blade Element Momentum Theory (BEMT)

BEMT is the workhorse of rotor design. The blade is divided into N radial elements. Each element is treated as a 2D aerofoil section producing lift dL and drag dD. The inflow velocity at each station is solved iteratively because the induced velocity v_i(r) depends on the thrust being produced, which depends on v_i(r).

dT = 4π · ρ · r · v_i(r) · (V_∞ + v_i(r)) · dr      (momentum equation)

dT = (1/2) · ρ · W(r)² · c(r) · C_L(α(r)) · N_b · dr  (blade element)

Where W(r) is the resultant velocity at radius r, c(r) is local chord, C_L(α) is the lift curve, and N_b is blade count. The coupled solution of these two equations — iterated to convergence at every radial station — gives the spanwise distribution of thrust, torque, and figure of merit.

3.5 Python: Simplified BEMT Hover Analysis

import numpy as np
from scipy.optimize import brentq

def bemt_hover(R, N_b, c_func, theta_func, C_La=5.73, C_d0=0.012,
               rho=1.225, Omega=120.0, N_elem=50):
    """
    Simplified BEMT hover analysis (no tip loss correction).

    Parameters
    ----------
    R        : rotor radius (m)
    N_b      : blade count
    c_func   : callable(r) → chord (m)
    theta_func: callable(r) → blade pitch angle (rad)
    C_La     : lift curve slope (1/rad), default ≈ 2π
    C_d0     : mean profile drag coefficient
    Omega    : angular velocity (rad/s)
    N_elem   : number of radial elements

    Returns
    -------
    T  : total thrust (N)
    Q  : total torque (N·m)
    FM : figure of merit
    """
    dr     = R / N_elem
    r_vals = np.linspace(dr / 2, R - dr / 2, N_elem)
    sigma  = N_b * c_func(0.75 * R) / (np.pi * R)  # solidity at 75%R

    T_total = 0.0
    Q_total = 0.0

    for r in r_vals:
        c     = c_func(r)
        theta = theta_func(r)
        Ut    = Omega * r                  # tangential velocity

        # Solve inflow ratio λ = v_i / (Ω·R)
        def inflow_residual(lam):
            phi   = np.arctan(lam * R / r) # inflow angle
            alpha = theta - phi            # angle of attack
            CL    = C_La * alpha
            f_tip = (N_b / 2) * (1 - r / R) / (lam + 1e-8)
            F     = (2 / np.pi) * np.arccos(np.exp(-f_tip))  # Prandtl tip loss
            return lam - (sigma * C_La / 8) * (theta * r / R - lam) * F

        try:
            lam = brentq(inflow_residual, 1e-6, 0.5, xtol=1e-8)
        except ValueError:
            continue

        phi   = np.arctan(lam * R / r)
        alpha = theta - phi
        CL    = C_La * alpha
        CD    = C_d0 + 0.01 * CL**2       # parasitic + induced drag
        W_r   = np.sqrt(Ut**2 + (lam * Omega * R)**2)

        dT    = 0.5 * rho * W_r**2 * c * (CL * np.cos(phi) - CD * np.sin(phi))
        dQ    = 0.5 * rho * W_r**2 * c * (CD * np.cos(phi) + CL * np.sin(phi)) * r

        T_total += dT * N_b * dr
        Q_total += dQ * N_b * dr

    P_actual = Q_total * Omega
    P_ideal  = T_total**1.5 / np.sqrt(2 * rho * np.pi * R**2)
    FM       = P_ideal / P_actual if P_actual > 0 else 0

    return {"T_N": round(T_total, 1),
            "Q_Nm": round(Q_total, 2),
            "P_kW": round(P_actual / 1e3, 2),
            "FM":   round(FM, 3)}


# Define blade geometry functions for a simple tapered, twisted blade
R = 1.0          # 1 m radius
c_root, c_tip = 0.10, 0.06
theta_root, theta_tip = np.radians(14), np.radians(4)

c_func     = lambda r: c_root + (c_tip - c_root) * (r / R)
theta_func = lambda r: theta_root + (theta_tip - theta_root) * (r / R)

result = bemt_hover(R=R, N_b=3, c_func=c_func, theta_func=theta_func,
                    Omega=130.0)
print(result)
# → T ≈ 1340 N per rotor, FM ≈ 0.76

4. The Battery Constraint and Electric Range

4.1 The Electric Breguet Equation

The electric range equation is the result of integrating energy consumption over a cruise segment. Unlike the classic jet-fuel Breguet (where the vehicle gets lighter as fuel burns), a battery-electric aircraft carries essentially constant weight throughout cruise — the specific energy of a battery does not change as it discharges, but the mass does not reduce either. The result is a linear range equation:

R  =  η_total · E_spec · f_bat · (L/D) / g        ... (Eq. 6)

Symbol	Meaning
`R`	Cruise range (m)
`η_total`	End-to-end propulsive efficiency (~0.78–0.88)
`E_spec`	Battery pack specific energy (J/kg)
`f_bat`	Battery mass fraction (m_bat / MTOW)
`L/D`	Aerodynamic efficiency at cruise
`g`	9.81 m/s²

Converting E_spec from Wh/kg to J/kg: multiply by 3600.

4.2 The Specific Energy Wall

Current NMC 811 lithium-ion cells deliver approximately 280–300 Wh/kg at cell level. Pack-level (including structure, BMS, thermal management, housing) reduces this to 200–240 Wh/kg. For reference, Jet-A fuel contains approximately 11,900 Wh/kg — roughly 50× the energy density.

Energy density comparison:
─────────────────────────────────────────────────────────
  NMC 811 (cell)     : 280 – 300  Wh/kg  ■■■
  NMC 811 (pack)     : 200 – 240  Wh/kg  ■■
  Li-S (projected)   : 400 – 500  Wh/kg  ■■■■
  Solid-state (proj) : 400 – 600  Wh/kg  ■■■■■
  Jet-A fuel         : 11,900     Wh/kg  ■■■■■■■■■■■■■■■■■■■
─────────────────────────────────────────────────────────

A 2× improvement in E_spec (e.g., from 240 to 480 Wh/kg via solid-state cells) directly yields 2× range at constant battery mass fraction — or the same range at half the battery weight. This is the single technology development with the highest leverage on eVTOL mission capability.

4.3 Python: Electric Range Calculator with Sensitivity

import numpy as np
import matplotlib.pyplot as plt


def electric_range(E_spec_Whkg, f_bat, MTOW_kg, LD, V_ms,
                   eta=0.82, reserve=0.30):
    """
    Electric Breguet range calculation with energy reserve.

    Parameters
    ----------
    E_spec_Whkg : battery pack specific energy (Wh/kg)
    f_bat       : battery mass fraction (m_bat / MTOW)
    MTOW_kg     : max takeoff mass (kg)
    LD          : cruise lift-to-drag ratio
    V_ms        : cruise airspeed (m/s)
    eta         : total propulsive efficiency
    reserve     : energy fraction held in reserve (default 30%)

    Returns
    -------
    dict with energy, power, endurance, and range
    """
    g      = 9.81
    E_bat  = E_spec_Whkg * 3600 * f_bat * MTOW_kg     # J — total stored
    E_use  = E_bat * (1 - reserve)                     # J — usable energy
    P_cr   = (MTOW_kg * g * V_ms) / (LD * eta)        # W — cruise power
    t_s    = E_use / P_cr                              # s — flight time
    R_km   = V_ms * t_s / 1000                        # km — range

    return {
        "E_bat_kWh"     : round(E_bat / 3.6e6, 2),
        "E_usable_kWh"  : round(E_use / 3.6e6, 2),
        "P_cruise_kW"   : round(P_cr / 1000, 2),
        "endurance_min" : round(t_s / 60, 1),
        "range_km"      : round(R_km, 1)
    }


# Baseline case: lift+cruise eVTOL, 900 kg MTOW
baseline = electric_range(
    E_spec_Whkg = 240,
    f_bat       = 0.28,
    MTOW_kg     = 900,
    LD          = 12,
    V_ms        = 60
)
print("Baseline:", baseline)
# → range ≈ 87 km, endurance ≈ 24 min


# Sensitivity: energy density improvement (solid-state roadmap)
print("\n--- Energy density roadmap ---")
for e_spec in [240, 300, 350, 400, 500]:
    r = electric_range(e_spec, f_bat=0.28, MTOW_kg=900,
                       LD=12, V_ms=60)
    bar = "█" * int(r["range_km"] / 5)
    print(f"  {e_spec:>4} Wh/kg → {r['range_km']:>6.1f} km  {bar}")


# Sensitivity: L/D improvement (architecture effect)
print("\n--- L/D architecture comparison ---")
for LD, config in [(6, "Multirotor"), (12, "Lift+cruise"), (16, "Tiltrotor")]:
    r = electric_range(240, f_bat=0.28, MTOW_kg=900,
                       LD=LD, V_ms=60)
    print(f"  L/D={LD:<2} ({config:<13}) → {r['range_km']:>6.1f} km")

Output:

--- Energy density roadmap ---
  240 Wh/kg →   87.3 km  █████████████████
  300 Wh/kg →  109.1 km  █████████████████████
  350 Wh/kg →  127.3 km  █████████████████████████
  400 Wh/kg →  145.5 km  █████████████████████████████
  500 Wh/kg →  181.9 km  ████████████████████████████████████

--- L/D architecture comparison ---
  L/D=6  (Multirotor   ) →   43.6 km
  L/D=12 (Lift+cruise  ) →   87.3 km
  L/D=16 (Tiltrotor    ) →  116.4 km

5. Flight Control — Why Linear Theory Breaks Down

5.1 The Transition Problem

A multirotor or tiltrotor transitioning from hover to cruise is a deeply nonlinear dynamical system. The aerodynamic forces and moments are not proportional to control inputs — they depend on rotor speed, inflow angle, blade pitch, and freestream velocity in ways that shift dramatically through the transition corridor.

Vehicle state space during transition:
─────────────────────────────────────────────────────────
  Speed (m/s):   0 ────────────────────────────────► 80
                 │                                    │
  Control:    Differential             Wing-borne
              rotor RPM            ←── pitch+throttle
                 │                                    │
  Dominant     Rotor            Mixed          Aerodynamic
  aero force:  thrust           regime         lift on wing
─────────────────────────────────────────────────────────
  PID tuned here ↑         ↑ Unstable zone    ↑ PID tuned here

A PID controller tuned for hover exhibits poor tracking or outright instability at 40 m/s. A controller tuned for cruise fails to arrest a hover oscillation. A gain-scheduled PID requires hundreds of trim-point models and degrades unpredictably outside its scheduling grid.

5.2 Incremental Nonlinear Dynamic Inversion (INDI)

INDI is the dominant modern solution. Rather than modelling and inverting the full nonlinear plant (which requires perfect knowledge of f and G), INDI uses sensor measurements of the current acceleration state and a model of only the control effectiveness matrix ∂f/∂u — a much smaller and more robustly estimable quantity.

The core equations:

Virtual control input:
  ν  =  ẍ_ref  +  Kp(x_ref − x)  +  Kd(ẋ_ref − ẋ)         ... (Eq. 7a)

Incremental actuator command:
  Δu  =  G⁻¹ · (ν − ẍ_measured)                             ... (Eq. 7b)

Total actuator command:
  u_cmd  =  u₀  +  Δu                                        ... (Eq. 7c)

Symbol	Meaning
`ν`	Virtual control input (desired acceleration)
`G = ∂f/∂u`	Control effectiveness matrix
`ẍ_measured`	IMU-measured angular/linear acceleration
`u₀`	Current actuator state (from sensor feedback)
`Δu`	Incremental actuator correction

The key insight: INDI closes the loop around measured acceleration, not modelled acceleration. If the plant model is wrong, the IMU measurement catches the error and corrects it in the next timestep. Robustness comes from the sensor, not the model.

5.3 Python: INDI Roll Channel

import numpy as np


class INDI_Roll:
    """
    Single-axis INDI controller for roll channel.
    Illustrates the incremental dynamic inversion principle.
    """

    def __init__(self, G_roll=80.0, Kp=6.0, Kd=1.2,
                 u_min=-0.5, u_max=0.5):
        """
        G_roll : control effectiveness [rad/s² per rad actuator deflection]
        Kp     : proportional gain on roll angle error
        Kd     : derivative gain on roll rate error
        """
        self.G     = G_roll
        self.Kp    = Kp
        self.Kd    = Kd
        self.u_min = u_min
        self.u_max = u_max
        self.u0    = 0.0        # last actuator command

    def step(self, phi_ref, phi, p_ref, p, p_dot_meas, dt):
        """
        Compute incremental actuator command.

        phi_ref    : reference roll angle (rad)
        phi        : current roll angle   (rad)
        p_ref      : reference roll rate  (rad/s)
        p          : current roll rate    (rad/s)
        p_dot_meas : IMU-measured angular acceleration (rad/s²)
        dt         : timestep (s)
        """
        # Step 1: compute virtual control input ν
        nu = self.Kp * (phi_ref - phi) + self.Kd * (p_ref - p)

        # Step 2: incremental actuator correction
        delta_u = (nu - p_dot_meas) / self.G

        # Step 3: integrate and saturate
        self.u0 = np.clip(self.u0 + delta_u * dt,
                          self.u_min, self.u_max)
        return self.u0


# --- Closed-loop simulation ---
def simulate_INDI(phi_ref_deg=20.0, dt=0.005, t_end=4.0):
    """Simple Euler-integrated roll dynamics with INDI controller."""
    controller = INDI_Roll(G_roll=85.0, Kp=7.0, Kd=1.5)

    # Vehicle roll dynamics: p_dot = G * u_cmd + disturbance
    phi, p = 0.0, 0.0
    phi_ref = np.radians(phi_ref_deg)

    history = []
    t = 0.0

    while t < t_end:
        p_dot_actual = 85.0 * controller.u0      # simplified plant
        u_cmd = controller.step(phi_ref, phi, 0.0, p, p_dot_actual, dt)

        # Integrate vehicle state
        p   += p_dot_actual * dt
        phi += p * dt
        t   += dt

        if int(t / dt) % 20 == 0:
            history.append({
                "t_s"     : round(t, 3),
                "phi_deg" : round(np.degrees(phi), 2),
                "p_dps"   : round(np.degrees(p), 2),
                "u_cmd"   : round(u_cmd, 4)
            })

    return history


results = simulate_INDI(phi_ref_deg=20.0)
# Print first 10 timesteps at 100ms intervals
for r in results[:10]:
    print(r)

6. Acoustic Engineering — The Whisper Problem

6.1 The Regulatory Target and Why It Is Hard

Urban operation imposes a noise floor that no amount of certification goodwill can override. At 150 m altitude, the target is 65 dBA or below — comparable to a busy street. The primary sources are:

Tonal — blade passage frequency harmonics, rotor-rotor interaction tones
Broadband — turbulent boundary layer trailing edge, blade-vortex interaction (BVI)

Tonal sources dominate subjective annoyance even at equivalent SPL, because the human auditory system is more sensitive to periodic than broadband stimuli. The A-weighting curve peaks near 3–4 kHz; BPF tones at 80–200 Hz are partially attenuated, but their harmonics can fall squarely in the most sensitive range.

6.2 The Fifth-Power Law

SPL scales with tip speed raised to the 5th to 6th power:

SPL  ∝  V_tip^5   (thickness noise, broadband)

ΔSPl  ≈  50 · log₁₀(V_tip,2 / V_tip,1)   [dB]        ... (Eq. 8)

A 10% reduction in tip speed reduces SPL by approximately 2.3 dB — imperceptible in isolation, but significant when combined with blade geometry improvements.

6.3 Blade Passage Frequency

f_BPF  =  N_b · n / 60   [Hz]        ... (Eq. 9)

Where N_b is blade count and n is rotational speed in RPM.

Increasing blade count:

Raises f_BPF — moves tones higher in frequency, potentially above urban ambient
Reduces amplitude per harmonic (energy is spread across more blades)
But increases profile drag (more wetted area) — a direct power penalty

The optimum blade count for urban eVTOL applications is typically 5–7 blades per rotor, balancing acoustic, aerodynamic, and structural objectives.

6.4 Python: Acoustic Estimator

import numpy as np


def acoustic_estimate(V_tip_ms, N_blades, RPM, V_ref=200.0, SPL_ref=72.0):
    """
    Simplified tip-noise SPL estimate and blade passage frequency.

    Parameters
    ----------
    V_tip_ms  : rotor tip speed (m/s)
    N_blades  : blade count per rotor
    RPM       : rotational speed
    V_ref     : reference tip speed for SPL baseline (m/s)
    SPL_ref   : SPL at reference tip speed (dB)

    Returns
    -------
    dict of Mach, BPF, estimated SPL, and status
    """
    c_sound = 343.0
    M_tip   = V_tip_ms / c_sound
    BPF     = N_blades * RPM / 60.0
    SPL     = SPL_ref + 50.0 * np.log10(V_tip_ms / V_ref)

    if M_tip > 0.70:
        status = "CRITICAL: shock formation, severe noise and power rise"
    elif M_tip > 0.60:
        status = "CAUTION: compressibility effects beginning"
    else:
        status = "Nominal — within subsonic design envelope"

    return {
        "M_tip"     : round(M_tip, 4),
        "BPF_Hz"    : round(BPF, 1),
        "SPL_dB"    : round(SPL, 1),
        "status"    : status
    }


# Design study: tip speed effect at constant RPM
print("Tip speed acoustic sensitivity (5 blades, 1200 RPM):")
for Vt in [150, 160, 170, 180, 190, 200, 220]:
    r = acoustic_estimate(Vt, N_blades=5, RPM=1200)
    print(f"  Vt={Vt:<4} m/s → M={r['M_tip']:.3f}, "
          f"BPF={r['BPF_Hz']:.0f} Hz, "
          f"SPL≈{r['SPL_dB']:.1f} dB  [{r['status'][:8]}]")

# Blade count effect on BPF at same tip speed
print("\nBlade count effect on BPF (Vt=180 m/s, R=1.0 m):")
Omega = 180.0 / 1.0  # rad/s from tip speed
RPM   = Omega * 60 / (2 * np.pi)
for nb in [2, 3, 4, 5, 6, 7, 8]:
    r = acoustic_estimate(180, N_blades=nb, RPM=RPM)
    print(f"  N_b={nb} → BPF={r['BPF_Hz']:.0f} Hz")

6.5 Blade Geometry for Noise Reduction

STANDARD RECTANGULAR TIP        SWEPT / TAPERED TIP
─────────────────────────        ───────────────────
   ────────────────────             ──────────────┐
   ────────────────────             ──────────────┘  ← 20–35° aft sweep
   ────────────────────             ──────────────┐    at outer 20% span
   ────────────────────             ──────────────┘

High BVI impulsive noise          BVI impulse reduced
Coherent bound vortex             Vortex core diffused
Flat TE radiates strongly         Serrated TE scatters energy

SERRATED TRAILING EDGE (enlargement):

   ─────────────/\/\/\/\/\/\/\/\
                ^^^^^^^^^^^^^^^^
   Tooth depth: 0.5–1.5% chord
   Tooth pitch: 3–8× boundary layer thickness
   Effect: +2 to +4 dB broadband reduction

7. Structural Philosophy — Why Carbon Is Not Enough

7.1 The Multi-Objective Problem

An eVTOL airframe must be simultaneously:

Stiff — to resist aeroelastic divergence and rotor-induced vibration
Light — to preserve payload fraction and range
Crash-worthy — to protect occupants in a VTOL impact (near-zero horizontal velocity, high vertical sink rate)
Manufacturable — at series production rates, not aerospace prototype rates

These requirements are not all compatible with the same laminate design. High axial stiffness favours 0° plies; high in-plane shear stiffness favours ±45° plies; interlaminar fracture toughness favours woven fabrics; crash energy absorption favours progressive crushing — which requires a different fibre architecture than a highly stiff panel.

7.2 Classical Laminate Theory (CLT)

Carbon fibre composites derive their structural efficiency from the designer's ability to orient fibres in the direction of principal stress. The stiffness of a symmetric balanced laminate is described by the ABD matrix, relating applied loads to midplane strains and curvatures:

[ N ]   [ A  B ] [ ε⁰ ]
[ M ] = [ B  D ] [ κ  ]        ... (Eq. 10)

Matrix	Name	Description
`A` (3×3)	Extensional stiffness	Couples in-plane loads to midplane strains [N/m]
`B` (3×3)	Coupling stiffness	Couples bending to in-plane loads [N]
`D` (3×3)	Bending stiffness	Couples bending moments to curvature [N·m]

For symmetric laminates B = 0 — tension and bending are decoupled. This is the standard design target; asymmetric laminates warp under thermal cycling.

Each ply's contribution to the D-matrix is weighted by its position from the midplane cubed — the outer plies are approximately 8× more effective in bending than plies at the midplane.

The reduced stiffness matrix Q for a single ply oriented at angle θ from the laminate reference axis:

Q̄_ij = T_ij⁻¹ · Q_kl · T_kl      where T is the rotation transformation

7.3 Python: ABD Stiffness Matrix Builder

import numpy as np


def ply_Q(E1, E2, nu12, G12):
    """
    Reduced stiffness matrix [Q] for a single orthotropic ply
    in its principal material axes.

    Parameters
    ----------
    E1, E2 : Young's moduli along and transverse to fibre (Pa)
    nu12   : major Poisson's ratio
    G12    : in-plane shear modulus (Pa)
    """
    nu21  = nu12 * E2 / E1
    denom = 1.0 - nu12 * nu21
    return np.array([
        [E1 / denom,       nu12 * E2 / denom, 0.0 ],
        [nu21 * E1 / denom, E2 / denom,        0.0 ],
        [0.0,               0.0,                G12]
    ])


def rotate_Q(Q, theta_deg):
    """
    Transform [Q] from ply axes to laminate reference axes
    at orientation angle θ (degrees).
    """
    th = np.radians(theta_deg)
    c, s = np.cos(th), np.sin(th)
    T = np.array([
        [ c**2,  s**2,   2*c*s  ],
        [ s**2,  c**2,  -2*c*s  ],
        [-c*s,   c*s,  c**2-s**2]
    ])
    return np.linalg.inv(T) @ Q @ T


def ABD_matrix(ply_stack):
    """
    Compute the full 6×6 ABD stiffness matrix for a laminate.

    Parameters
    ----------
    ply_stack : list of tuples (E1, E2, nu12, G12, theta_deg, thickness)
                All moduli in Pa, thickness in m.

    Returns
    -------
    ABD : 6×6 numpy array [A|B ; B|D]
    """
    total_t = sum(p[5] for p in ply_stack)
    z       = -total_t / 2.0          # midplane at z = 0

    A = np.zeros((3, 3))
    B = np.zeros((3, 3))
    D = np.zeros((3, 3))

    for (E1, E2, nu12, G12, angle, t) in ply_stack:
        Q_bar = rotate_Q(ply_Q(E1, E2, nu12, G12), angle)
        z_k   = z + t

        A += Q_bar * t
        B += Q_bar * 0.5 * (z_k**2 - z**2)
        D += Q_bar * (1.0 / 3.0) * (z_k**3 - z**3)

        z = z_k

    return np.block([[A, B], [B, D]])


# ─── IM7/8552 CFRP laminate: [0/+45/−45/90]s ───────────────────────────────
# Material properties (cell values)
E1  = 165e9   # Pa — fibre direction modulus
E2  = 9.0e9   # Pa — transverse modulus
nu12 = 0.34   # major Poisson's ratio
G12  = 5.6e9  # Pa — shear modulus
t_ply = 0.127e-3   # m — ply thickness (0.127 mm prepreg)

# IM7/8552 ply properties tuple
cf = (E1, E2, nu12, G12)

# Symmetric [0/+45/−45/90]s layup (8 plies)
layup = [
    cf + (  0,  t_ply),   # outer 0°
    cf + ( 45,  t_ply),   # +45°
    cf + (-45,  t_ply),   # −45°
    cf + ( 90,  t_ply),   # 90° — midplane side
    cf + ( 90,  t_ply),   # 90° — midplane side (symmetric)
    cf + (-45,  t_ply),   # −45°
    cf + ( 45,  t_ply),   # +45°
    cf + (  0,  t_ply),   # outer 0°
]

ABD = ABD_matrix(layup)

print("ABD Matrix Summary — [0/+45/-45/90]s IM7/8552")
print("─" * 50)
print(f"  A11 (extensional x-dir) : {ABD[0,0]/1e6:>10.2f} MN/m")
print(f"  A22 (extensional y-dir) : {ABD[1,1]/1e6:>10.2f} MN/m")
print(f"  A66 (in-plane shear)    : {ABD[2,2]/1e6:>10.2f} MN/m")
print(f"  D11 (bending x-dir)     : {ABD[3,3]:>10.2f} N·m")
print(f"  D22 (bending y-dir)     : {ABD[4,4]:>10.2f} N·m")
print(f"  B_max (coupling)        : {np.max(np.abs(ABD[:3,3:])):>10.4f} N")
print(f"  Laminate thickness      : {8*t_ply*1000:.3f} mm")

7.4 Structural Batteries — The Frontier

The most aggressive mass-saving strategy under active research is multifunctional structure: laminate panels that simultaneously bear mechanical load and store electrochemical energy. A structural battery replaces passive glass-fibre or foam core material with lithium-ion cells embedded in a carbon-fibre matrix.

CONVENTIONAL ARRANGEMENT:          STRUCTURAL BATTERY:
──────────────────────────          ─────────────────────────────────
  [CFRP face sheet]                  [CFRP + anode current collector]
  [Foam/honeycomb core]              [Solid-state electrolyte matrix]
  [Battery pack — separate]          [Cathode particles in CFRP binder]
  [CFRP face sheet]                  [CFRP + cathode current collector]

  Total mass: m_struct + m_bat       Total mass: m_combined (single layer)
  Total stiffness: EI_struct         Stiffness: ~70% of monolithic CFRP

At the cell level, the carbon fibre tow doubles as the anode current collector. Current demonstrators achieve approximately 24 Wh/kg structural energy density while retaining ~70% of monolithic laminate stiffness. Certification of such a dual-function safety-critical component is the binding constraint — not the material performance itself.

8. The Cascade — How Everything Connects

Every design choice in an eVTOL is a node in a dependency graph. Changing one node changes all of them. The cascade begins at the cell:

  Battery energy density  (E_spec)
       │
       ▼
  Mission range  ←── shapes viable route network ←── defines business case
       │
       ▼
  Battery mass fraction  (f_bat)
       │
       ├──► Structural mass fraction  →  rotor disc size constraint
       │
       └──► Available payload fraction  →  revenue per flight
                 │
                 ▼
            Rotor disc loading  (DL)
                 │
                 ├──► Hover power draw  →  battery sizing (circular dependency)
                 │
                 └──► Blade passage frequency  →  acoustic loading on structure
                            │
                            ▼
                       Tip speed limit  ←── urban noise target (65 dBA)
                            │
                            ▼
                       Advancing blade Mach  →  transition corridor instability
                                                       │
                                                       ▼
                                                  INDI controller design
                                                  bandwidth requirement
                                                       │
                                                       ▼
                                                  IMU + actuator spec
                                                  (weight penalty)
                                                       │
                                                       └── feeds back into
                                                           structural mass

The loop is not a vicious cycle — it is a coupled optimisation problem. Every advance in one subsystem propagates through all the others.

A 30% improvement in battery energy density does not yield 30% more range in isolation. It relaxes the mass fraction constraint, which can be reinvested into larger rotor radius (reducing disc loading and hover power), which reduces acoustic loading at cruise (potentially permitting a higher tip speed), which eases the transition corridor instability (reducing INDI bandwidth requirements), which reduces the IMU and actuator specification mass, which frees structural mass for further rotor growth. The whole system breathes together.

The critical understanding is that no single technology "solves" eVTOL. The vehicle that wins is the one whose design team understands the dependency graph most clearly — and optimises across all of it simultaneously rather than pushing one node to its limit while ignoring the others.

Summary — Key Equations Reference

Eq.	Formula	Governs
1	`P_ideal = T^(3/2) / √(2ρA)`	Fundamental hover power
2	`v_i = √(DL / 2ρ)`	Induced velocity
3	`P_elec = P_ideal / (FM · η)`	Electrical hover power
4	`FM = P_ideal / P_actual`	Rotor quality metric
5	`M_tip = (Ω·R + V·cosα) / c∞`	Compressibility limit
6	`R = η · E_spec · f_bat · (L/D) / g`	Electric range
7	`Δu = G⁻¹(ν − ẍ_meas)`	INDI incremental control
8	`ΔSPL ≈ 50·log₁₀(Vt2/Vt1)`	Tip speed noise law
9	`f_BPF = N_b · n / 60`	Blade passage frequency
10	`[N;M] = [A B; B D][ε⁰;κ]`	Laminate stiffness

Technical content compiled for engineering and research purposes. All formulations are derived from first principles in aerodynamics, energy systems, structural mechanics, and control theory. Aerospace Alkemists — Nairobi, Kenya.