Python Lab Reference

Every library, method and pattern you need for the Aircraft Engines labs.

49 entries across 9 libraries

🔢NumPy

The arithmetic of every lab: square roots for sound speeds and jet velocities, powers for isentropic ratios, trigonometry for velocity triangles. Angles are in radians.

import numpy as np

Maths and angles

np.sqrt / **

Roots and powers: speed of sound, jet velocity, isentropic ratios.

a = np.sqrt(1.4 * 287 * T)          # speed of sound
Tt_T = 1 + 0.2 * M**2
pt_p = Tt_T ** 3.5

Write exponents like (γ − 1)/γ as floats (0.4/1.4 or 2/7), never integer division.

np.radians / np.degrees

Convert between degrees (how angles are quoted) and radians (what trig functions take).

beta = np.radians(35)
Mn1 = 2.2 * np.sin(beta)
print(np.degrees(np.arctan(0.5)))

np.arctan2

np.arctan2(y, x)

The angle of a vector in the right quadrant — flow angles in velocity triangles.

beta1 = np.degrees(np.arctan2(U - c_theta, c_x))

np.hypot

np.hypot(a, b)

Length of a vector from its components: absolute and relative velocities.

W1 = np.hypot(c_x, U)

np.exp / np.log

Arrhenius rates, Breguet logarithms, entropy changes.

k = A * np.exp(-Ta / T)
R = u0 / (g * tsfc) * LD * np.log(mi / mf)

Arrays, sweeps and tables

np.linspace / np.arange

Sweep a pressure ratio, Mach number or speed.

pi_c = np.linspace(2, 40, 77)
machs = np.arange(0.0, 0.95, 0.05)

arange excludes its end point; linspace includes it.

np.argmax / np.max

Find the optimum of a sweep — the pressure ratio for maximum specific work.

i = np.argmax(work)
print(pi_c[i], work[i])

np.interp

np.interp(x, xp, fp)

Linear interpolation in a table — a schedule, a map, an atmosphere.

tt4_limit = np.interp(N, grid, table)

np.trapezoid

np.trapezoid(y, x)

Integrate sampled data — a mean hoop stress, a disk mass.

mass = np.trapezoid(rho * 2 * np.pi * r * h, r)

np.trapz was renamed np.trapezoid in NumPy 2.

np.linalg.solve / np.linalg.matrix_rank

Solve a linear system — bearing reactions, mount link forces — and check it is determinate.

A = np.array([[1.0, 1.0], [x1, x2]])
R = np.linalg.solve(A, [F_total, M_total])
print(np.linalg.matrix_rank(A6))   # 6 for a determinate mount

np.cross

Moment of a force about a point: r × F.

M = np.cross(p, d)   # column of the mount matrix

np.polyfit / np.roots / np.corrcoef

Fit a trend to engine data, solve a polynomial, measure a correlation.

slope, icpt = np.polyfit(year, bpr, 1)
r = np.corrcoef(bpr, tsfc)[0, 1]

🎯scipy.optimize

Most engine relations are implicit — Mach from an area ratio, Tt4 for a target thrust, a shock angle for a deflection. Root finders turn them around; optimisers find the best design.

from scipy.optimize import brentq, fsolve, minimize_scalar, minimize

Roots

brentq→ float

brentq(f, a, b, xtol=1e-12)

The root of f between a and b, where f(a) and f(b) have opposite signs. Fast and certain.

M_sup = brentq(lambda M: area_ratio(M) - 2.0, 1.0001, 5.0)
M_sub = brentq(lambda M: area_ratio(M) - 2.0, 0.01, 0.9999)

Choose the bracket to pick the branch: subsonic and supersonic Mach numbers share every area ratio.

fsolve→ ndarray

fsolve(F, x0, args=(...), full_output=True)

Several nonlinear equations in several unknowns — the off-design turbofan (Lab 8-2a).

x, info, ok, msg = fsolve(residuals, guess, args=(M0, T0, Tt4, ref), full_output=True)
if ok != 1:
    raise RuntimeError(msg)

It needs a good first guess. Step the throttle in small increments and start each solve from the last one (continuation).

Optimisation

minimize_scalar

minimize_scalar(f, bounds=(a, b), method='bounded')

The best value of one variable: optimum fan pressure ratio, insulation thickness.

res = minimize_scalar(lambda t: foam(t) + boiloff(t), bounds=(0.005, 0.5), method='bounded')
print(res.x, res.fun)

minimize

minimize(f, x0, bounds=..., method=...)

Several design variables at once — the angles of a two-ramp inlet.

res = minimize(lambda x: -recovery(*x), [8, 8], bounds=[(1, 20), (1, 20)])

