Periodic result using ode45 while simulating a parachute decent with 6DOF dynamics

So I've been trying to simulate a parachute descent (very basic model, almost a mass body with a coefficient of drag), from a height. While simulating only in the axial direction(z axis), the results coincide with kinematics equations. However, when I give a non-zero initial condition, the outputs go crazy (definitely not correct) with periodic solutions to related terms. I seem to have checked my equations, they seem to be right. Is there anything with 'odeset' I can work this out around.

Zero non-axial initial condition:

x0 = [0; 0; -100; 0; 0; 1; 0; 0; 0; 0; 0; 0];

Non-zero IC:

x0 = [0; 0; -100; 0; 1; 1; 0; 0; 0; 0; 0; 0];

0

this is the code i've used for these plots. Please feel free to ask any clarifications

%% 6DOF fucntion definition
function [dxdt, aoa_series] = six_DOF_dynamics(t,x)


M_para = 5; %mass of the parachute in Kgs
g = 9.81; %gravitational acceleration in m/s^2
I_x = 1; I_y = 1; I_z = 3; %Inertia matrix diagonal principal of the parachute in kg-m^2
CD_para = 0.9; %coefficient of drag of parachute
CL_para = 0; %coefficient of lift of parachute
rho = 1; %density of the atmosphere in Kg/m^3
x_com = 0; y_com = 0; z_com = -5; %position of COM of the paarchute from the integrated system's COM(body axis origin) in m
A = 10; %parachute area in m^2

    %state variables defining
    x_pos = x(1); y_pos = x(2); z_pos = x(3); %postion variables
    u = x(4); v = x(5); w = x(6); %velocity variables 
    phi = x(7); theta = x(8); psi = x(9); %attitude/euler angles variables
    p = x(10); q = x(11); r = x(12); %euler angular rate variables

    %rate of change of postion -> velocity
    dx_pos = u;
    dy_pos = v;
    dz_pos = w;

    %Gravity moment of the parachute along x,y and z components
    GravityMoment_para_x = M_para*g*cos(theta)*(y_com*cos(phi) - z_com*sin(phi));
    GravityMoment_para_y = M_para*g*(z_com*sin(theta) - x_com*cos(theta)*cos(phi));
    GravityMoment_para_z = M_para*g*(x_com*cos(theta)*sin(phi) + y_com*sin(theta));

    V_total = sqrt(u^2 + v^2 + w^2); %net speed
    
    persistent angle_of_attack %defining this variable to take it as a output outside ode45

    alpha = acos(w/V_total); %angle of attack calculated in radians

    angle_of_attack = [angle_of_attack; alpha]; %concatinating current iterated angle of attack in global variable

    %Componetize CD and CL into tangential and normal force coefficients
    CT = CD_para*cos(alpha) - CL_para*sin(alpha);
    CN = CL_para*cos(alpha) + CD_para*sin(alpha);

    % Aerodynamic force along x,y,z axes F = 0.5*rho*v^2*A
    if (u == 0 & v == 0) 
        Aerodynamic_force_x = 0;
        Aerodynamic_force_y = 0;
        Aerodynamic_force_z = -0.5*rho*CT*V_total*V_total*A;
    else 
        Aerodynamic_force_x = -0.5*rho*CN*V_total*V_total*A*(u/(u^2 + v^2)^0.5); 
        Aerodynamic_force_y = -0.5*rho*CN*V_total*V_total*A*(v/(u^2 + v^2)^0.5);
        Aerodynamic_force_z = -0.5*rho*CT*V_total*V_total*A;
    end

 
    %Moment caused by aerodynamic force along x,y,x
    Aerodynamic_moment_x = y_com*Aerodynamic_force_z - z_com*Aerodynamic_force_y;
    Aerodynamic_moment_y = z_com*Aerodynamic_force_x - x_com*Aerodynamic_force_z;
    Aerodynamic_moment_z = x_com*Aerodynamic_force_y - y_com*Aerodynamic_force_x;

    %rate change of velocity definition
    du = v*r - w*q + g*sin(theta) + Aerodynamic_force_x/M_para;
    dv = -u*r + w*p + g*cos(theta)*sin(phi) + Aerodynamic_force_y/M_para;
    dw = u*q - v*p + g*cos(theta)*cos(phi) + Aerodynamic_force_z/M_para;

    %rate change of euler angles definition
    dphi = p + (r*cos(phi) + q*sin(phi))*tan(theta);
    dtheta = q*cos(phi) - r*sin(phi);
    dpsi = (cos(theta)^-1)*(r*cos(phi) + q*sin(phi));

    %rate change of angular rates definition
    dp = (1/I_x)*(Aerodynamic_moment_x + GravityMoment_para_x);
    dq = (1/I_y)*(Aerodynamic_moment_y + GravityMoment_para_y);
    dr = (1/I_z)*(Aerodynamic_moment_z + GravityMoment_para_z);

    dxdt = [dx_pos; dy_pos; dz_pos; du; dv; dw; dphi; dtheta; dpsi; dp; dq; dr];

    %block to output angle of attack
    if nargout>1
        aoa_series = angle_of_attack;

    end
      
end

tspan = [0 3];

x0 = [0; 0; -100; 0; 1; 1; 0; 0; 0; 0; 0; 0];

[t, x] = ode45(@six_DOF_dynamics, tspan, x0);

[~, aoa_series] = six_DOF_dynamics(t,x);

%adding alpha as a part of x


% Plot results
figure;
subplot(4,1,1);
plot(t, x(:,1:3));
title('Position');
legend('x', 'y', 'z');

subplot(4,1,2);
plot(t, x(:,4:6));
title('Velocity');
legend('u', 'v', 'w');

subplot(4,1,3);
plot(t, x(:,7:9));
title('Euler Angles');
legend('phi', 'theta', 'psi');

subplot(4,1,4);
plot(t, x(:,10:12));
title('Euler Rates');
legend('p', 'q', 'r');

[Sam Chak]

Hi, Are you looking for a Mathematical Physicist to review equations in your code? If so, can you provide the derived mathematical model?

[Paul]

I didn't look too closely at the model, but I did notice that the equations for dx_pos, dy_pos, and dz_pos appear to be incorrect. I suspect that x_pos, y_pos, and z_pos are supposed to be positions of the mass in an inertial frame, in which case their derivatives need to be components of the velocity vector resolved in the inertial frame. But u, v, and w, are components of the velocity vector resolved in the body frame. I don't think this issue explains the results in question becaue x_pos, y_pos, and z_pos aren't used for anything but plotting, but it should be addressed.

Assuming all angles and angle rates are in rad and rad/sec (as suggested by the code), it looks like the periodicity in v and w is driven by the constant rate of change in roll angle. Investigation of the dp and dphi equations may be fruitful.

Also, p, q, and r are not Euler angle rates.

[mekg_10]

@Sam Chak its from an article named "A 6DOF mathematical model of parachute in Mars EDL". Link to article

[Sam Chak]

Hi @mekg_10

Thanks for sharing the paper.

I am not an expert in Mars EDL (Entry, Descent, and Landing); however, it is clear that in the simulation, the Mars EDL lacks any control variable, as the control variable cannot be identified in the numbered equations of the vehicle dynamics in the paper. In the paper, the authors claim to have used a conventional PID controller and Prof. Han's Active Disturbance Rejection Control (ADRC), although the specific control laws are not published in the 2015 paper.

Edit:

I'm NOT good at reading too many coupled equations with long variable names. But I have also included annotations in the comments to indicate some missing equations, and incorrect definition of Mpy.

%% 6-DOF function definition
function [dxdt, aoa_series] = six_DOF_dynamics(t,x)
    M_para  =  5;       % mass of the parachute in Kgs
    g       =  9.81;    % gravitational acceleration in m/s^2
    I_x     =  1; 
    I_y     =  1; 
    I_z     =  3;       % Inertia matrix diagonal principal of the parachute in kg-m^2
    CD_para =  0.9;     % coefficient of drag of parachute
    CL_para =  0;       % coefficient of lift of parachute
    rho     =  1;       % density of the atmosphere in Kg/m^3
    x_com   =  0; 
    y_com   =  0; 
    z_com   = -5;       % position of COM of the paarchute from the integrated system's COM(body axis origin) in m
    A       = 10;       % parachute area in m^2

    %state variables defining
    x_pos   = x(1); 
    y_pos   = x(2); 
    z_pos   = x(3);     % postion variables
    u       = x(4); 
    v       = x(5); 
    w       = x(6);     % velocity variables 
    phi     = x(7); 
    theta   = x(8); 
    psi     = x(9);     % attitude/euler angles variables
    p       = x(10); 
    q       = x(11); 
    r       = x(12);    % euler angular rate variables

    % rate of change of postion -> velocity
    dx_pos  = u;
    dy_pos  = v;
    dz_pos  = w;

    %% See first three sub-equations in Equation (7), but last three {Mvx, Mvy, Mvz} are unused
    GravityMoment_para_x = M_para*g*cos(theta)*(  y_com*cos(phi)            - z_com*sin(phi));              % Mpx
    GravityMoment_para_y = M_para*g*           (- z_com*sin(theta)          - x_com*cos(theta)*cos(phi));   % Mpy (originally incorrect)
    GravityMoment_para_z = M_para*g*           (  x_com*cos(theta)*sin(phi) + y_com*sin(theta));            % Mpz

    V_total = sqrt(u^2 + v^2 + w^2); % net speed
    
    persistent angle_of_attack %defining this variable to take it as a output outside ode45

    alpha   = acos(w/V_total); %angle of attack calculated in radians

    angle_of_attack = [angle_of_attack; alpha]; %concatinating current iterated angle of attack in global variable

    % Componetize CD and CL into tangential and normal force coefficients
    CT      = CD_para*cos(alpha) - CL_para*sin(alpha);
    CN      = CL_para*cos(alpha) + CD_para*sin(alpha);

    % Aerodynamic force along x,y,z axes F = 0.5*rho*v^2*A
    if (u == 0 & v == 0) 
        Aerodynamic_force_x =  0;
        Aerodynamic_force_y =  0;
        Aerodynamic_force_z = -0.5*rho*CT*V_total*V_total*A;
    else 
        %% See first three sub-equations in Equation (8)
        Aerodynamic_force_x = -0.5*rho*CN*V_total*V_total*A*(u/(u^2 + v^2)^0.5); 
        Aerodynamic_force_y = -0.5*rho*CN*V_total*V_total*A*(v/(u^2 + v^2)^0.5);
        Aerodynamic_force_z = -0.5*rho*CT*V_total*V_total*A;
    end

    %% See last three sub-equations in Equation (8)
    Aerodynamic_moment_x = y_com*Aerodynamic_force_z - z_com*Aerodynamic_force_y;
    Aerodynamic_moment_y = z_com*Aerodynamic_force_x - x_com*Aerodynamic_force_z;
    Aerodynamic_moment_z = x_com*Aerodynamic_force_y - y_com*Aerodynamic_force_x;

    %% See Equation (9)
    du      =   v*r - w*q + g*sin(theta)          + Aerodynamic_force_x/M_para;
    dv      = - u*r + w*p + g*cos(theta)*sin(phi) + Aerodynamic_force_y/M_para;
    dw      =   u*q - v*p + g*cos(theta)*cos(phi) + Aerodynamic_force_z/M_para;

    %% See Equation (12)
    dphi    = p + (r*cos(phi) + q*sin(phi))*tan(theta);
    dtheta  =      q*cos(phi) - r*sin(phi);
    dpsi    =     (r*cos(phi) + q*sin(phi))/cos(theta);

    %% See Equation (10), but Mvx, Mvy, Mvz from Eq. (7) are neglected
    dp      = (1/I_x)*(Aerodynamic_moment_x + GravityMoment_para_x);
    dq      = (1/I_y)*(Aerodynamic_moment_y + GravityMoment_para_y);
    dr      = (1/I_z)*(Aerodynamic_moment_z + GravityMoment_para_z);
    
    %% ODE vector
    dxdt    = [dx_pos; 
               dy_pos; 
               dz_pos; 
               du; 
               dv; 
               dw; 
               dphi; 
               dtheta; 
               dpsi; 
               dp; 
               dq; 
               dr];

    % block to output angle of attack
    if nargout>1
        aoa_series = angle_of_attack;
    end
