Trajectory of a 3d body

Problem Statement: -

I have a 3d body which i have launched from cruising aircraft travelling at 200m/sec at 10km altitude. The launch speed is 50m/sec and it is launched at 45degrees angle pointing downwards. I want to create the whole trajectory of the body using everything like change in pressure, density, drag, lift etc. I created a code but it is not giving the accurate answer.

the code: -

function aircraft_launch_trajectory_6DOF_advanced()
    % Constants
    g = 9.81;  % Acceleration due to gravity (m/s^2)
    
    % Convert given units to SI
    m = 49.964 * 0.0283495;  % Mass (kg) (1 ouncemass = 0.0283495 kg)
    A = 68.166 * 0.00064516;  % Area (m^2) (1 in^2 = 0.00064516 m^2)
    Izz = 19.58 * 0.0283495 * 0.0254^2;  % Moment of inertia (kg*m^2)
    
    % Body dimensions
    length = 1.80 * 0.0254;  % Length (m)
    width = 3.30 * 0.0254;   % Width (m)
    height = 8.00 * 0.0254;  % Height (m)
    
    % Center of mass location (assuming from the base of the body)
    cm_z = 1.708 * 0.0254;  % (m)
    
    % Assumed center of pressure location (you may need to adjust this)
    cp_z = 4.5 * 0.0254;  % (m)
    
    % Distance between center of mass and center of pressure
    d_cp_cm = cp_z - cm_z;

    % Initial conditions
    h0 = 10000;  % Initial height (m)
    v0_aircraft = 200;  % Aircraft velocity (m/s)
    v0_launch = 50;     % Launch velocity (m/s)
    angle = -45;        % Launch angle (degrees)
    omega0 = 0;         % Initial angular velocity (rad/s)

    % Calculate initial velocity components
    v0x = v0_aircraft - v0_launch * cosd(angle);
    v0y = -v0_launch * sind(angle);

    % Initial state vector [x, y, vx, vy, theta, omega]
    y0 = [0; h0; v0x; v0y; deg2rad(angle); omega0];

    % Time settings
    tspan = [0 300];  % Time span for integration

    % ODE solver options
    options = odeset('Events', @ground_impact, 'RelTol', 1e-6, 'AbsTol', 1e-6);

    % Solve ODE using ode45 (4th order Runge-Kutta)
    [t, y] = ode45(@(t, y) equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length), tspan, y0, options);

    % Extract results
    x = y(:, 1);
    h = y(:, 2);
    vx = y(:, 3);
    vy = y(:, 4);
    theta = y(:, 5);
    omega = y(:, 6);

    % Calculate velocity magnitude
    v = sqrt(vx.^2 + vy.^2);

    % Plot trajectory
    figure;
    plot(x/1000, h/1000);
    xlabel('Horizontal distance (km)');
    ylabel('Altitude (km)');
    title('2D Trajectory of Body (Advanced 6DOF Simulation)');
    grid on;

    % Plot velocity
    figure;
    plot(t, v);
    xlabel('Time (s)');
    ylabel('Velocity (m/s)');
    title('Velocity vs Time');
    grid on;

    % Plot angular position and velocity
    figure;
    subplot(2,1,1);
    plot(t, rad2deg(theta));
    ylabel('Angular Position (degrees)');
    title('Angular Motion');
    grid on;
    subplot(2,1,2);
    plot(t, rad2deg(omega));
    xlabel('Time (s)');
    ylabel('Angular Velocity (deg/s)');
    grid on;

    % Display results
    disp(['Time of flight: ', num2str(t(end)), ' s']);
    disp(['Horizontal distance traveled: ', num2str(x(end)/1000), ' km']);
    disp(['Impact velocity: ', num2str(v(end)), ' m/s']);
    disp(['Final angle: ', num2str(rad2deg(theta(end))), ' degrees']);
    disp(['Max altitude: ', num2str(max(h)/1000), ' km']);
    disp(['Max velocity: ', num2str(max(v)), ' m/s']);

    % Output state at specific times for debugging
    debug_times = [0, 5, 10, 30, 60, 120];
    for debug_time = debug_times
        idx = find(t >= debug_time, 1);
        if ~isempty(idx)
            disp(['State at t = ', num2str(t(idx)), ' s:']);
            disp(['  Position: (', num2str(x(idx)/1000), ' km, ', num2str(h(idx)/1000), ' km)']);
            disp(['  Velocity: ', num2str(v(idx)), ' m/s']);
            disp(['  Angle: ', num2str(rad2deg(theta(idx))), ' degrees']);
        end
    end
end

