Con lắc lò xo

Xem thêm: Chương trình mô phỏng chuyên nghiệp về con lắc lò xo viết bằng java.

Tại bài viết này, chúng ta phân tích hoạt động của con lắc lò xo trên nền Matlab. Với một con lắc nằm ngang, xét theo phương chuyển động, ta có phương trình động lực học đơn giản:

\(-kx=ma.\tag{1}\)

Tham số ban đầu của hệ do người dùng cài đặt:

1
2
3
4
5
6
7
%% INPUT DATA
l0 = 0.1;
A = 0.02;
m = 0.1;
k = 1;
dt = 0.01;
tmax = 100;

Trong vòng lặp tính toán sự tiến hoá của trạng thái theo thời gian, gia tốc \(a\) được gán theo quy luật của phương trình (1):

1
2
3
4
5
6
t = t+dt;
t_array = [t_array t];

a = -k*x/m;
v = v+a*dt;
x = x+v*dt+0.5*a*dt*dt;

Ngoài ra code chương trình cung cấp dưới đây còn giúp vẽ các đồ thị phụ thuộc vào thời gian của li độ, vận tốc, gia tốc. Đặc biệt quỹ đạo pha của chuyển động cũng được diễn tả qua hàm vẽ:

1
graph_x_v = plot(x_array,v_array,'k','linewidth',1);

Video minh hoạ

Code chương trình

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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
 

function conlacloxo
clc
close all
clear all

%% INPUT DATA
l0 = 0.1;
A = 0.02;
m = 0.1;
k = 1;
dt = 0.01;
tmax = 100;

%% Xu ly
x = A;
v = 0;
t = 0;
y = 0;
a = -k*x/m;

t_array = t;
x_array = x;
v_array = v;
a_array = a;

%% FIGURE
fig1 = figure('name','Con lac lo xo','color','black','numbertitle','off');
hold on
set(gcf,'Units','normalized');
set(gcf,'Position',[0 0.1 0.5 0.32]);
fig_loxo = veloxo(-l0,x,(l0+A)/20,10,[0.5 0.5 0.5],'',false);
fig_quanang = plot([x y],'ro','MarkerSize',30,'markerfacecolor','r');

l_vector_danhoi = l0/3;
ve_vec_danhoi = drawArrow('',[x 0],[x-50*k*x*l_vector_danhoi 0],'m',false);

l_vector_v = l0/3;
ve_vec_v = drawArrow('',[x 0],[x-20*l_vector_v*v 0],'y',false);

ht = title(sprintf('t = %0.2f s',t));
set(gca,'color','k','xcolor','w','ycolor','w')

axis equal
axis([-l0 3*A -0.15*l0 0.15*l0]);

fig_x = figure('name','Toa do','color','white','numbertitle','off');
set(gcf,'Units','normalized');
set(gcf,'Position',[0.5 0.1 0.3 0.38]);
graph_x = plot(t_array,x_array,'k','linewidth',2);
xlabel('t [s]');
ylabel('x [m]');

fig_v = figure('name','Van toc','color','white','numbertitle','off');
set(gcf,'Units','normalized');
set(gcf,'Position',[0.5 0.52 0.3 0.38]);
graph_v = plot(t_array,v_array,'y','linewidth',2);
xlabel('t [s]');
ylabel('v [m/s]');

fig_a = figure('name','Gia toc','color','white','numbertitle','off');
set(gcf,'Units','normalized');
set(gcf,'Position',[0.2 0.52 0.3 0.38]);
graph_a = plot(t_array,a_array,'m','linewidth',2);
xlabel('t [s]');
ylabel('a [m/s^2]');

fig_x_v = figure('name','Gia toc','color','white','numbertitle','off');
set(gcf,'Units','normalized');
set(gcf,'Position',[0 0.52 0.2 0.38]);
graph_x_v = plot(x_array,v_array,'k','linewidth',1);
xlabel('x [m]');
ylabel('v [m/s]');

%% CALCULATION

while t<tmax
t = t+dt;
t_array = [t_array t];

a = -k*x/m;
v = v+a*dt;
x = x+v*dt+0.5*a*dt*dt;

x_array = [x_array x];
v_array = [v_array v];
a_array = [a_array a];

figure(fig1);
veloxo(-l0,x,(l0+A)/20,10,[0.5 0.5 0.5],fig_loxo,true);
set(fig_quanang,'xdata',x,'ydata',y);
set(ht,'string',sprintf('t = %0.2f s',t));
drawArrow(ve_vec_danhoi,[x 0],[x-50*k*x*l_vector_danhoi 0],'m',true);
drawArrow(ve_vec_v,[x 0],[x+20*l_vector_v*v 0],'y',true);

figure(fig_x);
set(graph_x,'xdata',t_array,'ydata',x_array);

figure(fig_v);
set(graph_v,'xdata',t_array,'ydata',v_array);

figure(fig_a);
set(graph_a,'xdata',t_array,'ydata',a_array);

figure(fig_x_v);
set(graph_x_v,'xdata',x_array,'ydata',v_array);

pause(0.002);
end

end

%% Ham ve lo xo
function loxo = veloxo(x1,x2,A,q,color,fun,velai)
dis = abs(x2-x1);
lambda = dis/q;

N = (q+1)*2;
x_lx = linspace(x1-lambda/4,x2,N);
x_lx(1) = x_lx(1)+lambda/4;
x_lx(N) = x_lx(N)-lambda/4;
y_lx = zeros(1,N);

for i = 2:N-1
if mod(i,2)==0
y_lx(i) = A;
else
y_lx(i) = -A;
end
end

if velai == false
loxo = plot(x_lx,y_lx,'y','linewidth',2);
loxo.Color = color;
else
set(fun,'xdata',x_lx,'ydata',y_lx,'color',color);
end
end

%% Ham ve mui ten
function hArrow = drawArrow(fun,p0,p1,color,draw)
% drawArrow(p0,p1)
% draw: boolean
if nargin == 2
color = 'k';
end

% Parameters:
W1 = 0.08; % half width of the arrow head, normalized by length of arrow
W2 = 0.014; % half width of the arrow shaft
L1 = 0.18; % Length of the arrow head, normalized by length of arrow
L2 = 0.13; % Length of the arrow inset

% Unpack the tail and tip of the arrow
x0 = p0(1);
y0 = p0(2);
x1 = p1(1);
y1 = p1(2);

% Start by drawing an arrow from 0 to 1 on the x-axis
P = [...
0, (1-L2), (1-L1), 1, (1-L1), (1-L2), 0;
W2, W2, W1, 0, -W1, -W2, -W2
];

% Scale,rotate, shift and plot:
dx = x1-x0;
dy = y1-y0;
Length = sqrt(dx*dx + dy*dy);
Angle = atan2(-dy,dx);
P = Length*P; %Scale
P = [cos(Angle), sin(Angle); -sin(Angle), cos(Angle)]*P; %Rotate
P = p0(:)*ones(1,7) + P; %Shift

% Plot!
if draw == 0
hArrow = patch(P(1,:), P(2,:),color);
hArrow.EdgeColor = color;
else
set(fun,'xdata',P(1,:),'ydata',P(2,:),'EdgeColor',color);
end

end