end

tspan   = [0 3];
x0      = [0; 0; -100; 0; 1; 1; 0; 0; 0; 0; 0; 0];
[t, x]  = ode45(@six_DOF_dynamics, tspan, x0);
[~, aoa_series] = six_DOF_dynamics(t,x);

% adding alpha as a part of x

% Plot results
figure;
plot(t, x(:,1:3)), grid on
title('Position');
legend('x', 'y', 'z');

figure;
plot(t, x(:,4:6)), grid on
title('Velocity');
legend('u', 'v', 'w');

figure;
plot(t, x(:,7:9)), grid on
title('Euler Angles');
legend('phi', 'theta', 'psi');

figure;
plot(t, x(:,10:12)), grid on
title('Euler Rates');
legend('p', 'q', 'r');

[mekg_10]

@Sam Chak Thanks for the response. I see that you've performed the simulation from the mathematical model. But I do notice you've encountered the same problem as mine. I feel like the simulation results are wrong. Why? Because when you drop a mass from a height with some initial velocity( along both z-axis and x/y-axis), from the simulation results velocity v & w, we see its periodic. I do get that it's coupled and the prescence of a non-zero v can affect w but will it be to an extent where it goes from +3.2 to -3.2 in a high frequency fashion...This what I feel makes the simulations wrong

[Sam Chak]

Because the control variable in the Mars EDL dynamics is missing!

[mekg_10]

@Paul Thank you very much for pointing my mistake. I've used the tranformation matrix now but I seem to have another error...the frequency of the results are much lower now lol, but still periodic.... Yes you are right. I interchangebly plotted variables because some give a better idea in a frame of reference.

Let me look into the dp and dphi equations like you suggested. By any chance your background is Aerospace engineering?

[mekg_10]

@Sam Chak but I don't get how the missing control variable explains the weird phenomenon of the velocity of oscillating for a falling body, which practically cannot be right. From what I know, it should increase in magnitude and reach a terminal velocity

Answer 1

[Sam Chak]

Hi @mekg_10

At first, I thought you were trying to duplicate the results of the paper by attaching a parachute to the Mars EDL. However, after looking at the simulation parameters, it seems that you’re actually conducting a parachute simulation study on Earth. I’m not sure if the following simulation is correct or not from a purely mathematical perspective.

