Validating Simulations
When building a physics simulator, the most important milestone is not visualization or performance. It is verification. Before trusting a model to simulate real systems, the implementation must be validated against known reference solutions.
I recently found Ben Dickinson, a great resource on YouTube for learning GNC concepts! He created a whole series on creating and validating simulations for aerospace applications.
His implementation is in Python, so I followed along while writing everything in C++, keeping in mind that my project has different objectives. Lots of credit to him for helping me get started with my own simulation. Please go check him out!
I completed the first major validation step for my autonomous drone simulation framework by reproducing NESC Atmospheric Check Case 1 from the NASA Engineering and Safety Center 6-DOF Flight Simulation Check Cases. The results match the expected behavior, confirming that the core physics implementation is functioning correctly.
This post summarizes the validation setup, implementation approach, and results.
Why Validation Matters
A flight simulator contains many interacting pieces:
- rigid-body equations of motion
- coordinate transformations
- atmosphere models
- numerical integration
- aerodynamic force modeling
A small mistake in any of these components can produce simulations that look plausible but are physically incorrect.
The NESC Check Cases provide reference scenarios specifically designed to detect implementation errors in flight dynamics models. If your simulator reproduces the expected results, it is strong evidence that the core physics pipeline is implemented correctly.
Simulator Architecture
The simulator is written entirely in C++ and structured as a modular physics engine.
I decided to use C++ as the language because (1) I enjoy C++ and want to practice it more and (2) I want to eventually do real-time Hardware-in-the-Loop (HIL) simulations. I thought that I could write more performant code in C++ than in other languages.
The design separates:
- Rigid-body dynamics (EOM)
- Flight condition computation
- Aerodynamic load models
- Numerical integration
At the core is a generic rigid-body derivative function:
namespace sim::physics {
template <typename Vehicle>
types::State derivative(Vehicle const& vehicle,
types::State const& x,
types::ExternalLoads const& loads,
types::AuxData* aux = nullptr);
}
This function computes the time derivative of the 12-state flat-earth rigid-body model using:
- body forces
- body moments
- vehicle inertia
- gravitational acceleration
External forces and moments are generated by a separate load model object, allowing different validation cases and vehicles to reuse the same dynamics implementation.
namespace sim::loads {
template <typename Vehicle> struct LoadModel {
types::ExternalLoads compute(Vehicle const& vehicle, types::State const& x,
types::FlightCondition const& fc,
types::AuxData* aux = nullptr) const {
types::ExternalLoads loads{};
// External forces and moment equations here
return loads;
}
};
}
The numerical integration layer uses a templated 4th-Order Runge-Kutta (RK4) integrator. This design supports both:
- offline trajectory simulation
- real-time simulation loops
Atmosphere Model
The simulator uses a tabulated implementation of the US Standard Atmosphere 1976 model. Density, gravity, and speed of sound are obtained via interpolation tables.
double rho = physics::ussa1976_density(h_m);
double a = physics::ussa1976_v_sound(h_m);
double g = physics::gravity_lookup(h_m);
From these values the simulator computes standard aerodynamic flight conditions:
- true airspeed
- dynamic pressure
- Mach number
- angle of attack
- sideslip
These quantities are packaged in a FlightCondition struct and passed to the aerodynamic load model.
NESC Atmospheric Check Case 1
Check Case 1 represents a constant-mass sphere dropped from an initial position located 30,000 ft over the Equator/Prime Meridian intersection. The sphere has no interaction with the atmosphere, meaning there is no drag or aerodynamic damping.
Check Case 1, as described in NESC’s Technical Assessment Report, uses the ellipsoidal Earth model while my simulation assumes the flat Earth model. This means that there will be some discrepancies in the results.
This is due to the rotation of the Earth (or lack thereof in my case) which would cause the sphere to have a small amount of roll rate with respect to the rotating Earth. This is observed in the discrepancies in roll angle and Y velocities in the NED and inertial reference frames.
Key properties of the case:
- spherical vehicle
- constant mass
- flat-earth dynamics
- no drag
- no external moments
The load model to compute the external forces and moments is defined as:
namespace sim::loads {
template <typename Vehicle> struct LoadModel {
types::ExternalLoads compute(Vehicle const& vehicle, types::State const& x,
types::FlightCondition const& fc,
types::AuxData* aux = nullptr) const {
types::ExternalLoads loads{};
loads.fx_b_N = 0.0;
loads.fy_b_N = 0.0;
loads.fz_b_N = 0.0;
loads.mx_b_Nm = 0.0;
loads.my_b_Nm = 0.0;
loads.mz_b_Nm = 0.0;
return loads;
}
};
}
For Check Case 1, all forces and moments are equal to 0. Other loading cases can be simulated by incorporating the proper models for forces and moments for each axis.
Numerical Integration
The system is integrated using a fixed-step RK4 integrator. Euler is quick and simple, but RK4 is the one you trust when the answer actually matters.
png image Svchb; svg image tobi; English translation by Profywld
namespace sim::numerical_methods {
template <typename Model>
inline void rk4_step(Model const& model, types::State& x, double h_s,
types::AuxData* aux = nullptr) {
types::State const k1 = h_s * model.derivative(x, nullptr);
types::State const k2 = h_s * model.derivative(x + 0.5 * k1, nullptr);
types::State const k3 = h_s * model.derivative(x + 0.5 * k2, nullptr);
types::State const k4 = h_s * model.derivative(x + k3, nullptr);
x += (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;
if (aux != nullptr) {
(void)model.derivative(x, aux);
}
}
}
For the offline simulation, a wrapper around the rk4_step is used to store the entire vehicle state. The results are exported to CSV to be post-processed and analyzed.
namespace sim::numerical_methods {
template <typename Model>
inline void rk4(Model const& model,
Eigen::Matrix<double, 12, Eigen::Dynamic>& x_hist,
std::vector<double> const& t_s, double h_s,
std::vector<types::AuxData>* aux_log = nullptr) {
std::size_t const nt_s = t_s.size();
types::State x = x_hist.col(0);
for (std::size_t i = 1; i < nt_s; ++i) {
types::AuxData aux{};
rk4_step(model, x, h_s, aux_log ? &aux : nullptr);
x_hist.col(i) = x;
if (aux_log != nullptr) {
aux_log->push_back(aux);
}
}
}
}
Validation Results
The simulation reproduces the expected behavior of the NESC Check Case 1 reference values. Check Case 1 is primarily designed to validate gravitation and translational equations of motion.
I was also able to verify my implementation of the USSA 1976 Atmospheric Model and RK4 numerical integration methods. The plots are generated using Python.
Discussion of Results
Some key things to note are:
- Slight discrepancies in the atmospheric data e.g., speed of sound and air density.
- Discrepancy in the roll and pitch Euler angles.
- Discrepancy in the Y body velocity in the NED frame.
Regarding (1), these are likely small errors in the interpolation of the data. My implementation of the USSA 1976 Atmosphere model is a lookup table as a function of altitude with explicit values calculated 100 meters between 0-86 kilometers. I interpolate for values that do not fall exactly on the breakpoints.
However, this should not be an issue for the actual drone simulation, since the vehicle will operate at relatively constant altitudes where atmospheric properties such as air density change very little.
Items (2) and (3) can be explained by the fact that the NESC reference data uses the WGS-84 rotating Earth model, while my simulation uses a flat-Earth model.
The other difference is the pitch angle, which is the one discrepancy I cannot fully explain yet.
The NESC report details initial conditions of 0 degrees for initial pitch angle while the provided Initial_Conditions.xlsx provide -90 degrees for initial pitch.
Further investigation is needed to determine the validity of the Euler angle values due to discrepancies in the simulation results and contradictory initial condition values provided by the NESC.
Lessons Learned
Several implementation details were important during validation.
Separating Physics Layers
Keeping these components independent simplified debugging:
- flight condition computation
- aerodynamic load model
- rigid-body dynamics
Each part could be verified independently.
Avoiding Duplicate Computation
Flight conditions such as $ \alpha $, $ \beta $, and $\bar{q}$ are computed once and reused by the load model rather than recomputed inside the dynamics.
Designing for Real-Time Simulation
The integrator was structured around a single-state step function, allowing the same integrator to be used in both offline simulations and real-time loops.
Next Steps
With the first validation case complete, the next milestones are:
- implementing NESC Check Case 2 and 3
- extending aerodynamic models for state-dependent coefficients
Once the remaining Check Cases are verified, the simulator will serve as the physics backbone for the autonomous drone stack.
Closing Thoughts
Validation is an often overlooked but essential step when building physics-based systems. Passing the NESC Check Cases provides confidence that the simulator’s fundamental equations are correct before moving on to more complex vehicle models.
Once the baseline physics are validated, the focus can shift toward vehicle-specific dynamics and control development.
Enjoy Reading This Article?
Here are some more articles you might like to read next: