Aeronautical engineer: MQL recipes

Recipes

Standard atmosphere at FL350
Cruise lift, drag and minimum drag speed
Best glide of a sailplane from its measured polar
Parabolic drag polar from wind tunnel data
Take-off ground roll of a twin-jet
Wind triangle of a flight plan leg and runway wind components
Load factor, radius and rate of a level turn
Breguet range, endurance and trip fuel of an airliner
Mach, true and calibrated airspeed from the air data computer
Lift of a Joukowski airfoil with complex numbers

executed on MentDB

Standard atmosphere at FL350

At 35000 ft, compute the ISA temperature, pressure, density and speed of sound, then the Mach number of an airliner flying at 480 kn true airspeed.

Pipeline

01 · convertFeet to metresunit convert→02 · isaT, p, rho, amath calc→03 · machTAS in m/s / aunit convert

Script

admin

#1. Flight level 350 in metres;
-> "[h]" (unit convert 35000 "ft" "m");
#2. ISA troposphere: temperature, pressure, density, speed of sound;
-> "[T]" (math calc "288.15-0.0065*[h]");
-> "[p]" (math calc "101325*([T]/288.15)^5.25588");
-> "[rho]" (math calc "[p]/(287.053*[T])");
-> "[a]" (math calc "sqrt(1.4*287.053*[T])");
#3. Mach number for a true airspeed of 480 kn;
-> "[M]" (math calc (concat (unit convert 480 "kn" "m/s") "/[a]"));
concat "h = " (math round [h] 1) " m, T = " (math round (unit convert [T] "K" "degC") 2) " degC, p = " (math round (unit convert [p] "Pa" "hPa") 1) " hPa, rho = " (math round [rho] 4) " kg/m3, a = " (math round (unit convert [a] "m/s" "kn") 1) " kn, Mach " (math round [M] 3)

Result

mentdb

h = 10668 m, T = -54.34 degC, p = 238.4 hPa, rho = 0.3796 kg/m3, a = 576.4 kn, Mach 0.833

Go further

Add the stratosphere above 11 km with an exponential pressure law
Get the calibrated airspeed from the Mach number

executed on MentDB

Cruise lift, drag and minimum drag speed

For an A320 class airliner at FL350, compute the lift coefficient, the lift-to-drag ratio and the thrust required, then find the speed of minimum drag on the parabolic polar.

Pipeline

01 · liftCL for L = Wmath calc→02 · polarDrag as a function of Vconcat→03 · optimiseMinimum drag speednum minimize→04 · convertm/s to knotsunit convert

Script

admin

#1. A320 class aircraft at FL350: mass, wing area, ISA density, Mach 0.78 with a = 295.07 m/s;
-> "[W]" (math calc "64000*_g0");
-> "[S]" 122.6;
-> "[rho]" 0.3796;
-> "[V]" (math calc "0.78*295.07");
#2. Lift coefficient for L = W, drag polar CD = 0.024 + 0.0375 CL^2;
-> "[CL]" (math calc "2*[W]/([rho]*[V]^2*[S])");
-> "[CD]" (math calc "0.024+0.0375*[CL]^2");
#3. Speed of minimum drag D = W*CD/CL as a function of V;
-> "[D]" (concat "0.5*" [rho] "*V^2*" [S] "*(0.024+0.0375*(2*" [W] "/(" [rho] "*V^2*" [S] "))^2)");
-> "[best]" (num minimize [D] "V" 100 300);
json load "rx_drag" [best];
concat "CL = " (math round [CL] 3) ", L/D = " (math round (math calc "[CL]/[CD]") 2) ", thrust required = " (math round (math calc "[W]*[CD]/[CL]/1000") 1) " kN, min drag speed = " (math round (unit convert (json select "rx_drag" /x) "m/s" "kn") 0) " kn TAS"

Result

mentdb

CL = 0.509, L/D = 15.1, thrust required = 41.6 kN, min drag speed = 357 kn TAS

Go further

Find the speed of best range: maximise V/D
Plot the drag curve with plot function

executed on MentDB

Best glide of a sailplane from its measured polar

Fit a quadratic polar to the sink rates measured in flight, find the speed of best glide and the glide ratio, then the still-air range from 1000 m.

Pipeline

01 · fitQuadratic polarstat fit_model→02 · optimiseBest glide speednum maximize→03 · rangeRange in km and nmiunit convert

Script

admin

