Two-Stage Rocket: Gravity-Turn Ascent to Escape Velocity

Project.toml

@@ -1,33 +1,49 @@
 name = "JuliaCon2026Components"
 uuid = "cb0cf210-d080-4a98-89fc-0e4ae9a4086e"
-authors = [" <>"]
 version = "0.1.0"
+authors = [" <>"]
 
 [deps]
-Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
+BlockComponents = "1ef5d832-be8e-447f-9f2c-8b5d17d15fd3"
+DyadData = "cab12561-e36c-4498-ae70-c9a5e79fef2a"
 DyadInterface = "99806f68-afab-45ca-9d8c-ceff6bc61f54"
+ElectricalComponents = "3fc4c519-469d-4ff5-98c9-d4d06b1514a2"
+HydraulicComponents = "b6ec4bdb-baef-4449-8ab3-f7a9cc7f1342"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Moshi = "2e0e35c7-a2e4-4343-998d-7ef72827ed2d"
+MultibodyComponents = "01883e52-22cd-4538-b14d-b44f958a131d"
 OrdinaryDiffEqDefault = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+RotationalComponents = "824fa3ed-d7dc-4365-8bdf-ac77a1ac2011"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+ThermalComponents = "084a29d9-4f48-4925-95b4-b7d4c79725f4"
+TranslationalComponents = "42acc370-819b-4877-bdd6-eb3294461f5d"
 
 [compat]
-Markdown = "1"
+BlockComponents = "4.5.1"
+DyadData = "2.0.1"
+ElectricalComponents = "2.2.1"
+HydraulicComponents = "2.0.3"
 JSON = "0.21, 1"
+Markdown = "1"
+MultibodyComponents = "0.2.3"
+RotationalComponents = "2.5.4"
 TOML = "1"
+ThermalComponents = "2.0.5"
+TranslationalComponents = "2.5.0"
+
+[dyad]
+kernel = "3.3.0"
 
 [extras]
-DyadEcosystemDependencies = "7bc808db-8006-421e-b546-062440d520b7"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DyadEcosystemDependencies = "7bc808db-8006-421e-b546-062440d520b7"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DyadEcosystemDependencies","CSV","DataFrames","Plots","Test"]
-
-[dyad]
-kernel = "3.3.0"
+test = ["DyadEcosystemDependencies", "CSV", "DataFrames", "Plots", "Test"]

dyad/rocket.dyad

