Validating Simulations (NESC Check Cases 2 & 3)
Following the successful validation of translational dynamics in Check Case 1, the next step is to validate rotational EOM and inertial coupling effects.
NESC Check Cases 2 and 3 are designed specifically to expose errors in:
- angular rate propagation
- moment application
- inertia tensor handling
- cross-axis coupling effects
These are the parts of a simulator that can silently fail while still producing “reasonable-looking” trajectories.
In Case 1, we dropped a sphere with constant mass without considering drag.
For Cases 2 and 3, we drop a brick with initial angular rates. Case 2 will isolate rigid body rotation dynamics so no aerodynamic drag or damping effects will be considered. Case 3 builds on this by introducing aerodynamic damping, making inertial coupling effects more observable under dissipative dynamics.
Load Model for Case 2 and 3
The load models for Case 2 and Case 3 share the same structure, but differ in their coefficients. In the codebase, I have separated them to make it explicitly clear whether we are using Case 2 or Case 3 loads. The example code shown below is the general CaseLoad struct for both Case 2 and Case 3.
- Case 2: all damping coefficients are zero = no external moments
- Case 3: non-zero damping coefficients = rate-dependent aerodynamic moments
Forces are zero in both cases, so translational motion is unaffected by aerodynamics.
namespace sim::loads {
template <typename Vehicle> struct CaseLoads {
types::ExternalLoads compute(Vehicle const& vehicle, types::State const& x,
types::FlightCondition const& fc,
types::AuxData* aux = nullptr) const {
types::ExternalLoads loads{};
double const alpha_rad = fc.alpha_rad;
double const beta_rad = fc.beta_rad;
double const qbar = fc.dynamic_pressure;
double const Aref = vehicle.A_ref_m2;
double const c_alpha = std::cos(alpha_rad);
double const s_alpha = std::sin(alpha_rad);
double const c_beta = std::cos(beta_rad);
double const s_beta = std::sin(beta_rad);
double C_l =
vehicle.Clp * x(3) * vehicle.b_m / (2 * fc.true_airspeed_mps)
+ vehicle.Clr * x(5) * vehicle.b_m / (2 * fc.true_airspeed_mps);
double C_m =
vehicle.Cmq * x(4) * vehicle.c_m / (2 * fc.true_airspeed_mps);
double C_n =
vehicle.Cnp * x(3) * vehicle.b_m / (2 * fc.true_airspeed_mps)
+ vehicle.Cnr * x(5) * vehicle.b_m / (2 * fc.true_airspeed_mps);
loads.fx_b_N = 0.0;
loads.fy_b_N = 0.0;
loads.fz_b_N = 0.0;
loads.mx_b_Nm = C_l * qbar * Aref * vehicle.b_m;
loads.my_b_Nm = C_m * qbar * Aref * vehicle.c_m;
loads.mz_b_Nm = C_n * qbar * Aref * vehicle.b_m;
return loads;
}
};
}
You may be wondering how Case 2 and Case 3 are different if the load models are identical. The difference comes from the vehicle model. Case 2 does not consider aerodynamic drag nor damping, while Case 3 considers aerodynamic damping.
We can adjust the damping coefficients in the vehicle model to control whether rate-dependent aerodynamic moments are present.
To keep this post from getting too long with code, I will just show the vehicle model for Case 2. We can set $C_{l_p} = C_{m_q} = C_{n_r} = -1$ to switch to the vehicle model used for Case 3.
To reduce the chance of errors, I will separate the vehicle models used for Case 2 and Case 3 into their own files. Case 2 will use BrickVehicle and Case 3 will use BrickVehicleDamped.
namespace sim::vehicles {
struct BrickVehicle {
static constexpr double in2m = 0.0254;
static constexpr double slug2kg = 14.5939;
static constexpr double kg2slug = 1 / slug2kg;
static constexpr double ft2m = 0.304878;
static constexpr double m_brick_slug = 0.1554048;
static constexpr double m_kg = kg2slug * m_brick_slug;
static constexpr double Jxx_slugft2 = 0.00189422;
static constexpr double Jxx_kgm2 = slug2kg * (ft2m * ft2m) * Jxx_slugft2;
static constexpr double Jyy_slugft2 = 0.00621102;
static constexpr double Jyy_kgm2 = slug2kg * (ft2m * ft2m) * Jyy_slugft2;
static constexpr double Jzz_slugft2 = 0.00719467;
static constexpr double Jzz_kgm2 = slug2kg * (ft2m * ft2m) * Jzz_slugft2;
static constexpr double Jxz_slugft2 = 0;
static constexpr double Jxz_kgm2 = slug2kg * (ft2m * ft2m) * Jxz_slugft2;
static constexpr double length_brick_m = 8 * in2m;
static constexpr double width_brick_m = 4 * in2m;
static constexpr double A_ref_m2 = length_brick_m * width_brick_m;
static constexpr double b_m = 0.33333 * ft2m; // Reference wing span
static constexpr double c_m = 0.66667 * ft2m; // Reference wing chord
static constexpr double Clp = 0.0; // Roll damping from roll rate
static constexpr double Clr = 0.0; // Roll damping from yaw rate
static constexpr double Cmq = 0.0; // Pitch damping from pitch rate
static constexpr double Cnp = 0.0; // Yaw damping from roll rate
static constexpr double Cnr = 0.0; // Yaw damping from yaw rate
};
}
What Should The Simulation Results Look Like?
We can build an intuition about what should physically happen, if we were to perform this experiment in real life.
Case 2 (no damping, no drag)
With no external moments, the motion is governed entirely by rigid-body dynamics.
Due to inertial coupling, angular rates do not remain constant. Instead, rotation about one axis induces motion in others, producing oscillatory behavior while conserving rotational energy.
Case 3 (damped, no drag)
When damping is introduced, the angular rates decay over time.
The system still exhibits cross-axis coupling, but the motion is gradually damped as energy is dissipated through the aerodynamic moment terms.
Simulation Results
Case 2 Results:
Case 3 Results:
The simulator matches the expected NESC trends across angular rates, Euler angles, and translational velocities.
No artificial energy growth or numerical instability was observed, indicating that the integration and state propagation remain stable under coupled rotational dynamics.
Discussion
Like before, the discrepancies in the y-velocity in the NED frame can safely be ignored, as this is most likely due to the difference in the geodesy (i.e., flat-earth vs WGS-84)
At first glance, the x-velocity in the NED frame for Case 2 looks like it is diverging from the reference data. But you will see that the vertical scale is on the order of magnitude of $10^{-8}$ which is practically 0, thus we can attribute this to numerical precision limits (floating-point error).
Case 3 builds on Case 2 by introducing rate-dependent damping, making inertial coupling effects more observable.
This case is particularly useful for validating:
- correct implementation of cross-axis coupling terms
- proper application of aerodynamic moment coefficients
- stability of the rotational dynamics under dissipative effects
Discrepancies here are harder to diagnose because errors can originate from multiple layers. The agreement with reference data suggests:
- aerodynamic moment terms are computed correctly
- coupling between translation and rotation is correct
- no major frame or sign errors exist
Next Steps
With Check Cases 2 and 3 validated, the simulator now has:
- verified translational dynamics
- verified rotational dynamics
- verified inertial coupling effects
The next steps are:
- implement the drone dynamic model
- develop a rendering pipeline to visualize the simulation
Passing Check Cases 2 and 3 is a major milestone.
At this point, the simulator is reproducing physically correct rigid-body dynamics under both inertial and aerodynamic conditions.
This forms a solid foundation for building higher-level autonomy algorithms with confidence that the underlying physics are sound!
Enjoy Reading This Article?
Here are some more articles you might like to read next: