clear; clc;
%% Parameters
L = 1; % Beam length (m)
EI = 2.1*10^5; % Flexural rigidity (Pa.m^4)
% rhoA = 2.5*2450*10^(-4); % Mass per unit length (kg/m)
C = .5; % Stability constant
p=100;
omega=2*pi;
Nx = 10; % Number of spatial points
x = linspace(0, L, Nx); % Spatial discretization
dx = x(2) - x(1); % Spatial step size
k = 0; % Stiffness coefficient
c = 0; % Damping coefficient
m = 78; % rhoA Mass coefficient
% Calculate time step size for stability
dt = C * (dx^2) / sqrt(EI / m); % Time step for stability
T = 1.0; % Total time
Nt = floor(T/dt); % Number of time steps
t = linspace(0,T, Nt ); %t values
% Beam initial conditions
w= zeros(Nt, Nx); % Displacement over time
v = zeros(1,Nx); % Initial velocity
%% Initial deflection (example: sine wave)
w(1, :) = sin(pi*x/L); % Initial condition for beam deflection
w( 2,:) = w(1,:);
% Load function
f = zeros(Nt , Nx );
for n = 1:Nt
f(n, :) =p* cos(omega * t(n));
end
%% Finite Difference Method
for n = 2:Nt-1
for i = 2:Nx-1 % Spatial points excluding boundaries
% Boundary conditions
if(i == 1 || i == Nx) % Outside Boundary
w(n,i) = 0;
elseif(i==2) % Near left boundary
w(n+1,i) = ((-EI / dx^4) * (w(n,i+2) - 4*w(n,i+1) + 6*w(n,i) - 4*w(n,i-1) - w(n,i)) - ...
2*k*w(n,i) + ...
(c/dt)*w(n-1,i) + ...
(m / dt^2) * (2*w(n,i) - w(n-1,i)) + ...
p * cos(omega * t(n))) / ...
((m / dt^2) + (2*c / (2*dt)));
elseif i == Nx-1 % Near right boundary
w(n+1,i) = ((-EI / dx^4) * (-w(n,i) - 4*w(n,i+1) + 6*w(n,i) - 4*w(n,i-1) + w(n,i-2)) - ...
2*k*w(n,i) + ...
(c/dt)*w(n-1,i) + ...
(m / dt^2) * (2*w(n,i) - w(n-1,i)) + ...
p * cos(omega * t(n))) / ...
((m / dt^2) + (2*c / (2*dt)));
else % Interior points
w(n+1,i) = ((-EI / dx^4) * (w(n,i+2) - 4*w(n,i+1) + 6*w(n,i) - 4*w(n,i-1) + w(n,i-2)) - ...
2*k*w(n,i) + ...
(c/dt)*w(n-1,i) + ...
(m / dt^2) * (2*w(n,i) - w(n-1,i)) + ...
p * cos(omega * t(n))) / ...
((m / dt^2) + (2*c / (2*dt)));
end
end
end
%% Plot results
figure;
for n = 1:100:Nt
plot(x, w(n,:));
title(['Time = ', num2str(n*dt), ' sec']);
xlabel('x (m)'); ylabel('Deflection (m)');
grid on;
pause(0.1);
end