function dydt = equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length)
    % Unpack state vector
    x = y(1);
    h = y(2);
    vx = y(3);
    vy = y(4);
    theta = y(5);
    omega = y(6);

    % Calculate velocity magnitude and angle of attack
    v = sqrt(vx^2 + vy^2);
    alpha = atan2(vy, vx) - theta;

    % Calculate air density and speed of sound using improved atmospheric model
    [rho, a] = atmosphere_model(h);

    % Calculate Mach number
    M = v / a;

    % Get aerodynamic coefficients based on angle of attack and Mach number
    [Cd, Cl, Cm] = aero_coefficients(alpha, M);

    % Calculate aerodynamic forces and moment
    q = 0.5 * rho * v^2;  % Dynamic pressure
    Fd = q * Cd * A;      % Drag force
    Fl = q * Cl * A;      % Lift force
    M_aero = q * Cm * A * length;  % Aerodynamic moment

    % Add moment due to center of pressure offset
    M_total = M_aero + Fl * d_cp_cm;

    % Apply initial boost for first 7 seconds
    boost_duration = 7;  % seconds
    if t <= boost_duration
        boost_force = m * 200 / boost_duration;  % Force to achieve 200 m/s in 7 seconds
        boost_fx = boost_force * cos(theta);
        boost_fy = boost_force * sin(theta);
    else
        boost_fx = 0;
        boost_fy = 0;
    end

    % Calculate accelerations
    ax = (-Fd * cos(theta) + Fl * sin(theta) + boost_fx) / m;
    ay = -9.81 + (-Fd * sin(theta) - Fl * cos(theta) + boost_fy) / m;
    alpha_dot = M_total / Izz;  % Angular acceleration

    % Return state derivatives
    dydt = [vx; vy; ax; ay; omega; alpha_dot];
end

function [rho, a] = atmosphere_model(h)
    % Constants
    R = 287.05;  % Specific gas constant for air (J/(kg*K))
    g0 = 9.80665;  % Standard gravity (m/s^2)
    gamma = 1.4;  % Ratio of specific heats for air

    % Sea level conditions
    T0 = 288.15;  % Temperature at sea level (K)
    P0 = 101325;  % Pressure at sea level (Pa)

    if h < 11000  % Troposphere
        T = T0 - 0.0065 * h;
        P = P0 * (T / T0)^(g0 / (R * 0.0065));
    elseif h < 20000  % Lower Stratosphere
        T = 216.65;
        P = 22632.1 * exp(-g0 * (h - 11000) / (R * T));
    else  % Upper Stratosphere
        T = 216.65 + 0.001 * (h - 20000);
        P = 5474.89 * (216.65 / T)^(34.1632);
    end

    rho = P / (R * T);
    a = sqrt(gamma * R * T);  % Speed of sound
end

function [Cd, Cl, Cm] = aero_coefficients(alpha, M)
    % This function returns aerodynamic coefficients based on angle of attack and Mach number
    % Note: These are simplified models and should be replaced with actual data or more complex models

    % Convert alpha to degrees for easier calculations
    alpha_deg = rad2deg(alpha);

    % Drag coefficient
    Cd0 = 0.1 + 0.1 * M;  % Base drag increases with Mach number
    Cdi = 0.01 * alpha_deg^2;  % Induced drag
    Cd = Cd0 + Cdi;

    % Lift coefficient
    Cl_slope = 0.1;  % Lift curve slope
    Cl = Cl_slope * alpha_deg * (1 - M^2/4)^(-0.5);  % Prandtl-Glauert correction for compressibility

    % Moment coefficient
    Cm0 = -0.05;  % Zero-lift moment coefficient
    Cm_alpha = -0.01;  % Moment curve slope
    Cm = Cm0 + Cm_alpha * alpha_deg;

    % Limit coefficients to realistic ranges
    Cd = max(0, min(Cd, 2));
    Cl = max(-2, min(Cl, 2));
    Cm = max(-1, min(Cm, 1));
end

function [position, isterminal, direction] = ground_impact(t, y)
    position = y(2);  % Height
    isterminal = 1;   % Stop integration when ground is reached
    direction = -1;   % Negative direction (approaching ground)
end

Answer 1

[Ashok]

Hi Rakshit,

In the code shared, I am assuming that 'theta' is the pitch angle of the projectile,'alpha' is the angle of attack, the positive angles indicate pitch up and negative angles indicate pitch down. Based on this, you should check the correctness of the following relationships:

1. The initial velocity components in the ‘aircraft_launch_trajectory_6DOF_advanced’ function.

    v0x = v0_aircraft - v0_launch * cosd(angle);
    v0y = -v0_launch * sind(angle);

2. Relations in the ‘equations_of_motion_6DOF_advanced’ function for 'alpha', total moment, and acceleration.

    alpha = atan2(vy, vx) - theta;
    M_total = M_aero + Fl * d_cp_cm;
    ax = (-Fd * cos(theta) + Fl * sin(theta) + boost_fx) / m;
    ay = -9.81 + (-Fd * sin(theta) - Fl * cos(theta) + boost_fy) / m;

Additionally, the drag coefficient calculated in the ‘aero_coefficients’ function and the characteristic length in the moment relation, ‘M_aero = q * Cm * A * length;’ seem quite high.

Since the ‘aero_coefficients’ function doesn’t account for the stall region, it might be better to use a smaller launch angle. Instead of ‘angle = -45;’ consider using a launch angle in the range -15 to 15 degrees.

To improve the visualization of the projectile’s trajectory, you can make the X and Y axes have the same scale by using the command ‘axis equal’.

You can find more details about this command in the link below:

https://www.mathworks.com/help/matlab/ref/axis.html#buk989s-1-style

I believe this will assist you!

[Rakshit Patil]

Thank you very much.