#1. Measured polar of a sailplane: speed in km/h, sink rate in m/s;
-> "[v]" "[70, 80, 90, 100, 120, 140, 160]";
-> "[w]" "[0.62, 0.60, 0.63, 0.70, 0.93, 1.25, 1.68]";
#2. Quadratic polar fit w = a*x^2+b*x+c;
-> "[fit]" (stat fit_model "a*x^2+b*x+c" "x" [v] [w] "{\"a\": 0.0001, \"b\": -0.02, \"c\": 1.5}");
json load "rx_polar" [fit];
-> "[polar]" (json select "rx_polar" /expression);
#3. Best glide: maximise the glide ratio x/3.6/w;
-> "[best]" (num maximize (concat "x/3.6/(" [polar] ")") "x" 60 160);
json load "rx_best" [best];
-> "[LD]" (json select "rx_best" /value);
#4. Still air range from 1000 m;
concat "best glide at " (math round (json select "rx_best" /x) 1) " km/h, L/D = " (math round [LD] 1) ", range from 1000 m = " (math round (unit convert (math calc "[LD]*1000") "m" "km") 1) " km = " (math round (unit convert (math calc "[LD]*1000") "m" "nmi") 1) " nmi"

Result

mentdb

best glide at 97.8 km/h, L/D = 39.4, range from 1000 m = 39.4 km = 21.3 nmi

Go further

Add a headwind: maximise (V - Vw)/w
MacCready speed to fly with a thermal climb rate

executed on MentDB

Parabolic drag polar from wind tunnel data

Wind tunnel points of a wing (aspect ratio 8) are fitted with CD = CD0 + k CL^2. The script gives CD0 and k with their standard errors, the Oswald efficiency, the best lift-to-drag ratio and its CL, and the extrapolated drag at CL = 1.

Pipeline

01 · loadLoad CL^2, CD pairspa rl load_from_json→02 · fitCD0 and kpa rl intercept→03 · errorsStandard errorspa rl slope_std_err→04 · predictCD at CL = 1pa rl predict

Script

admin

#1. Wind tunnel drag polar of a wing AR = 8: points CL^2, CD;
pa rl load_from_json "rx_polar" "[[0.04, 0.0216], [0.09, 0.0242], [0.16, 0.0271], [0.25, 0.0314], [0.36, 0.0360], [0.49, 0.0423], [0.64, 0.0487]]";
#2. Parabolic polar CD = CD0 + k CL^2;
-> "[CD0]" (pa rl intercept "rx_polar");
-> "[k]" (pa rl slope "rx_polar");
#3. Oswald factor, best lift to drag ratio and its CL;
-> "[e]" (math calc "1/(pi*8*[k])");
-> "[LDmax]" (math calc "1/(2*sqrt([CD0]*[k]))");
#4. Extrapolated drag at CL = 1.0 and uncertainty of the fit;
-> "[CD1]" (pa rl predict "rx_polar" 1.0);
-> "[r]" (concat "{\"exists\": " (pa rl exist "rx_polar") ", \"CD0\": " (math round [CD0] 5) ", \"CD0_std_err\": " (math round_sig (pa rl intercept_std_err "rx_polar") 2) ", \"k\": " (math round [k] 4) ", \"k_std_err\": " (math round_sig (pa rl slope_std_err "rx_polar") 2) ", \"oswald_e\": " (math round [e] 3) ", \"L_D_max\": " (math round [LDmax] 1) ", \"CL_best\": " (math round (math calc "sqrt([CD0]/[k])") 3) ", \"CD_at_CL1\": " (math round [CD1] 4) "}");
pa rl close "rx_polar";
[r]

Result

mentdb

{"exists": 1, "CD0": 0.01996, "CD0_std_err": 0.00013, "k": 0.0451, "k_std_err": 0.00038, "oswald_e": 0.882, "L_D_max": 16.7, "CL_best": 0.665, "CD_at_CL1": 0.0651}

Go further

Convert the best CL to a best glide speed with unit calc and the wing loading.
Check the residuals: a polar that curves up at high CL shows the onset of separation.

executed on MentDB

Take-off ground roll of a twin-jet

A 70 t twin-jet takes off at sea level with 210 kN of thrust. Compute the stall and lift-off speeds, then integrate the equation of motion with aerodynamic drag, lift and rolling friction to get the ground roll distance and time.

Pipeline