@@ -0,0 +1,183 @@
+"""
+Free-flying point body with time-varying mass in an inverse-square (point)
+gravity field. Translational dynamics only: orientation is held fixed at the
+identity so no rotational states are introduced.
+
+The instantaneous total mass is supplied through the `m_in` signal, so the
+staging / propellant logic can live in the enclosing model. External forces
+(thrust, drag) enter through `frame_a`; gravity is computed internally from the
+absolute position.
+
+Newton's law in the frame (with R = I, r_cm = 0) mirrors the standard
+`MultibodyComponents.Body`: `frame_a.f = m * (a_0 - g_0)`.
+"""
+component VariableMassBody
+  "Gravitational parameter of the central body [m^3/s^2]"
+  parameter mu::Real(units="m3/s2") = 3.986004418e14
+  frame_a = Frame3D()
+  "Instantaneous total mass"
+  m_in = RealInput()
+  "Absolute position resolved in the world frame"
+  variable r_0::Position[3]
+  "Absolute velocity resolved in the world frame"
+  variable v_0::Velocity[3]
+  "Absolute acceleration resolved in the world frame"
+  variable a_0::Acceleration[3]
+  "Gravitational acceleration resolved in the world frame"
+  variable g_0::Acceleration[3]
+  "Distance from the central body centre"
+  variable rmag::Length
+relations
+  r_0 = frame_a.r_0
+  v_0 = der(r_0)
+  a_0 = der(v_0)
+  # Fixed orientation: identity rotation, no torque
+  frame_a.R = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
+  frame_a.tau = [0, 0, 0]
+  # Inverse-square gravity pointing toward the origin
+  rmag = sqrt(r_0[1]^2 + r_0[2]^2 + r_0[3]^2)
+  g_0 = (-mu / rmag^3) * r_0
+  # Newton's second law with time-varying mass
+  frame_a.f = m_in * (a_0 - g_0)
+end
+
+"""
+Two-stage launch vehicle flying a gravity-turn ascent to escape velocity.
+
+The rocket is a `VariableMassBody` in a point-gravity field. Thrust and
+aerodynamic drag are applied through `MultibodyComponents.WorldForce` elements
+attached to the body frame. Propellant burn, staging (jettison of the spent
+launch vehicle when stage-1 propellant is depleted), the open-loop pitch
+program, and the exponential-atmosphere drag law are all expressed here.
+
+Coordinates: launch from `[0, R_body, 0]`, downrange toward `+x`, motion in the
+x-y plane. Escape is checked via the specific orbital energy
+`spec_energy = v^2/2 - mu/r >= 0`.
+"""
+component TwoStageRocket
+  # ---- Environment ----
+  "Gravitational parameter [m^3/s^2]"
+  parameter mu::Real = 3.986004418e14
+  "Central body radius [m]"
+  parameter R_body::Real = 6.371e6
+  "Sea-level atmospheric density [kg/m^3]"
+  parameter rho0::Real = 1.225
+  "Atmospheric scale height [m]"
+  parameter Hscale::Real = 8500.0
+  # ---- Stage 1 (launch vehicle) ----
+  "Stage-1 structural (dry) mass [kg]"
+  parameter m1_dry::Real = 18000.0
+  "Stage-1 propellant mass [kg]"
+  parameter m1_prop::Real = 320000.0
+  "Stage-1 exhaust velocity [m/s]"
+  parameter ve1::Real = 3100.0
+  "Stage-1 propellant mass-flow rate [kg/s]"
+  parameter mdot1::Real = 2200.0
+  "Full-stack reference (frontal) area [m^2]"
+  parameter A_stack::Real = 11.0
+  # ---- Stage 2 (orbiter) ----
+  "Stage-2 structural (dry) mass [kg]"
+  parameter m2_dry::Real = 3500.0
+  "Stage-2 propellant mass [kg]"
+  parameter m2_prop::Real = 80000.0
+  "Stage-2 exhaust velocity [m/s]"
+  parameter ve2::Real = 4400.0
+  "Stage-2 propellant mass-flow rate [kg/s]"
+  parameter mdot2::Real = 260.0
+  "Orbiter reference (frontal) area [m^2]"
+  parameter A_orbiter::Real = 5.0
+  # ---- Payload & aerodynamics ----
+  "Payload mass carried to the end [kg]"
+  parameter payload::Real = 1500.0
+  "Drag coefficient"
+  parameter Cd::Real = 0.4
+  # ---- Pitch program (open-loop gravity turn) ----
+  "Duration of the initial vertical rise [s]"
+  parameter t_vert::Real = 12.0
+  "Time at which the pitch-over ramp ends [s]"
+  parameter t_pitch_end::Real = 140.0
+  "Held pitch angle from local vertical [rad]"
+  parameter phi_max::Real = 1.05
+  # ---- Multibody components ----
+  world = MultibodyComponents.World()
+  body = VariableMassBody(mu = mu, r_0 = initial [0, R_body, 0], v_0 = initial [0, 0, 0])
+  thrust = MultibodyComponents.WorldForce(resolve_in_frame = MultibodyComponents.ResolveInFrame.World())
+  drag = MultibodyComponents.WorldForce(resolve_in_frame = MultibodyComponents.ResolveInFrame.World())
+  # ---- Propellant states ----
+  variable mp1::Real
+  variable mp2::Real
+  # ---- Diagnostics / intermediate quantities ----
+  variable m::Real
+  variable Fthrust::Real
+  variable phi::Real
+  variable rmag::Real
+  variable vmag::Real
+  variable altitude::Real
+  variable rho::Real
+  variable q::Real
+  variable A::Real
+  variable stage::Real
+  variable radial_hat::Real[3]
+  variable downrange_hat::Real[3]
+  variable thrust_dir::Real[3]
+  variable v_esc::Real
+  variable spec_energy::Real
+relations
+  initial mp1 = m1_prop
+  initial mp2 = m2_prop
+  # Propellant burn; staging by depletion of stage-1 propellant
+  der(mp1) = ifelse(mp1 > 0.0, -mdot1, 0.0)
+  der(mp2) = ifelse(mp1 > 0.0, 0.0, ifelse(mp2 > 0.0, -mdot2, 0.0))
+  # Thrust magnitude: stage-1, then stage-2, then coast
+  Fthrust = ifelse(mp1 > 0.0, mdot1 * ve1, ifelse(mp2 > 0.0, mdot2 * ve2, 0.0))
+  # Total mass; jettison the spent launch vehicle at staging
+  m = ifelse(mp1 > 0.0, m1_dry + mp1 + m2_dry + mp2 + payload, m2_dry + mp2 + payload)
+  # Reference area switches from full stack to orbiter at staging
+  A = ifelse(mp1 > 0.0, A_stack, A_orbiter)
+  stage = ifelse(mp1 > 0.0, 1.0, ifelse(mp2 > 0.0, 2.0, 0.0))
+  # Kinematics from the body
+  rmag = sqrt(body.r_0[1]^2 + body.r_0[2]^2 + body.r_0[3]^2)
+  vmag = sqrt(body.v_0[1]^2 + body.v_0[2]^2 + body.v_0[3]^2)
+  altitude = rmag - R_body
+  # Guidance: pitch program relative to the local vertical
+  radial_hat = body.r_0 / rmag
+  downrange_hat = [body.r_0[2], -body.r_0[1], 0] / rmag
+  phi = ifelse(time < t_vert, 0.0,
+        ifelse(time < t_pitch_end, phi_max * (time - t_vert) / (t_pitch_end - t_vert), phi_max))
+  thrust_dir = cos(phi) * radial_hat + sin(phi) * downrange_hat
+  # Atmosphere and dynamic pressure
+  rho = rho0 * exp(-altitude / Hscale)
+  q = 0.5 * rho * vmag^2
+  # Drive the multibody force elements and the body mass
+  body.m_in = m
+  thrust.force_x = Fthrust * thrust_dir[1]
+  thrust.force_y = Fthrust * thrust_dir[2]
+  thrust.force_z = Fthrust * thrust_dir[3]
+  drag.force_x = -Cd * 0.5 * rho * A * vmag * body.v_0[1]
+  drag.force_y = -Cd * 0.5 * rho * A * vmag * body.v_0[2]
+  drag.force_z = -Cd * 0.5 * rho * A * vmag * body.v_0[3]
+  connect(thrust.frame_b, body.frame_a)
+  connect(drag.frame_b, body.frame_a)
+  # Escape diagnostics
+  v_esc = sqrt(2 * mu / rmag)
+  spec_energy = 0.5 * vmag^2 - mu / rmag
+end
+
+analysis LaunchToEscape
+  extends TransientAnalysis(stop=2000, abstol=1e-6, reltol=1e-6)
+  model = TwoStageRocket()
+end
+
+"Orbit sanity check: constant-mass body on a circular orbit must keep constant radius."
+test component TestOrbit
+  world = MultibodyComponents.World()
+  mass = BlockComponents.Sources.Constant(k=1.0)
+  body = VariableMassBody(r_0 = initial [7.0e6, 0, 0], v_0 = initial [0, 7546.05, 0])
+relations
+  connect(mass.y, body.m_in)
+end
+
+analysis TestOrbitTransient
+  extends TransientAnalysis(stop=6000, abstol=1e-8, reltol=1e-8)
+  model = TestOrbit()
+end

