-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathbernoullieuler_beamvibration_withoutcontrol.m
133 lines (112 loc) · 4.85 KB
/
bernoullieuler_beamvibration_withoutcontrol.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
% Bernoulli-Euler Beam Vibration with FDM and Adjusted Robust Boundary Control for Stability
% Parameters
L = 41.4*10^-2; % Beam length in meters
wd = 14.88*10^-3; % Beam width in meters
t = 1*10^-3; % Beam thickness in meters
E = 70*10^9; % Modulus of elasticity in N/m^2
ro = 5186.5; % Linear mass density in kg/m^2
I = ((wd*(t^3))/12); % Moment of inertia in kgm^2
A = wd*t; % Cross-sectional area in m^2
lamda = A*ro;
ksi = 0.2; % Damping ratio
P = 0.05; % Tip payload in N
nt = 2*10^5; % Length of time domain
nx = 15; % Length of space domain
tf = 1; % Time of simulation in seconds
g = 9.81;
% Adjusted Control Gains for Stability
%k1 = 0.01; % Reduced proportional control gain
%k2 = 0.02; % Reduced derivative control gain
%k = 100; % Reduced feedback gain
k1=0;
k2=0;
k=0;
% Compute the mesh spacing and time step
dx = L/nx; % Space step
dt = tf/nt; % Time step
% Create mesh for drawing
x = linspace(0, L, nx+1); % Space mesh
time = linspace(0, tf, nt+1); % Time mesh
% Preallocate matrices for efficiency
w = zeros(nx+1, nt+1); % Displacement
d = zeros(nt+1, 1); % Boundary disturbance
f = zeros(nx+1, nt+1); % External disturbance
% Natural frequencies and damping coefficient
k_vals = [1.875104069, 4.694091133, 7.854757438, 10.99554073, 14.13716839, 17.27875953] / L;
w_n = sqrt((k_vals.^4) * E * I / lamda);
ra = (2 * ksi * w_n(1) * ro * A);
% Boundary and External disturbance definitions
for j = 1:nt+1
d(j) = 0.1 + 0.1*sin(pi*(j-1)*dt) + 0.1*sin(2*pi*(j-1)*dt) + 0.1*sin(3*pi*(j-1)*dt);
end
for i = 1:nx+1
for j = 1:nt+1
f(i,j) = (1 + sin(0.1*pi*(i-1)*dx*(j-1)*dt) + sin(0.2*pi*(i-1)*dx*(j-1)*dt) + sin(0.3*pi*(i-1)*dx*(j-1)*dt)) * (i-1)*dx / (20*L);
end
end
% Initial Conditions
scaling_factor = 0.5; % Adjust this factor to scale down the initial displacement
for i = 1:nx+1
initial_deflection = -(P*L/(E*I)) * (((x(i).^2)./2) - ((x(i).^3)./(6*L)));
w(i,1) = scaling_factor * initial_deflection;
w(i,2) = w(i,1); % Assuming zero initial velocity
end
% Damping coefficient for time integration
%alpha = (dt * ra / (ro * A));
alpha=0;
% Constants for boundary conditions
%T = 2; % Tension
T=0;
% Control point
%control_point = round(nx / 2); % Middle of the beam
control_point=3;
% Main loop for finite difference method
for j = 3:nt+1
for i = 3:nx-1
wxxxx = (w(i+2,j-1) - 4*w(i+1,j-1) + 6*w(i,j-1) - 4*w(i-1,j-1) + w(i-2,j-1)) / dx^4;
w(i,j) = ((2 + alpha) * w(i,j-1) - w(i,j-2) + (-E*I*wxxxx + f(i,j)) * (dt^2 / (ro * A))) / (1 + alpha);
end
w(nx+1,j) = 2*w(nx,j) - w(nx-1,j);
% Apply control law at the control point
%wx_control = (w(control_point+1,j-1) - w(control_point,j-1)) / dx;
%wxxx_control = (w(control_point+1,j-1) - 3*w(control_point,j-1) + 3*w(control_point-1,j-1) - w(control_point-2,j-1)) / dx^3;
%wt_control = (w(control_point+1,j-1) - w(control_point+1,j-2)) / dt;
%wtx_control = (wt_control - (w(control_point,j-1) - w(control_point,j-2)) / dt) / dx;
%wtxxx_control = (wtx_control - ((w(control_point-1,j-1) - w(control_point-1,j-2)) / dt - (w(control_point-2,j-1) - w(control_point-2,j-2)) / dt) / dx^2) / dx;
%u = -E*I*wxxx_control + T*wx_control - (P/g) * (k1*wtx_control - k2*wtxxx_control) - k*(wt_control - k2*wxxx_control + k1*wx_control) - d(j-1);
%w(control_point,j) = w(control_point,j-1) + dt * u; % Control action applied at the control point
% Fixed end boundary conditions at the starting point
w(1,j) = 0; % Displacement at the fixed end is zero
w(2,j) = 0; % Slope at the fixed end is zero
end
% Plotting the displacement
figure(1)
mesh(x, time, w');
title('Displacement of Beam without Control');
ylabel('Time (s)', 'FontSize', 12);
xlabel('x (m)', 'FontSize', 12);
zlabel('w(x,t) (m)', 'FontSize', 12);
view([60 45]);
% Set up the figure for animation
h = figure;
axis tight manual; % Ensure that getframe() returns a consistent size
filename = 'beam_vibration_without_control_animation.avi'; % Name of the video file
v = VideoWriter(filename); % Create a VideoWriter object
open(v);
frame_skip = 500; % Capture every 500th frame for the animation
short_pause = 0.1; % Short pause time, adjust as needed
% Animation of displacement with frame skipping and reduced pause time
for j = 1:frame_skip:length(time) % Loop with frame skipping
plot(x, w(:,j));
xlabel('x (m)');
ylabel('Displacement (m)');
title(sprintf('Beam Displacement at Time = %.2f s', time(j)));
axis([0 L min(w(:)) max(w(:))]);
drawnow;
pause(short_pause); % Short pause for animation speed control
% Capture the plot as a frame and write to the video
frame = getframe(h); % Capture the current plot
writeVideo(v, frame);
end
% Close the video file
close(v);