How to plot horizontal vs vertical amplitude at any node of a system

close all;
clear all;
R = 0.015 / 2;       % Radius of the shaft 
rhoe = 7860;         % Density of the shaft material (kg/m^3)
d = 0.015;           % Shaft diameter (m)
I_ps=(d^4*pi)/32;    % polar moment of inertia of the shaft
Kxx = 5 * 10^7;      % Define stiffness values
Kyy = 5 * 10^7; 
Cxx = 2 * 10^2;      % Define damping coefficient value
Cyy = 2 * 10^2; 
Ae = (pi*d^2)/4;     % Cross-sectional area (m^2) of shaft
L = 0.75;            % Length of shaft
d1=0.138;            % dia of disk
le = L / 24;         % Length of each element (m)
Ie = (pi*d^4)/64;    % Moment of inertia of the shaft cross-section (m^4)
md = 0.610;          %  mass of the disk in kg
mke=7.308*10^-4;
E = 2.1*10^11;       % Youngs modulus
             
I_p = (pi*(d1^4-0))/32;  % Polar moment of inertia of the disk
I_d=I_p/2;
mu = 0;            % crack                 
gima = sqrt(mu*sqrt(2 - mu));                                                                          
         
omega=205;
theta = 0; 
n_nodes = 25;        % Number of nodes
dofs_per_node = 4;   % DOFs per node (assuming 4 DOFs per node: displacements and rotations)
total_dofs = n_nodes * dofs_per_node; % Total number of degrees of freedom
%% with crack
Ixc = ((pi * R^4) / 8) - ((R^4 / 4) * (((1 - mu) * (((2 * mu^2) -( 4 * mu) + 1) * gima)) + asin(1 - mu)));                                                       
Iyc = (R^4 / 8) * (((1 - mu) * (((2 * mu^2) -( 4 * mu) - 3) * gima)) + (3 * asin( gima)));                                

% Transformed moments of area
Ix = (Ixc * cos(theta) * cos(theta)) + (Iyc * sin(theta) * sin(theta));                                    
Iy = (Iyc * cos(theta) * cos(theta) )+ (Ixc * sin(theta) * sin(theta));                                                              
Ixy = (Iyc - Ixc) * sin(theta) * cos(theta);

% Factor ft
ft = 0.5 + ((13 / 18) * cos(theta)) - ((2 / 9) * cos(theta) * cos(theta) * cos(theta));

% Cracked Stiffness matrix calculation
K_crack = (E / le^3) * [
    12*Ix, -Ixy, 6*le*Ixy, 6*le*Ix, -12*Ix, Ixy, 6 * le * Ixy, 6 * le * Ix;
    -Ixy, 12*Iy, -6 * le * Iy, -6 * le * Ixy, Ixy, -12*Iy, -6 * le * Iy, -6 * le * Ixy;
    6 * le * Ixy, -6 * le * Iy, 4 * le^2 * Iy, 4 * le^2 * Ixy, -6 * le * Ixy, 6 * le * Iy, 2 * le^2 * Iy, 2 * le^2 * Ixy;
    6 * le * Ix, -6 * le * Ixy, 4 * le^2 * Ixy, 4 * le^2 * Ix, -6 * le * Ix, 6 * le * Ixy, 2 * le^2 * Ixy, 2 * le^2 * Ix;
    -12*Ix, Ixy, -6 * le * Ixy, -6 * le * Ix, 12*Ix, -Ixy, -6 * le * Ixy, -6 * le * Ix;
    Ixy, -12*Iy, 6 * le * Iy, 6 * le * Ixy, -Ixy, 12*Iy, 6 * le * Iy, 6 * le * Ixy;
    6 * le * Ixy, -6 * le * Iy, 2 * le^2 * Iy, 2 * le^2 * Ixy, -6 * le * Ixy, 6 * le * Iy, 4 * le^2 * Iy, 4 * le^2 * Ixy;
    6 * le * Ix, -6 * le * Ixy, 2 * le^2 * Ixy, 2 * le^2 * Ix, -6 * le * Ix, 6 * le * Ixy, 4 * le^2 * Ixy, 4 * le^2 * Ix
];

% Scaling the stiffness matrix by ft
K_crs = ft * K_crack;

% Initialize the global stiffness matrix
K_global = zeros(total_dofs, total_dofs);

% Initialize the 100x100 matrix with zeros
K_global_stiffness_matrix = zeros(100, 100);

% Determine the starting index for element 8
start_idx = (7 * 4) + 1;

% Place the 8x8 element matrix in the global matrix
K_global_stiffness_matrix(start_idx:start_idx+7, start_idx:start_idx+7) = K_crs;

% Display the global stiffness matrix
disp(K_global_stiffness_matrix);

%% Total shaft stiffness matrix with out crack 

% Initialize global stiffness matrix
K_g = zeros(total_dofs, total_dofs);

% Define the local stiffness matrix for an element
K_e = ((E * Ie) / le^3) * [12, 0, 0, 6*le, -12, 0, 0, 6*le;
                           0, 12, -6*le, 0, 0, -12, -6*le, 0;
                           0, -6*le, 4*le^2, 0, 0, 6*le, 2*le^2, 0;
                           6*le, 0, 0, 4*le^2, -6*le, 0, 0, 2*le^2;
                           -12, 0, 0, -6*le, 12, 0, 0, -6*le;
                           0, -12, 6*le, 0, 0, 12, 6*le, 0;
                           0, -6*le, 2*le^2, 0, 0, 6*le, 4*le^2, 0;
                           6*le, 0, 0, 2*le^2, -6*le, 0, 0, 4*le^2];

for i = 1:(n_nodes - 1)
    start_idx = (i - 1) * dofs_per_node + 1;
    end_idx = start_idx + dofs_per_node * 2 - 1;
    K_g(start_idx:end_idx, start_idx:end_idx) = K_g(start_idx:end_idx, start_idx:end_idx) + K_e;
end

% Display the global stiffness matrix
disp(K_g);
K_g1 = zeros(total_dofs, total_dofs);
% Define the stiffness matrix at node 1
Kb1 = [Kxx, 0, 0, 0;
       0, Kyy, 0, 0;
       0, 0, 0, 0;
       0, 0, 0, 0];

% Assign Kb1 to the global stiffness matrix at node 1
start_idx = 1;
end_idx = start_idx + dofs_per_node - 1;
K_g1(start_idx:end_idx, start_idx:end_idx) = Kb1;

% at node 7
% Initialize global stiffness matrix
K_g2 = zeros(total_dofs, total_dofs);

% Define the stiffness matrix at node 25
Kb2 = [Kxx, 0, 0, 0;
       0, Kyy, 0, 0;
       0, 0, 0, 0;
       0, 0, 0, 0];

% Calculate start and end indices for node 7
start_idx = (n_nodes - 1) * dofs_per_node + 1;
end_idx = start_idx + dofs_per_node - 1;

% Assign Kb1 to the global stiffness matrix at node 25
K_g2(start_idx:end_idx, start_idx:end_idx) = Kb2;

K_total=K_g+K_g1+K_g2;

%% Mass matrix
M_g1 = zeros(total_dofs, total_dofs); % Translation mass matrix with shear


M_e_t = ((rhoe * Ae * le) / 420) * [
    156, 0, 0, 22 * le, 54, 0, 0, -13 * le;
    0, 156, -22 * le, 0, 0, 54, 13 * le, 0;
    0, -22 * le, 4 * le^2, 0, 0, -13 * le, -3 * le^2, 0;
    22 * le, 0, 0, 4 * le^2, 13 * le, 0, 0, -3 * le^2;
    54, 0, 0, 13 * le, 156, 0, 0, -22 * le;
    0, 54, -13 * le, 0, 0, 156, 22 * le, 0;
    0, 13 * le, -3 * le^2, 0, 0,22 * le,4 * le^2,0;
    -13 * le, 0, 0, -3 * le^2, -22 * le, 0, 0, 4 * le^2];

% Loop through each element and assemble the global translation mass matrix (M_g1)
for i = 1:(n_nodes - 1)
    start_idx = (i - 1) * dofs_per_node + 1;
    end_idx = start_idx + dofs_per_node * 2 - 1;
    M_g1(start_idx:end_idx, start_idx:end_idx) = M_g1(start_idx:end_idx, start_idx:end_idx) + M_e_t;
end

M_g2 = zeros(total_dofs, total_dofs); % Rotary mass matrix with shear

M_e_r = ((rhoe * Ie) / (le * 30)) * [
    36, 0, 0, 3 * le, -36, 0, 0, 3 * le;
    0, 36, -3 * le, 0, 0, -36, -3 * le, 0;
    0, -3 * le, 4 * le^2, 0, 0, 3 * le, -le^2, 0;
    3 * le, 0, 0, 4 * le^2, -3 * le, 0, 0, -le^2;
    -36, 0, 0, -3 * le, 36, 0, 0, -3 * le;
    0, -36, 3 * le, 0, 0, 36, 3 * le, 0;
    0, -3 * le, -le^2, 0, 0, 3 * le, 4 * le^2, 0;
    3 * le, 0, 0, -le^2, -3 * le, 0, 0, 4 * le^2
];
% Loop through each element and assemble the global translation mass matrix (M_g2)
for i = 1:(n_nodes - 1)
    start_idx = (i - 1) * dofs_per_node + 1;
    end_idx = start_idx + dofs_per_node * 2 - 1;
    M_g2(start_idx:end_idx, start_idx:end_idx) = M_g2(start_idx:end_idx, start_idx:end_idx) + M_e_r;
