Aerospace starter - Aircraft performance and drag fit

MODEL "Aerospace starter - Aircraft performance and drag fit"
DESCRIPTION "Parabolic drag, analytic and bounded numerical optima, and a separate synthetic drag-polar fit."
VERSION "1.0.0"
TAGS "aerospace", "starter", "aircraft", "optimization"

strict

input Mass_kg = 1200
  VALIDATE BETWEEN 500 AND 2000 MESSAGE "Use 500 to 2000 kg in this synthetic light-aircraft case."
input Density_kg_m3 = 1.225
  VALIDATE BETWEEN 0.5 AND 1.3 MESSAGE "Use density from 0.5 to 1.3 kg/m^3."
input Speed_m_s = 50
  VALIDATE BETWEEN 20 AND 150 MESSAGE "Use 20 to 150 m/s; check the separate stall flag."
WingArea_m2 = 16.2
AspectRatio = 7.5
OswaldEfficiency = 0.8
ZeroLiftDrag = 0.025
MaxLiftCoefficient = 1.5
Weight_N = Mass_kg * 9.80665
InducedFactor = 1 / (PI() * AspectRatio * OswaldEfficiency)
DragAtSpeed = (speed) => 0.5 * Density_kg_m3 * WingArea_m2 * ZeroLiftDrag * speed^2 + 2 * InducedFactor * Weight_N^2 / (Density_kg_m3 * WingArea_m2 * speed^2)

## Stall estimate from the selected density, weight and fixed CLmax.
output StallSpeed_m_s = SQRT(2 * Weight_N / (Density_kg_m3 * WingArea_m2 * MaxLiftCoefficient))
output Drag_N = DragAtSpeed(Speed_m_s)
output Power_kW = Drag_N * Speed_m_s / 1000
output CruiseAboveStall = Speed_m_s >= StallSpeed_m_s
## Unconstrained analytic minimum-drag speed; assess the stall flag separately.
output AnalyticBestSpeed_m_s = SQRT(2 * Weight_N / (Density_kg_m3 * WingArea_m2) * SQRT(InducedFactor / ZeroLiftDrag))
## Scalar minimization on 20 to 150 m/s; no lift or envelope constraint is added by MINIMIZE.
output SolvedBestSpeed_m_s = MINIMIZE(DragAtSpeed, 20, 150, 1e-9, 200)
output BestSpeedAboveStall = SolvedBestSpeed_m_s >= StallSpeed_m_s
output MaximumLiftDrag = 1 / (2 * SQRT(ZeroLiftDrag * InducedFactor))

# The fit data are synthetic and are independent of the flight assumptions above.
LiftSamples = [0.2, 0.4, 0.6, 0.8, 1.0]
DragSamples = [0.0273, 0.0330, 0.0444, 0.0585, 0.0783]
FittedPolar = CURVE_FIT((lift, offset, induced) => offset + induced * lift^2, LiftSamples, DragSamples, [0.02, 0.04])
output FittedZeroLiftDrag = INDEX(FittedPolar, 1, 1)
output FittedInducedFactor = INDEX(FittedPolar, 1, 2)
DragCurve = MAP([30, 40, 50, 60, 70, 80, 90], (speed) => {speed: speed, drag: DragAtSpeed(speed)})
output Assumptions = "Steady level flight; fixed parabolic drag polar, density and wing geometry; no compressibility, power limit, trim, gust or maneuver analysis. Values use kg, m, s and N. Synthetic fit is not experimental aircraft validation."
output ScenarioGuide = "Save Baseline at 1200 kg and 50 m/s: drag 916.3266173 N; optimum about 41.56522388 m/s. Set Speed_m_s to 60 and inspect the updated drag and power. Select DragCurve and choose Plot selection, X=speed and Y=drag. Save, restore Baseline and reopen the .gridoc archive."

Charts_1!type = "charts"
Charts_1!config = <toml>
  version = 1
  kind = "charts"
  title = "Aircraft drag curve"
  [layout]
  x_axis = "Speed (m/s)"
  y_axis = "Drag (N)"
  [[charts]]
  type = "line"
  smooth = false
  bind = "DragCurve"
  x = "speed"
  y = "drag"
  label = "Seven computed samples; straight segments are presentation"
</toml>

END MODEL