01 · constantStandard gravityconst get→02 · speedsVS and VLOF = 1.15 VSmath calc→03 · integrates = integral of m V / F dVnum integrate→04 · convertSpeeds in knotsunit convert

Script

admin

#1. Twin-jet at MTOW, sea level ISA: mass, average thrust, wing area;
-> "[m]" 70000; -> "[T]" 210000; -> "[S]" 122.6; -> "[rho]" 1.225;
-> "[g]" (const get "g0");
#2. Ground roll coefficients: CL and CD in take-off configuration, rolling friction;
-> "[CL]" 0.8; -> "[CD]" 0.06; -> "[mu]" 0.02;
#3. Lift-off speed VLOF = 1.15 VS with CLmax = 2.4;
-> "[VS]" (math calc "sqrt(2*[m]*[g]/([rho]*[S]*2.4))");
-> "[VLOF]" (math calc "1.15*[VS]");
#4. Ground roll s = integral of m V / F dV and time t = integral of m / F dV;
-> "[q]" (math calc "0.5*[rho]*[S]");
-> "[F]" (concat "(" [T] "-" [q] "*V^2*" [CD] "-" [mu] "*(" (* [m] [g]) "-" [q] "*V^2*" [CL] "))");
-> "[s]" (num integrate (concat [m] "*V/" [F]) "V" 0 [VLOF]);
-> "[time]" (num integrate (concat [m] "/" [F]) "V" 0 [VLOF]);
concat "VS = " (math round (unit convert [VS] "m/s" "kn") 1) " kt, VLOF = " (math round (unit convert [VLOF] "m/s" "kn") 1) " kt | ground roll = " (math round [s] 0) " m in " (math round [time] 1) " s | mean acceleration = " (math round (/ [VLOF] [time]) 2) " m/s2"

Result

mentdb

VS = 120 kt, VLOF = 138 kt | ground roll = 939 m in 26.1 s | mean acceleration = 2.72 m/s2

Go further

Make the thrust decrease with the speed, T = T0 - k V, and compare.
Solve with num solve the runway length needed for a hot and high airport (lower rho).

executed on MentDB

Wind triangle of a flight plan leg and runway wind components

The pilot knows the true course, the true airspeed and the forecast wind. The script solves the wind triangle for the wind correction angle and the heading, checks the ground speed and the track by adding the air and wind vectors, then splits the wind into headwind and crosswind on the destination runway.

Pipeline

01 · correctionWind correction anglemath calc→02 · speedGround speedmath hypot→03 · trackTrack checkmath atan2→04 · runwayHead and cross windmath abs→05 · sideLeft or rightif

Script

admin

#1. Leg of a flight plan: true course 045, TAS 120 kt, wind from 310 at 25 kt;
-> "[tc]" (math deg_to_rad 45); -> "[tas]" 120; -> "[wd]" (math deg_to_rad 310); -> "[ws]" 25;
#2. Wind correction angle and true heading;
-> "[wca]" (math calc "asin([ws]/[tas]*sin([wd]-[tc]))");
-> "[hdg]" (math calc "[tc]+[wca]");
#3. Check: air vector plus wind vector, the wind blows toward wd + 180;
-> "[n]" (math calc "[tas]*cos([hdg])+[ws]*cos([wd]+pi)");
-> "[e]" (math calc "[tas]*sin([hdg])+[ws]*sin([wd]+pi)");
-> "[gs]" (math hypot [e] [n]);
-> "[trk]" (math rad_to_deg (math atan2 [e] [n]));
#4. Components on runway 36 at the destination;
-> "[head]" (math calc "[ws]*cos([wd])"); -> "[cross]" (math calc "[ws]*sin([wd])");
concat "WCA = " (math round (math rad_to_deg [wca]) 1) " deg, heading = " (math round (math rad_to_deg [hdg]) 1) ", ground speed = " (math round [gs] 1) " kt, track check = " (math round [trk] 2) ", runway 36: headwind " (math round [head] 1) " kt, crosswind " (math round (math abs [cross]) 1) " kt from the " (if (> [cross] 0) {"right"} {"left"})

Result

mentdb

WCA = -12 deg, heading = 33, ground speed = 119.6 kt, track check = 45, runway 36: headwind 16.1 kt, crosswind 19.2 kt from the left

Go further

Loop over the legs of the flight plan stored in a JSON array with json load and sum the times.
Compare the crosswind with the demonstrated crosswind of the aircraft.

executed on MentDB

Load factor, radius and rate of a level turn

