在先导1中,我们对保守系统的双球系统进行了仿真,现在我们先从单个小球摆动入手,考虑加入了阻尼后,如果通过拉个朗日分析来得到拉格朗日函数进而得到小球的运动方程。
首先回顾一下,先导1中的拉个朗日方程使用的是保守场中的形式,也就是说所有主动力都是保守力,即保守力的做工与路径无关,通过保守力的性质我们可以得到拉格朗日量L = T - V;但是现在的情况是,出现了一个与速度相关的阻尼力,这个阻尼力不是保守力,所以不能通过拉个拉格朗日的保守场形式,需要做一下变动。
考虑单个小球单摆下,受到重力和阻尼力,重力任然可以归结为势能处理,阻尼力需要通过一般形式拉格朗日方程中的广义力来计算求得。公式推导如下所示:
这里有一个要注意的地方,就是重力我还是将其归结为势能,归类到保守场的拉格朗日方程中,如果没有阻尼力的话,保守场拉格朗日方程右边应该为0,但是现在还有一个非保守力阻尼力,所以现在拉格朗日方程的右边就是这个阻尼力对应的广义力,求解出来带入就得到系统的运行方程了。matlab代码如下所示:
%基于ode45的单小球摆动带阻尼的运动方程
%此处是基于拉个朗日函数分析法给出的微分方程求解
close all;clearvars;clc;
dbstop if error;
global m l g k;
m = 5;
l = 1;
g = 9.8;
k = 2;
[t , y] = ode45(@vdp1,[0 25],[30*pi/180 0]);
x1 = y(: , 1) .* 180/pi;
h1 = figure;
for i = 1 : 1 : length(t)
set(0,'CurrentFigure',h1);
title(num2str(i));
scatter(0 , 0 , '*');hold on;
x = l*sind(x1(i));
y = -1*l*cosd(x1(i));
scatter(x , y);hold on;
plot([0 x(1)] , [0 y(1)] , 'r');hold on;
xlim([-3 3]);
ylim([-3 3]);
pause(0.05);
hold off;
end
function dydt = vdp1(t,y) %y = [phi1 phi1']输入为弧度
global m l g k;
dydt = zeros(2 , 1);
dydt(1) = y(2);
dydt(2) = (-2*k*l/m)*y(2) - (g/l)*sin(y(1));
end