Aerospace starter - Vibration and flight events

MODEL "Aerospace starter - Vibration and flight events"
DESCRIPTION "Linear damped structural oscillator and a separate drag-free vertical trajectory. SI units; analytic verification."
VERSION "1.0.0"
AUTHOR "Aerospace evaluation"
TAGS "aerospace", "starter", "dynamics"

strict

# A representative elastic mode with 2 Hz natural frequency and 5% damping.
input Frequency_Hz = 2
  VALIDATE BETWEEN 0.1 AND 10 MESSAGE "Use a positive natural frequency from 0.1 to 10 Hz."
input DampingRatio = 0.05
  VALIDATE BETWEEN 0 AND 0.9 MESSAGE "This analytic reference is underdamped: use 0 to 0.9."
input InitialDisplacement_m = 0.01
  VALIDATE BETWEEN -0.1 AND 0.1 MESSAGE "Use an initial displacement between -0.1 and 0.1 m."
input EndTime_s = 2
  VALIDATE BETWEEN 0.01 AND 10 MESSAGE "Use an end time from 0.01 to 10 s; 101 output samples are retained."
Omega = 2 * PI() * Frequency_Hz
DampedOmega = Omega * SQRT(1 - DampingRatio^2)
Trajectory = ODE_SOLVE((time, state) => [INDEX(state, 1, 2), 0 - 2 * DampingRatio * Omega * INDEX(state, 1, 2) - Omega^2 * INDEX(state, 1, 1)], [InitialDisplacement_m, 0], 0, EndTime_s, 101, 1e-9, 1e-12)
AnalyticEnd = InitialDisplacement_m * EXP(0 - DampingRatio * Omega * EndTime_s) * (COS(DampedOmega * EndTime_s) + DampingRatio / SQRT(1 - DampingRatio^2) * SIN(DampedOmega * EndTime_s))
DisplacementCurve = MAP(SEQUENCE(101), (row) => {time: INDEX(Trajectory, row, 1), displacement: INDEX(Trajectory, row, 2)})
MaximumError = MAX(MAP(SEQUENCE(101), (row) => ABS(INDEX(Trajectory, row, 2) - InitialDisplacement_m * EXP(0 - DampingRatio * Omega * INDEX(Trajectory, row, 1)) * (COS(DampedOmega * INDEX(Trajectory, row, 1)) + DampingRatio / SQRT(1 - DampingRatio^2) * SIN(DampedOmega * INDEX(Trajectory, row, 1))))))

# Separate vertical flight: find the descending ground crossing.
input InitialHeight_m = 100
  VALIDATE BETWEEN 1 AND 1000 MESSAGE "Use 1 to 1000 m initial height."
input InitialVelocity_m_s = 50
  VALIDATE BETWEEN 0 AND 50 MESSAGE "Use an initial upward velocity from 0 to 50 m/s; the 30 s event interval brackets ground impact."
Gravity = 9.80665
GroundEvent = ODE_EVENT((time, state) => [INDEX(state, 1, 2), 0 - Gravity], [InitialHeight_m, InitialVelocity_m_s], 0, 30, (time, state) => INDEX(state, 1, 1), {event_direction: "falling", rtol: 1e-9, atol: 1e-12})
AnalyticImpact_s = (InitialVelocity_m_s + SQRT(InitialVelocity_m_s^2 + 2 * Gravity * InitialHeight_m)) / Gravity

# Numerical answers and validation quantities.
A1 = "Displacement at selected end time (m)"
output FinalDisplacement_m = INDEX(Trajectory, 101, 2)
A2 = "Velocity at selected end time (m/s)"
output FinalVelocity_m_s = INDEX(Trajectory, 101, 3)
A3 = "Analytic displacement (m)"
output AnalyticDisplacement_m = AnalyticEnd
A4 = "Maximum trajectory error (m)"
output MaximumDisplacementError_m = MaximumError
A5 = "Impact time (s)"
output ImpactTime_s = INDEX(GroundEvent, 1, 1)
A6 = "Impact height residual (m)"
output ImpactHeightResidual_m = INDEX(GroundEvent, 1, 2)
A7 = "Impact velocity (m/s)"
output ImpactVelocity_m_s = INDEX(GroundEvent, 1, 3)
A8 = "Analytic impact time (s)"
output AnalyticImpactTime_s = AnalyticImpact_s
A9 = "Impact time error (s)"
output ImpactTimeError_s = ABS(ImpactTime_s - AnalyticImpactTime_s)
A10 = "Maximum ballistic height (m)"
output MaximumBallisticHeight_m = InitialHeight_m + InitialVelocity_m_s^2 / (2 * Gravity)
output Assumptions = "Linear single-degree oscillator with zero initial velocity, constant natural frequency and damping, no forcing; rtol=1e-9, atol=1e-12, 101 output samples. The separate vertical flight uses constant Earth gravity and no drag. Agreement here does not validate aerodynamics or nonlinear structures. SI numerical working units."
output ScenarioGuide = "Save Baseline at 2 Hz and damping 0.05. Set DampingRatio to 0.1 and inspect DisplacementCurve. Plot selection: X=time, Y=displacement, with unsmoothed sample segments. Restore Baseline; maximum analytic displacement error should be below 1e-9 m and impact time about 11.909590765 s. Save and reopen .gridoc."

Charts_1!type = "charts"
Charts_1!config = <toml>
  version = 1
  kind = "charts"
  title = "Damped vibration"
  [layout]
  x_axis = "Time (s)"
  y_axis = "Displacement (m)"
  [[charts]]
  type = "line"
  smooth = false
  bind = "DisplacementCurve"
  x = "time"
  y = "displacement"
  label = "101 ODE output samples"
</toml>

END MODEL