For a light aircraft in a 45 deg banked turn, the script gets the load factor with the secant of the bank angle, the stall speed in the turn, the turn radius and rate, and the bank angle of a standard rate one turn at the same speed.

Pipeline

01 · speedKnots to m/sunit convert→02 · loadn = sec phimath sec→03 · radiusR = v^2 / g tan phimath tan→04 · standardRate one bankmath atan

Script

admin

#1. Level turn at 140 kt TAS with 45 deg of bank, stall speed 61 kt in level flight;
-> "[v]" (unit convert 140 "kn" "m/s"); -> "[phi]" (math deg_to_rad 45); -> "[g]" 9.80665;
#2. Load factor n = sec phi and stall speed in the turn;
-> "[n]" (math sec [phi]);
-> "[vs]" (math calc "61*sqrt([n])");
#3. Radius v^2 / g tan phi, rate of turn and time for 360 deg;
-> "[R]" (math calc (concat "[v]^2/([g]*" (math tan [phi]) ")"));
-> "[rate]" (math rad_to_deg (math calc "[v]/[R]"));
#4. Bank angle of a rate one turn, 3 deg/s, at the same speed;
-> "[bank1]" (math rad_to_deg (math atan (math calc "[v]*3*pi/180/[g]")));
concat "n = " (math round [n] 3) " g, stall in turn = " (math round [vs] 1) " kt, radius = " (math round [R] 0) " m = " (math round (unit convert [R] "m" "nmi") 2) " nmi, rate = " (math round [rate] 2) " deg/s, 360 deg in " (math round (math calc "360/[rate]") 0) " s, rate one bank = " (math round [bank1] 1) " deg"

Result

mentdb

n = 1.414 g, stall in turn = 72.5 kt, radius = 529 m = 0.29 nmi, rate = 7.8 deg/s, 360 deg in 46 s, rate one bank = 21 deg

Go further

Tabulate n and the stall speed for banks from 0 to 70 deg with num table.
Add the limit load factor 3.8 g and find the bank where the stall speed reaches Va.

executed on MentDB

Breguet range, endurance and trip fuel of an airliner

With the cruise speed, the lift-to-drag ratio and the specific fuel consumption of a single-aisle jet, the script gets the Breguet range and endurance (ln W0/W1 as log1p, accurate for small fuel fractions) and the inverse problem: the fuel burnt on a 3000 km leg with expm1.

Pipeline

01 · massesWeights at start and endmath calc→02 · rangeln W0/W1 with log1pmath log1p→03 · fuelTrip fuel with expm1math expm1→04 · unitsnmi and hoursunit convert

Script

admin

#1. Single-aisle jet at Mach 0.78 and FL350: TAS 230 m/s, L/D 17, TSFC 0.60 per hour;
-> "[V]" 230; -> "[LD]" 17; -> "[c]" (math calc "0.60/3600");
#2. Masses: take-off 73500 kg, 18000 kg of fuel, 2500 kg kept as reserve;
-> "[W1]" (math calc "73500-18000+2500");
#3. Breguet range and endurance, ln W0/W1 computed as log1p of burnt fuel / W1;
-> "[ln]" (math log1p (math calc "(18000-2500)/[W1]"));
-> "[R]" (math calc "[V]/[c]*[LD]*[ln]");
-> "[E]" (math calc "[LD]/[c]*[ln]");
#4. Fuel burnt on a 3000 km leg: W1 expm1 of R g c / V L/D;
-> "[fuel]" (math calc (concat "[W1]*" (math expm1 (math calc "3000e3*[c]/([V]*[LD])"))));
concat "range = " (math round (math calc "[R]/1000") 0) " km = " (math round (unit convert [R] "m" "nmi") 0) " nmi, endurance = " (math round (unit convert [E] "s" "h") 2) " h, fuel for 3000 km = " (math round [fuel] 0) " kg"

Result

mentdb

range = 5556 km = 3000 nmi, endurance = 6.71 h, fuel for 3000 km = 7912 kg

Go further

Find with num maximize the speed of maximum range for a parabolic drag polar.
Compare with the step climb: L/D and TSFC change with the altitude.

executed on MentDB

Mach, true and calibrated airspeed from the air data computer

From the impact pressure, the static pressure and the total air temperature measured in cruise, the script gets the Mach number with the isentropic pitot formula, the static temperature, the true and calibrated airspeeds, and solves the Rayleigh pitot formula for a supersonic case.