%% Parachute descent dynamics
function [dstate, aoa_series] = Parachute_dynamics(t, state)
    %% parameters
    m       =  5;           % mass of the parachute in Kgs
    g       =  9.81;        % gravitational acceleration in m/s^2
    I_x     =  1; 
    I_y     =  1; 
    I_z     =  3;           % Inertia matrix diagonal principal of the parachute in kg-m^2
    CD_para =  0.9;         % coefficient of drag of parachute
    CL_para =  0;           % coefficient of lift of parachute
    rho     =  1;           % density of the atmosphere in Kg/m^3
    x_com   =  0; 
    y_com   =  0; 
    z_com   = -5;           % position of COM of the paarchute from the integrated system's COM(body axis origin) in m
    A       = 10;           % parachute area in m^2

    %% state variables defining
    x       = state(1); 
    y       = state(2); 
    z       =-state(3);     % postion variables
    u       = state(4); 
    v       = state(5); 
    w       = state(6);     % velocity variables 
    phi     = state(7); 
    theta   = state(8); 
    psi     = state(9);     % attitude/euler angles variables
    p       = state(10); 
    q       = state(11); 
    r       = state(12);    % angular rate velocity components

    %% Moments due to gravitational acceleration
    Mg_x    = m*g*cos(theta)*(  y*cos(phi)            - z*sin(phi));              % Mpx
    Mg_y    = m*g*           (- z*sin(theta)          - x*cos(theta)*cos(phi));   % Mpy (originally incorrect)
    Mg_z    = m*g*           (  x*cos(theta)*sin(phi) + y*sin(theta));            % Mpz

    V_total = sqrt(u^2 + v^2 + w^2); % net speed
    
    persistent angle_of_attack %defining this variable to take it as a output outside ode45

    alpha   = acos(w/V_total); %angle of attack calculated in radians

    angle_of_attack = [angle_of_attack; alpha]; %concatinating current iterated angle of attack in global variable

    % Componetize CD and CL into tangential and normal force coefficients
    CT      = CD_para*cos(alpha) - CL_para*sin(alpha);
    CN      = CL_para*cos(alpha) + CD_para*sin(alpha);

    %% Aerodynamic forces along x, y, z axes, F = 0.5*rho*v^2*A
    if (u == 0 & v == 0) 
        Fa_x =  0;
        Fa_y =  0;
        Fa_z = -0.5*rho*CT*V_total*V_total*A;
    else 
        Fa_x = -0.5*rho*CN*V_total*V_total*A*(u/(u^2 + v^2)^0.5); 
        Fa_y = -0.5*rho*CN*V_total*V_total*A*(v/(u^2 + v^2)^0.5);
        Fa_z = -0.5*rho*CT*V_total*V_total*A;
    end

    %% Aerodynamic moments
    Ma_x    = y*Fa_z - z*Fa_y;
    Ma_y    = z*Fa_x - x*Fa_z;
    Ma_z    = x*Fa_y - y*Fa_x;

    %% Translational kinematics w.r.t. Earth
    dx_pos  = u*(cos(theta)*cos(psi)) + v*(sin(phi)*sin(theta)*cos(psi) - cos(phi)*sin(psi)) + w*(cos(phi)*sin(theta)*cos(psi) + sin(phi)*sin(psi));
    dy_pos  = u*(cos(theta)*sin(psi)) + v*(sin(phi)*sin(theta)*sin(psi) + cos(phi)*cos(psi)) + w*(cos(phi)*sin(theta)*sin(psi) - sin(phi)*cos(psi));
    dh_pos  = u* sin(theta)           - v* sin(phi)*cos(theta)                               - w* cos(phi)*cos(theta);
    
    %% Rotational kinematics
    du      =   v*r - w*q + g*sin(theta)          + Fa_x/m;
    dv      = - u*r + w*p + g*cos(theta)*sin(phi) + Fa_y/m;
    dw      =   u*q - v*p + g*cos(theta)*cos(phi) + Fa_z/m;
    
    %% Force equations
    dphi    = p + (r*cos(phi) + q*sin(phi))*tan(theta);
    dtheta  =      q*cos(phi) - r*sin(phi);
    dpsi    =     (r*cos(phi) + q*sin(phi))/cos(theta);

    %% Moment equations
    dp      = (1/I_x)*(Ma_x + Mg_x + (I_y - I_z)*q*r);
    dq      = (1/I_y)*(Ma_y + Mg_y + (I_z - I_x)*r*p);
    dr      = (1/I_z)*(Ma_z + Mg_z + (I_x - I_y)*p*q);
    
    %% ODE vector
    dstate  = [dx_pos; 
               dy_pos; 
               dh_pos; 
               du; 
               dv; 
               dw; 
               dphi; 
               dtheta; 
               dpsi; 
               dp; 
               dq; 
               dr];

    % block to output angle of attack
    if nargout>1
        aoa_series = angle_of_attack;
    end
end

%% Touchdown Event: The moment at which the parachute make contact with the ground during descending
function [position, isterminal, direction] = touchdown(~, state)
    position    = state(3); % The height that we want to be zero
    isterminal  =  1;       % Halt integration 
    direction   = -1;       % height is in decreasing direction
end

%% Call ode45 solver
tspan   = [0 40];
x0      = [0; 0; 100; 0; 1; 1; 0; 0; 0; 0; 0; 0];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-9, 'Events', @touchdown);
[t, x]  = ode45(@Parachute_dynamics, tspan, x0, options);

%% Plot result
figure;
plot(t, x(:,3)), grid on
xlabel('Time (t) / s'), ylabel('Height (h) / m')
title('Parachute height');