end
% Initialize the disk mass matrix
disk_mass_matrix_1 = zeros(total_dofs, total_dofs);

% Add the disk mass matrix to the global mass matrix (M_g3)
mid_node = 8 * dofs_per_node + 1;
disk_mass_matrix_1(mid_node:mid_node + dofs_per_node - 1, mid_node:mid_node + dofs_per_node - 1) = ...
    md * [
        1, 0, 0, 0;
        0, 1, 0, 0;
        0, 0, 0, 0;
        0, 0, 0, 0
    ] + I_d * [
        0, 0, 0, 0;
        0, 0, 0, 0;
        0, 0, 1, 0;
        0, 0, 0, 1
    ];
disp(disk_mass_matrix_1);

disk_mass_matrix_2 = zeros(total_dofs, total_dofs);

% Add the disk mass matrix to the global mass matrix (M_g3)
mid_node = 16 * dofs_per_node + 1;
disk_mass_matrix_2(mid_node:mid_node + dofs_per_node - 1, mid_node:mid_node + dofs_per_node - 1) = ...
    md * [
        1, 0, 0, 0;
        0, 1, 0, 0;
        0, 0, 0, 0;
        0, 0, 0, 0
    ] + I_d * [
        0, 0, 0, 0;
        0, 0, 0, 0;
        0, 0, 1, 0;
        0, 0, 0, 1
    ];
disp(disk_mass_matrix_2);
M_global=M_g1+M_g2+disk_mass_matrix_1+disk_mass_matrix_2;
Y=-inv(M_global);

%% Gyroscopoic matrix
M_g = zeros(total_dofs, total_dofs); 

% Element mass matrix (gyroscopic)
M_z =(((rhoe * I_ps*omega)) / (15 * le)) * [
    0, 36, 3 * le, 0, 0, 36, -3 * le, 0;
    -36, 0, 0, -3 * le, 36, 0, 0, -3 * le;
    3 * le, 0,0, 4 * le^2, -3 * le,  0, 0, -le^2;
    0, 3 * le, -4 * le^2, 0, 0, -3 * le, le^2, 0;
    0, -36, 3 * le, 0, 0, 36, 3 * le, 0;
    36, 0, 0, 3 * le, -36, 0, 0, 3 * le;
    3 * le, 0, 0, -le^2, -3 * le, 0, 0, 4 * le^2;
    0, 3 * le, le^2, 0, 0, -3 * le, -4 * le^2, 0
];
for i = 1:(n_nodes - 1)
    start_idx = (i - 1) * dofs_per_node + 1;
    end_idx = start_idx + dofs_per_node * 2 - 1;
    M_g(start_idx:end_idx, start_idx:end_idx) = M_g(start_idx:end_idx, start_idx:end_idx) + M_z;
end
disk_gyroscope_matrix_1 = zeros(total_dofs, total_dofs);

% Add the disk mass matrix_1 to the global mass matrix (M_g3)
mid_node = 8 * dofs_per_node + 1;
disk_gyroscope_matrix_1(mid_node:mid_node + dofs_per_node - 1, mid_node:mid_node + dofs_per_node - 1) = ...
    I_p * [
        0, 0, 0, 0;
        0, 0, 0, 0;
        0, 0, 0, 1;
        0, 0, -1, 0
    ];

disk_gyroscope_matrix_2 = zeros(total_dofs, total_dofs);

% Add the disk mass matrix_2 to the global mass matrix (M_g3)
mid_node = 16 * dofs_per_node + 1;
disk_gyroscope_matrix_2(mid_node:mid_node + dofs_per_node - 1, mid_node:mid_node + dofs_per_node - 1) = ...
    I_p * [
        0, 0, 0, 0;
        0, 0, 0, 0;
        0, 0, 0, 1;
        0, 0, -1, 0
    ];
G_global=M_g+disk_gyroscope_matrix_1+disk_gyroscope_matrix_2;

C_g1 = zeros(total_dofs, total_dofs);
% Define the stiffness matrix at node 1
Cb1 = [Cxx, 0, 0, 0;
       0, Cyy, 0, 0;
       0, 0, 0, 0;
       0, 0, 0, 0];

% Assign Kb1 to the global stiffness matrix at node 1
start_idx = 1;
end_idx = start_idx + dofs_per_node - 1;
C_g1(start_idx:end_idx, start_idx:end_idx) = Cb1;

% at node 7
% Initialize global stiffness matrix
C_g2 = zeros(total_dofs, total_dofs);