Pipeline

01 · machIsentropic pitot formulamath power→02 · temperatureSAT from TATmath calc→03 · casCalibrated airspeedunit convert→04 · supersonicRayleigh pitot formulanum solve

Script

admin

#1. Air data computer: impact pressure 118 hPa at a static pressure of 238.4 hPa, FL350, total air temperature 249 K;
-> "[p]" 238.4; -> "[qc]" 118; -> "[tat]" 249;
#2. Subsonic Mach from the isentropic pitot formula, power 2/7;
-> "[M]" (math calc (concat "sqrt(5*(" (math power (math calc "[qc]/[p]+1") (math calc "2/7")) "-1))"));
#3. Static temperature, speed of sound and true airspeed;
-> "[sat]" (math calc "[tat]/(1+0.2*[M]^2)");
-> "[tas]" (math calc "[M]*sqrt(1.4*287.053*[sat])");
#4. Calibrated airspeed: same formula with the sea level pressure 1013.25 hPa;
-> "[cas]" (math calc (concat "340.294*sqrt(5*(" (math power (math calc "[qc]/1013.25+1") (math calc "2/7")) "-1))"));
#5. Supersonic case: a pitot ratio of 4.5 behind the normal shock, Rayleigh formula solved for M;
-> "[M2]" (num solve "(5.76*x^2/(5.6*x^2-0.8))^3.5*(2.8*x^2-0.4)/2.4-4.5" "x" 1 3);
concat "Mach = " (math round [M] 4) ", SAT = " (math round (unit convert [sat] "K" "degC") 1) " degC, TAS = " (math round (unit convert [tas] "m/s" "kn") 1) " kt, CAS = " (math round (unit convert [cas] "m/s" "kn") 1) " kt, supersonic Mach for pt2/p1 = 4.5: " (math round [M2] 4)

Result

mentdb

Mach = 0.7802, SAT = -51.2 degC, TAS = 453 kt, CAS = 264.5 kt, supersonic Mach for pt2/p1 = 4.5: 1.7627

Go further

Add the probe recovery factor of about 0.98 in the SAT formula.
Get the altitude from the static pressure with the ISA model, then the equivalent airspeed.

executed on MentDB

Lift of a Joukowski airfoil with complex numbers

A Joukowski airfoil is the image of a circle by zeta = z + 1/z. The script works in complex numbers to get the circle radius, the zero lift angle, the leading edge and the chord of the airfoil, then the lift coefficient from the Kutta condition and the lift slope.

Pipeline

01 · circlea = abs 1 - mucomplex abs→02 · angleZero lift anglecomplex arg→03 · edgeImage of the leading pointcomplex div→04 · liftKutta conditionmath calc

Script

admin

#1. Joukowski airfoil: circle of centre mu = -0.08+0.10i through the point z = 1, angle of attack 4 deg;
-> "[mu]" "-0.08+0.10i"; -> "[alpha]" (math deg_to_rad 4);
#2. Radius a = abs 1 - mu and zero lift angle beta = -arg 1 - mu;
-> "[d]" (complex sub 1 [mu]);
-> "[a]" (complex abs [d]);
-> "[beta]" (math calc (concat "-1*" (complex arg [d])));
#3. Leading edge: image of the circle point mu - a by zeta = z + 1/z;
-> "[zle]" (complex sub [mu] [a]);
-> "[le]" (complex re (complex add [zle] (complex div 1 [zle])));
-> "[chord]" (math calc (concat "2-(" [le] ")"));
#4. Kutta condition: Gamma = 4 pi V a sin alpha + beta, lift coefficient CL = 8 pi a sin alpha + beta / chord;
-> "[CL]" (math calc "8*pi*[a]*sin([alpha]+[beta])/[chord]");
concat "a = " (math round [a] 4) ", zero lift angle = " (math round (math rad_to_deg [beta]) 2) " deg, chord = " (math round [chord] 4) ", CL = " (math round [CL] 4) ", lift slope = " (math round (math calc "8*pi*[a]*cos([alpha]+[beta])/[chord]*pi/180") 4) " per deg"

Result

mentdb

a = 1.0846, zero lift angle = 5.29 deg, chord = 4.017, CL = 1.0955, lift slope = 0.1169 per deg

Go further

Draw the airfoil: map 100 points of the circle with complex from_polar and plot them with plot xy.
Compare the lift slope with the thin airfoil theory, 2 pi per radian.