minimize finds minima: maximise by minimising the negative.

∫scipy.integrate, interpolate, linalg

Time histories (a spool accelerating), smooth curves through rig data (a compressor map) and the natural frequencies of a rotor.

from scipy.integrate import solve_ivp, quad
from scipy.interpolate import CubicSpline
from scipy.linalg import eigh

Differential equations

solve_ivp→ OdeResult (.t, .y, .t_events)

solve_ivp(rhs, (t0, t1), y0, events=..., dense_output=True, max_step=...)

Integrate dy/dt = rhs(t, y) — spool speed (Lesson 8-3), a rotating disk outward from its bore (Lesson 9-2).

sol = solve_ivp(lambda t, y: [dNdt(y[0])], (0, 30), [0.6], max_step=0.05)
N_end = sol.y[0][-1]

events

Stop the integration when something happens: 95 % thrust reached, a surge margin crossing zero.

def hit(t, y):
    return thrust(y[0]) - 0.95 * F_rated
hit.terminal = True
sol = solve_ivp(rhs, (0, 20), [0.6], events=hit)
print(sol.t_events[0][0])

quad

quad(f, a, b)

A definite integral of a function.

I, err = quad(lambda x: x**2, 0, 1)

Interpolation, eigenvalues, rotations

CubicSpline

CubicSpline(x, y)

A smooth curve through measured points — a compressor speed line.

line = CubicSpline(m_points, pr_points)
pr = line(0.95)

eigh

eigh(K, M)

The generalised eigenproblem K φ = ω² M φ for symmetric matrices: a rotor’s natural frequencies and mode shapes.

w2, modes = eigh(K, M)
rpm = np.sqrt(w2[:3]) * 30 / np.pi

Rotation

from scipy.spatial.transform import Rotation as R

Compose rotations — a vectoring nozzle’s pitch and yaw.

r = R.from_euler('yz', [10, 10], degrees=True)
thrust_dir = r.apply([1, 0, 0])

∂SymPy

Let the computer do the algebra: balance a reaction, derive the running line, check a derivation from the theory.

import sympy as sp

Symbols and equations

sp.symbols

sp.symbols('a b c', positive=True)

Create symbols. Declaring them positive lets SymPy simplify square roots.

pi, dc, dt = sp.symbols('pi dc dt', positive=True)

Use positive helper symbols (dc = τc − 1) rather than expressions that may look negative.

sp.Eq / sp.solve

Write an equation and solve it.

a = sp.symbols('a')
sp.solve(sp.Eq(2 * a, 2 * 12 + sp.Rational(23, 2)), a)   # [71/4]

sp.simplify / sp.factor / sp.expand

Tidy an expression.

running = sp.simplify(C * pi / sp.sqrt(T42_work))

sp.Rational

Exact fractions — no rounding in a stoichiometric balance.

sp.Rational(23, 2)   # 23/2 exactly

sp.lambdify

sp.lambdify(x, expr, 'numpy')

Turn a derived expression into a fast numerical function.

m_run = sp.lambdify(pi, m_d * ratio, 'numpy')
print(m_run(10.2))

sp.diff

Differentiate — find where a specific work or thrust is greatest.

tau = sp.symbols('tau', positive=True)
sp.solve(sp.diff(work(tau), tau), tau)

🐼pandas

Engine data, station tables, LTO cycles and missions: tables with named rows and columns.

import pandas as pd

Tables

pd.DataFrame

Build a table from columns.

df = pd.DataFrame({'mode': ['take-off', 'climb-out'], 'thrust %': [100, 85], 'minutes': [0.7, 2.2]})

.set_index / .loc

Name the rows and pick one out.

df = df.set_index('mode')
print(df.loc['take-off', 'minutes'])

column arithmetic

Whole columns at once — fuel per mode, emissions per mode.

df['fuel (kg)'] = df['fuel flow'] * df['minutes'] * 60
total = df['fuel (kg)'].sum()

pd.read_csv(io.StringIO(...))

Read a small dataset written into the lab.

import io
data = '''engine,year,bpr\nJT9D-7R4,1982,4.8\nGE90-115B,2004,9.0'''
df = pd.read_csv(io.StringIO(data))

.round / .to_string

Print a table neatly.

print(df.round({'fuel (kg)': 1}).to_string())

📏pint

Quantities that carry their units: mixing lbf and newtons, or kilowatts and joules, becomes an error instead of a silent mistake.

import pint
ureg = pint.UnitRegistry()

Quantities

ureg.Quantity / multiplication

Attach a unit to a number.

F = 105000 * ureg.lbf
m_dot = 0.685 * ureg.kg / ureg.s