generated/LaunchToEscape_definition.jl

@@ -0,0 +1,52 @@
+### DO NOT EDIT THIS FILE
+### This file is auto-generated by the Dyad command-line compiler.
+### If you edit this code it is likely to get overwritten.
+### Instead, update the Dyad source code and regenerate this file
+
+
+using DyadInterface
+using DyadInterface: ODEAlg, DEVerbosity, OptimizationLevel
+using ModelingToolkit: SymbolicT, toggle_namespacing
+using DyadInterface: AbstractTransientAnalysisSpec, TransientAnalysisSpec
+@kwdef mutable struct LaunchToEscapeSpec <: AbstractTransientAnalysisSpec
+  name::Symbol = :LaunchToEscape
+  var"alg"::ODEAlg.Type = ODEAlg.Auto()
+  var"start"::Float64 = 0
+  var"stop"::Float64 = 2000
+  var"abstol"::Float64 = 0.000001
+  var"reltol"::Float64 = 0.000001
+  var"saveat"::Float64 = 0
+  var"dtmax"::Float64 = 0
+  var"tstops"::Array{Float64, 1} = []
+  var"automatic_discontinuity_detection"::Bool = false
+  var"optimize"::OptimizationLevel.Type = OptimizationLevel.Aggressive()
+  var"progress"::Bool = true
+  var"respecialize"::Bool = false
+  var"verbose"::DEVerbosity.Type = DEVerbosity.Standard()
+  var"log_file"::String = ""
+  # Two-stage launch vehicle flying a gravity-turn ascent to escape velocity.
+  # 
+  # The rocket is a `VariableMassBody` in a point-gravity field. Thrust and
+  # aerodynamic drag are applied through `MultibodyComponents.WorldForce` elements
+  # attached to the body frame. Propellant burn, staging (jettison of the spent
+  # launch vehicle when stage-1 propellant is depleted), the open-loop pitch
+  # program, and the exponential-atmosphere drag law are all expressed here.
+  # 
+  # Coordinates: launch from `[0, R_body, 0]`, downrange toward `+x`, motion in the
+  # x-y plane. Escape is checked via the specific orbital energy
+  # `spec_energy = v^2/2 - mu/r >= 0`.
+  var"model"::Union{Nothing, System} = JuliaCon2026Components.TwoStageRocket(; name=:TwoStageRocket)
+end
+
+function DyadInterface.run_analysis(spec::LaunchToEscapeSpec)
+  overrides = Dict{SymbolicT, SymbolicT}()
+  no_namespace_model = toggle_namespacing(spec.model, false)
+
+  base_spec = TransientAnalysisSpec(;
+    name=:TransientAnalysis, overrides, alg=spec.alg, start=spec.start, stop=spec.stop, abstol=spec.abstol, reltol=spec.reltol, saveat=spec.saveat, dtmax=spec.dtmax, tstops=spec.tstops, automatic_discontinuity_detection=spec.automatic_discontinuity_detection, optimize=spec.optimize, progress=spec.progress, respecialize=spec.respecialize, verbose=spec.verbose, log_file=spec.log_file, model=spec.model
+  )
+  run_analysis(base_spec)
+end
+
+LaunchToEscape(;kwargs...) = run_analysis(LaunchToEscapeSpec(;kwargs...))
+export LaunchToEscape, LaunchToEscapeSpec