% Define the stiffness matrix at node 25
Cb2 = [Cxx, 0, 0, 0;
       0, Cyy, 0, 0;
       0, 0, 0, 0;
       0, 0, 0, 0];

% Calculate start and end indices for node 7
start_idx = (n_nodes - 1) * dofs_per_node + 1;
end_idx = start_idx + dofs_per_node - 1;

% Assign Kb1 to the global stiffness matrix at node 25
C_g2(start_idx:end_idx, start_idx:end_idx) = Cb2;
C_total=C_g1+C_g2;

%% force
% Initialize the overall force matrix
F_overall_9 = zeros(100, 1);
F_overall_17 = zeros(100, 1);

% Given 4x1 force matrix for node 9
F_node_9 = [mke*omega^2*cos(theta+0);mke*omega^2*sin(theta+0); 0; 0];

F_node_17 = [mke*omega^2*cos(theta+0);mke*omega^2*sin(theta+0); 0; 0];

% Node 9 corresponds to positions 33 to 36
F_overall_9(33:36) = F_node_9;
% Node 17 corresponds to positions 65 to 68
F_overall_17(65:68) = F_node_17;
F_total=F_overall_9+F_overall_17;



% Define initial conditions
X0 = zeros(100, 1); % Initial displacement
dX0 = zeros(100, 1); % Initial velocity

% Time span
tspan = 0:1:100; 
A = [zeros(100) eye(100); -M_global\K_global -M_global\(C_total+(omega*G_global))];
B = [zeros(100,100); M_global\eye(100)];

% Define the ODE function including external force
odefun = @(t, Y) A * Y + B * F_total;
% Solve the system of equations
Y0 = [X0; dX0]; % Initial state vector
[t, Y] = ode45(odefun, tspan, Y0);

% Extract the displacement components
X = Y(:, 1:100);

% Select node of interest
node = 9; % Example node (1 to 24)
dof_per_node = 4;
u_idx = (node-1)*dof_per_node + 1; % Horizontal displacement index
v_idx = (node-1)*dof_per_node + 2; % Vertical displacement index

u = X(:, u_idx); % Horizontal displacement
v = X(:, v_idx); % Vertical displacement

% Plot horizontal amplitude vs vertical amplitude
figure;
plot(t, u);
xlabel('Horizontal Amplitude (u)');
ylabel('Vertical Amplitude (v)');
title(['Horizontal vs Vertical Amplitude at Node ' num2str(node)]);
grid on;
 

i have to plot vertical vs horizontal amplitude for any node using ode45. for a system having 24 element and 25 node and each node having 4 degree of freedom two transalational and two rotational. and two unbalance disk are placed at node 9 and node 17

Answer 1

[Umar]

Hi @ prabin kumar meher, very impressive job with your code providing insights into the behavior of the structural system under dynamic loading conditions, showcasing how displacements evolve over time at a specific node. However, to address your comments regarding, “plot the vertical vs horizontal amplitude for any node in a system with 24 elements, 25 nodes, and each node having 4 degrees of freedom (2 translational and 2 rotational) with two unbalanced disks placed at node 9 and node 17 using `ode45`”, Define the system parameters and matrices as provided in your content then set the initial conditions for displacement and velocity. Afterwards, define the ODE function that includes external force based on the system matrices and use `ode45` to solve the system of equations over a specified time span. Finally, extract the displacement components from the solution. If you want to plot to analyze your results, select a specific node of interest (e.g., Node 9),obtain the horizontal and vertical displacements for that node and plot the horizontal amplitude vs vertical amplitude. I will illustrate it with example,

% Set initial conditions

X0 = zeros(100, 1); % Initial displacement

dX0 = zeros(100, 1); % Initial velocity

% Time span

tspan = 0:1:100;

% Define system matrices A and B

% Define the ODE function with external force

odefun = @(t, Y) A * Y + B * F_total;

% Solve ODEs

Y0 = [X0; dX0]; % Initial state vector

[t, Y] = ode45(odefun, tspan, Y0);

% Extract displacement components

X = Y(:, 1:100);

% Select node of interest (e.g., Node 9)

node = 9;

dof_per_node = 4;

u_idx = (node-1)*dof_per_node + 1; % Horizontal displacement index

v_idx = (node-1)*dof_per_node + 2; % Vertical displacement index

u = X(:, u_idx); % Horizontal displacement

v = X(:, v_idx); % Vertical displacement

% Plot horizontal amplitude vs vertical amplitude

figure;

plot(u, v);

xlabel('Horizontal Amplitude (u)');

ylabel('Vertical Amplitude (v)');

title(['Horizontal vs Vertical Amplitude at Node ' num2str(node)]);

grid on;

Hope this helps.

[prabin kumar meher]

thanks umar