Build a 5-DOF Rocket Simulation

Aerodynamics Subsystem

1. Feed the 3-axis body velocity Vb into a MATLAB Function.
2. The function outputs individual components (Vx, Vy, Vz) plus magnitude:
   
   

   function [Vx, Vy, Vz, Vmag, VxAbs] = mag(Vb)
   
   Vx = Vb(1);
   Vy = Vb(2);
   Vz = Vb(3);
   
   VxAbs = abs(Vx);
   Vmag = norm(Vb);
   
   end

3. Use two If blocks to guard against negative Vx/Vy when computing AoA and Beta.
4. Each branch contains a Gain and a Merge.
5. Convert AoA and Beta from radians to degrees (x180/pi) for intuition.
6. Bundle AoA, Beta, Vb and body rates w into the main aerodynamic-coefficients MATLAB Function block
7. Deliver Fb and Mb to the calling level through Outports.

Lift, Side-Force & Drag Function

1. The following MATLAB function computes aerodynamic forces (Flift, Fside, Fdrag) and moments (Mb).
2. It uses inputs: angle-of-attack (AoA), sideslip angle (Beta), body-axis velocity (Vb), and angular-rates (w).
3. Lift and side-force act perpendicularly, aligned to body frame axes.
4. Forces are summed and inverted to match coordinate conventions.

function [Fb, Mb, Flift, Fside, Fdrag, Cn, Cs, CmpDamp, CmyDamp] = Lift(AoA, Beta, Vb, w)
% LIFT calculates aerodynamic forces and moments based on AoA, sideslip,
% body-frame velocity, and angular rates. Used in Aero block of 6DOF rocket sim.

% Geometry and constants
Ab = 0.003504635;            % Reference area (body cross-section)
Afin = 0.00965;              % Fin area
Aplanform = 0.0835;          % Total lifting planform area
l = 0.202;                   % Lever arm from CoG to CoP
Mroll = 0;                   % Roll moment (not modeled here)
rho = 1.225;                 % Air density
Cdx = 0.65;                  % Base drag coefficient
d = 0.068;                   % Reference length (rocket diameter)
Vbx = Vb(1); Vby = Vb(2); Vbz = Vb(3); % Velocity components
wy = w(2); wz = w(3);        % Angular rates (pitch & yaw)
span = 0.1;                  % Fin span
midchord = (90 - 85.134)*pi/180; % Mid-chord sweep angle (rad)
K = 1.1;                     % Empirical lift correction factor

% Flow properties
AoAThreshold = 5*pi/180;     % Linear vs nonlinear lift region (5 deg)
Vmag = norm(Vb);             % Total airspeed
q = 0.5 * rho * Vmag^2;      % Dynamic pressure
M = Vmag / 340;              % Mach number
r = d / 2;                   % Body radius

if M > 1, M = 0.99; end      % Avoid sqrt(1-M^2) becoming complex

% Compressibility-adjusted aspect ratio factor (used in lift curve slope)
rho = span^2 * sqrt(1 - M^2) / (Afin * cos(midchord));

%% --- AoA-dependent Lift Coefficient (Cn) ---
if abs(AoA) < AoAThreshold
    % Linear region: only nose and fin effects
    CnNose = 2*sin(AoA)/Ab * pi*0.229^2;
    Cnfins = 3*(2*pi*span^2/Ab)/(1 + sqrt(1 + rho^2));
    CnAoA = CnNose + Cnfins;
    Cn = CnAoA * -sin(AoA);   % Negated for coordinate convention
else
    % Nonlinear region: includes body lift ~ sin^2(AoA)
    CnNose = 2*sin(AoA)/Ab * pi*0.229^2;
    Cnfins = 3*(2*pi*span^2/Ab)/(1 + sqrt(1 + rho^2));
    CnBodyLift = K * Aplanform / Ab * sin(AoA)^2;  % Nonlinear lift
    CnAoA = CnNose + Cnfins + CnBodyLift;
    Cn = CnAoA * -sin(AoA);   % Overall normal force coefficient
end

%% --- Beta-dependent Side Force Coefficient (Cs) ---
if abs(Beta) < AoAThreshold
    % Linear sideslip region
    CsNose = 2*sin(Beta)/Ab * pi*0.229^2;
    Csfins = 3*(2*pi*span^2/Ab)/(1 + sqrt(1 + rho^2));
    CsBeta = CsNose + Csfins;
    Cs = CsBeta * sin(Beta);
else
    % Nonlinear sideslip response: adds body lift effect
    CsNose = 2*sin(Beta)/Ab * pi*0.229^2;
    Csfins = 3*(2*pi*span^2/Ab)/(1 + sqrt(1 + rho^2));
    CsBodyLift = K * Aplanform / Ab * sin(Beta)^2;
    CsBeta = CsNose + Csfins + CsBodyLift;
    Cs = CsBeta * sin(Beta);
end

%% --- Pitch and Yaw Moment Coefficients ---
Cmp = (Cn * l) / d;          % Pitch moment from Cn
Cmy = (Cs * l) / d;          % Yaw moment from Cs

% Angular rate damping (dynamic stability): increases with rate^2
CmpDamp = -(1.25^4 * r * -wy^2) / (Ab * d * Vmag^2); % pitch rate damping
CmyDamp =  (1.25^4 * r * -wz^2) / (Ab * d * Vmag^2); % yaw rate damping

% Combine static and damping moments
Cmp = Cmp + CmpDamp;
Cmy = Cmy + CmyDamp;

%% --- Aerodynamic Forces ---
if Vmag > 1e-3
    % Drag force based on total velocity direction
    Fdragmagx = q * Ab * Cdx;
    Fdragx = Fdragmagx * (Vbx / Vmag * sign(Vbx));  % x-component
    Fdragyz = Fdragmagx * ([Vby * sign(Vby); Vbz * sign(Vbz)] / Vmag);  % y/z

    Fdrag = [Fdragx; Fdragyz];       % Total drag vector (body frame)
    Flift = q * Ab * Cn;             % Lift from AoA
    Fside = q * Ab * Cs;             % Side force from Beta
    Mpitch = q * Ab * d * Cmp;       % Pitching moment
    Myaw = q * Ab * d * Cmy;         % Yawing moment
else
    % Near-zero speed: no aerodynamic force
    Fdrag = [0;0;0]; Flift = 0; Fside = 0; Mpitch = 0; Myaw = 0;
end

%% --- Outputs ---
% Combine lift, side-force, and drag in body frame
Fb = [0; -Fside; Flift] - Fdrag;     % Net aerodynamic force (body frame)
Mb = [Mroll; -Mpitch; -Myaw];        % Net aerodynamic moment
end

Attitude Kinematics Subsystem


1. Bring the quaternion Q in through an Inport.
2. Square-and-sum it using a Math Function |u|² block.
3. A Compare to Zero block detects if the norm is zero and opens a bypass path to avoid NaNs.
4. The validated quaternion is passed into a MATLAB Function block titled q [R₁,R₂,R₃] (Rotation Order XYZ) which outputs Euler angles.
5. Insert a small Gain −1 to flip the sign convention before sending to the Outport Euler angles (feeding the 6dof simulation block which works best option for quick simulation results).
6. In parallel, feed Q into another MATLAB Function block that converts quaternion → axis-angle, outputting via axis angles Outport.
7. Also send Q into a Quaternion to DCM block and export the Direction Cosine Matrix through Outport DCM_be.

Integration / Stop-Condition Subsystem



Top path (attitude):

1. Integrate angular-rate derivative wd using an Integrator to get body-rates w.
2. Pass w into a second Integrator to accumulate Euler angle theta.
3. Attach Scope blocks to both integrators to visualise w and theta.

Bottom path (translation):

1. Integrate acceleration pdd once to obtain Velocity.
2. Integrate Velocity again to obtain Dis (displacement).
3. Route Vx and Xx into a nested Stop Condition subsystem.
4. Use two Relational Operator blocks to check the thresholds: “< −3 m/s” and “≤ −1.5 m”.
5. The results feed into an AND block that triggers a Stop Simulation block.
6. Multiplex the position components using a Mux block.
7. Apply a Gain 10 to scale the position signal for the visualiser (producing myTranslation. (I never got this working well).
8. If your inertial frame requires it, flip the sign of Xe using a simple Gain −1.

Download the Full Simulink Model

Get the complete 5-DOF rocket simulation system - ready to run in Simulink.

⬇ Download .SLX File