generated/TestOrbitTransient_definition.jl

@@ -0,0 +1,42 @@
+### DO NOT EDIT THIS FILE
+### This file is auto-generated by the Dyad command-line compiler.
+### If you edit this code it is likely to get overwritten.
+### Instead, update the Dyad source code and regenerate this file
+
+
+using DyadInterface
+using DyadInterface: ODEAlg, DEVerbosity, OptimizationLevel
+using ModelingToolkit: SymbolicT, toggle_namespacing
+using DyadInterface: AbstractTransientAnalysisSpec, TransientAnalysisSpec
+@kwdef mutable struct TestOrbitTransientSpec <: AbstractTransientAnalysisSpec
+  name::Symbol = :TestOrbitTransient
+  var"alg"::ODEAlg.Type = ODEAlg.Auto()
+  var"start"::Float64 = 0
+  var"stop"::Float64 = 6000
+  var"abstol"::Float64 = 1e-8
+  var"reltol"::Float64 = 1e-8
+  var"saveat"::Float64 = 0
+  var"dtmax"::Float64 = 0
+  var"tstops"::Array{Float64, 1} = []
+  var"automatic_discontinuity_detection"::Bool = false
+  var"optimize"::OptimizationLevel.Type = OptimizationLevel.Aggressive()
+  var"progress"::Bool = true
+  var"respecialize"::Bool = false
+  var"verbose"::DEVerbosity.Type = DEVerbosity.Standard()
+  var"log_file"::String = ""
+  # Orbit sanity check: constant-mass body on a circular orbit must keep constant radius.
+  var"model"::Union{Nothing, System} = JuliaCon2026Components.TestOrbit(; name=:TestOrbit)
+end
+
+function DyadInterface.run_analysis(spec::TestOrbitTransientSpec)
+  overrides = Dict{SymbolicT, SymbolicT}()
+  no_namespace_model = toggle_namespacing(spec.model, false)
+
+  base_spec = TransientAnalysisSpec(;
+    name=:TransientAnalysis, overrides, alg=spec.alg, start=spec.start, stop=spec.stop, abstol=spec.abstol, reltol=spec.reltol, saveat=spec.saveat, dtmax=spec.dtmax, tstops=spec.tstops, automatic_discontinuity_detection=spec.automatic_discontinuity_detection, optimize=spec.optimize, progress=spec.progress, respecialize=spec.respecialize, verbose=spec.verbose, log_file=spec.log_file, model=spec.model
+  )
+  run_analysis(base_spec)
+end
+
+TestOrbitTransient(;kwargs...) = run_analysis(TestOrbitTransientSpec(;kwargs...))
+export TestOrbitTransient, TestOrbitTransientSpec