.to()

Convert, checking the dimensions agree.

print(F.to('kN'))                 # 467.1 kilonewton
print((m_dot * cp * dT).to('kW'))

.magnitude / f"{q:~P}"

The bare number, or a compact printout.

print(f"{F.to('kN'):.1f~P}")   # 467.1 kN

🌡️CoolProp

Real-gas properties of air and other fluids: specific heats that change with temperature, enthalpies and entropies from reference equations of state.

import CoolProp.CoolProp as CP

Properties

CP.PropsSI

CP.PropsSI(output, 'T', T, 'P', p, 'Air')

One property from two others, in SI units.

cp = CP.PropsSI('C', 'T', 1200, 'P', 101325, 'Air')
h = CP.PropsSI('H', 'T', 800, 'P', 3e6, 'Air')

Outputs: H (J/kg), S (J/(kg·K)), C (cp), CVMASS, D (kg/m³), T, P.

AbstractState

AS = CP.AbstractState('HEOS', 'Air')

The fast way inside loops: update the state, then read properties.

AS.update(CP.PT_INPUTS, 3e6, 800)
h = AS.hmass(); s = AS.smass()

Calling PropsSI thousands of times in a loop is slow; AbstractState is several times faster.

🔥Cantera

Chemical equilibrium and kinetics: flame temperatures with dissociation, ignition delays, thermal NO — the chemistry of Units 4 and 7.

import cantera as ct

Equilibrium

a Jet A phase

Build an ideal gas from the NASA species database, including the Jet-A(g) surrogate.

names = ['N2', 'O2', 'Ar', 'CO2', 'H2O', 'CO', 'H2', 'OH', 'O', 'H', 'NO', 'Jet-A(g)']
species = [s for s in ct.Species.list_from_file('nasa_gas.yaml') if s.name in names]
gas = ct.Solution(thermo='ideal-gas', species=species)

set_equivalence_ratio / TP

Set the mixture and its state.

AIR = 'O2:0.2095, N2:0.7809, Ar:0.0093, CO2:0.0004'
gas.set_equivalence_ratio(1.0, 'Jet-A(g):1', AIR)
gas.TP = 850, 40e5

equilibrate('HP')

Burn at constant enthalpy and pressure: the adiabatic flame temperature.

gas.equilibrate('HP')
print(gas.T)

Kinetics

ct.Solution('gri30.yaml')

Detailed natural-gas chemistry (GRI-Mech 3.0) — the course’s stand-in for kerosene kinetics.

gas = ct.Solution('gri30.yaml')
gas.set_equivalence_ratio(1.0, 'CH4:1', 'O2:1, N2:3.76')

IdealGasConstPressureReactor / ReactorNet

Let a mixture react in time: ignition delay, NO formation.

r = ct.IdealGasConstPressureReactor(gas, clone=False)
net = ct.ReactorNet([r])
net.advance(2e-3)
print(r.thermo['NO'].X)

Pass clone=False so the reactor works on your gas object.

✈️Engine patterns

The small pieces of code the labs build again and again.

Gas dynamics

isentropic ratios

Stagnation-to-static ratios for any Mach number.

def t_ratio(M, g=1.4):
    return 1 + (g - 1) / 2 * M**2
def p_ratio(M, g=1.4):
    return t_ratio(M, g) ** (g / (g - 1))

area–Mach inversion

Mach number from A/A*, on the branch you choose.

def area_ratio(M, g=1.4):
    return (1 / M) * ((2 / (g + 1)) * t_ratio(M, g)) ** ((g + 1) / (2 * (g - 1)))
M = brentq(lambda M: area_ratio(M) - AR, 1.0001, 6)

normal shock

Downstream Mach number and stagnation-pressure ratio.

M2 = np.sqrt((1 + 0.2 * M1**2) / (1.4 * M1**2 - 0.2))

Cycles and matching

station table

March through the engine station by station, then print the table.

st = {}
st['2'] = (Tt2, pt2)
st['3'] = (Tt2 * tau_c, pt2 * pi_c)
for k, (T, p) in st.items():
    print(f'{k:>3}  {T:7.1f} K  {p / 1e3:8.1f} kPa')

continuation

Solve a nonlinear system along a path, each answer the next guess.

g = None
for Tt4 in np.arange(1750, 949, -25):
    p = off_design(M0, T0, p0, Tt4, ref, g)
    g = p["x"]

shooting

Turn a boundary-value problem into a root: guess the unknown start, integrate, correct.

u0 = brentq(lambda u0: shoot(u0).y[1, -1] - sigma_rim, 0.0, 1e-2)