连续体振动讲解

Vibration of Countinuous System
1. Euler —Bernoulli Beam
Euler —Bernoulli 运动过程中没有考虑剪切效应的影响。

Euler —Bernoull Beam ：变形前垂直于梁中心线的截面在变形后仍保持垂直于梁的中心线。

Timoshenko Beam ：Euler —Bernoull 梁中并没有考虑梁的剪切变形，在实际工程中，会存在梁的剪切变形，变形后截面与中心线存在一个夹角，截面的转角变为
y x
θγ∂=
-∂
中性面定义：+++
变形的几何关系：
假设距离中性层距离为h 的层为b b -。

根据平面假设，单元体d x 变形后层面b b -为b b ''-
其中，d θ为变形的角度
b
b
o o
h
C
C '
d θ
o
o '
o '
b 'b 'M
M
R
h
A '
B '
A ϕB
ϕ
M
C
应变的表达式为
()d -d d R h R h
R R
θθεθ+=
=
弯矩的表达式为
2d d d A
A
A E EI
M h A h E A h A R R
σε===
=⎰⎰⎰
其中，I 为截面的惯性矩。

转角的表达式，A 点的转角为A ϕ，B 点的转角为B ϕ
A y x
ϕ∂=
∂
对于B 点，假定转角对位置坐标线性变化，有
2
2
d d A
B A x x
y y x x x ϕϕϕ∂=+
∂∂∂=+∂∂
因此，弯曲的角度d θ表示为
22d d B A y
x x
θϕϕ∂=-=∂
由于梁弯曲变形为小变形，有如下
d d R x θ=
得到
221y R x
∂=∂
得到弯矩的表达式
2M EI x
=∂
1.1 Newton 2th law Equation.
y
x
取长度为L 的梁中的微元体研究，单元体的长度为d x 。

假定受到与位置坐标x 相关的载荷()p x 的作用，考虑到变截面梁，假定截面面积为()A x 。

梁的密度表示为位置坐标x 的函数()x ρ。

单元体受力情况如图++所示。

其中M 为弯矩(bending)，V 为剪力(shear force)
这里认为弯矩M 与剪力V 在位置坐标系下随x 线性变化，因此单元体右端的弯矩和剪力分别表示为如下形式
d right M
M M x x
∂=+
∂ (1)
d right V V x x
=+
∂ (2)
假设梁在y 方向的位移(挠度)为时间t 与位置坐标x 的函数，写成(,)y x t 的表达形式。

在y 方向合力为零，即满足0y F =∑，有
2
2d ()d ()()d Q y Q Q x p x x x A x x x t ρ∂∂⎛
⎫-++= ⎪∂∂⎝
⎭
(3)
式中，()d p x x 为作用在单元体上的载荷；等号右端表示惯性力项。

公式(3)略去d x ，化简为
22()()()y Q x A x p x t x
ρ∂∂+=∂∂
(4)
考虑到Euler —Bernoulli Beam 假设，忽略剪切力的影响，弯矩平衡。

假设单元体左端固定，有弯矩平衡关系
d d d d ()d 02M Q x M x M Q x x p x x x x ∂∂⎛⎫⎛
⎫+--++= ⎪ ⎪
∂∂⎝⎭⎝⎭
(5)
整理公式(5)，有
2d d d d ()d 02
M Q x
x Q x x p x x x x ∂∂--+=∂∂ (6)
略去2d x 项有
M
Q x
∂=∂ (7)
将式(7)带入式(4)有
2222()()()y M
x A x p x t x
ρ∂∂+=∂∂
(8)
弯矩M 与挠度y 的关系为
22
(,)
()y x t M EI x x ∂=∂
(9)
式中，E 为杨氏模量；()I x 为梁截面对x 轴的惯性矩，对于矩形截面梁，表达式
为3
()12
bh I x =。

将弯矩表达式(9)带入式(8)，有
22
222
2(,)()()()()y y x t x A x EI x p x t x
x ρ⎡⎤
∂∂∂+=⎢⎥∂∂∂⎣⎦
(10)
公式(10)即为梁强迫振动的微分方程形式，若梁受到的外力同时随时间t 与位置坐标x 的影响，方程右侧的力可以写成(,)p x t 的形式。

方程(10)总阶次为6，其中包含对空间坐标y 的4阶导数和对时间t 的2阶导数，求解需要6个条件。

在动力学领域，采用边界条件提供求解条件。

提供四个边界条件和两个初始条件，可以确定方程(10)的解。

对于求解梁的固有频率问题和自由振动问题，方程(10)右端的力()0p x =。

对于方程(10)的求解采用一些列正交基函数序列进行求解，即所谓的Garlerkin 法。

求解微分方程采用变量分离法，基于这种方法，将式（10）的解写成如下形式
(,)()()y x t Y x q t = (11)
其中，()Y x 表示与位置坐标x 有关的函数；()q t 表示与时间相关的函数。

将式(11)带入式(10)并令()0p x =，有
44
d ()
()()()()()()0d Y x x A x Y x q t EI x q t x
ρ+= (12)
式(12)可以写成如下形式
44d ()
()()d ()()()()
Y x q t EI x x q t x A x Y x ρ=- (13)
公式（13）中，方程左侧与时间t 有关，方程右侧与位置坐标x 有关。

若要求方程成立，方程两端等于同一常数，令常数为2ω-。

则方程(13)有
2
442()()
d ()
()d ()()()
q t q t Y x EI x x x A x Y x ωωρ⎧=-⎪⎪
⎨⎪-
=-⎪⎩ (14)
进一步写成如下形式
2424
()()0
()d ()()()()()0()d q t q t a Y x EI x x A x Y x b x
ωωρ+=-=
(15)
观察方程（15-a ），发现方程的解应为三角函数组合形式
12()sin sin q t A t A t ωω=+
(16)
这里进行分析时，假定梁是各项同性材料，密度为常数，梁截面积不发生变化，即(),()x A x A ρρ==，方程（15-b ）可以写成如下形式
444
d ()
()0d Y x Y x x
β-= (17)
其中，
2
β=
方程(17)的解可以表达为
()sx Y x e =
(18)
带入方程（17）得到特征方程为
440s β-=
(19)
求解方程(19)得到4个特征根为
1,23,4,s s i ββ=±=±
(20)
因此，方程(17)的解的表达式为
()1234x x i x i x Y x C e C e C e C e ββββ--=+++
(21)
根据Euler 公式如下关系
cosh sinh cos sin x i x
e x x e
x i x
ββββββ±±=±=± (22)
将欧拉公式（22）带入方程（21）得到梁模型的振型函数为
()cosh sinh cos sin Y x A x B x C x D x ββββ=+++
(23)
式中，,,,A B C D 表示常数，可以通过几何边界条件确定。

因此，梁模型的振动方程的解的形式为
()()12(,)cosh sinh cos sin sin sin y x t A x B x C x D x A t A t ββββωω=++++ 可以通过方程（16）和方程（23）确定梁的固有频率和振型。

1.2 梁模型的几何边界条件 Boundary condition
常见的梁模型的支撑方式有：两端固支、两端自由、两端简支、固支自由、固支简支等，这些支撑方式可以简化为几何边界条件，梁模型常见的几何边界条件有：固定边界条件、自由边界条件、铰支边界条件等。

(1) 固定端：
位移和转角为零0x or x L ==
(,)0
(,)
0y x t y x t x
=∂=∂
(2) 自由端：
弯矩和剪力等于零0x or x L ==
22
2
2(,)()0(,)()0y x t EI x x
y x t EI x x x ∂=∂⎡⎤∂∂=⎢⎥∂∂⎣⎦
(3) 简支端：
位移和弯矩为零0x or x L ==
22
(,)()0
(,)0
y x t EI x x y x t ∂=∂=
1.3 悬臂梁cantilever beam
悬臂梁是指一端固定，一端自由支撑的梁模型，如图所示
cantilever beam
Boundary
图 悬臂梁示意图
悬臂梁的边界条件为：
固定端位移和转角为零，自由端弯矩和剪力为零。

假设惯性矩()I x I =，有边界条件表达式
23
23
(,)
Camp :
(,)0,
0(,)(,)Free:
0,()0
y x t y x t x
y x t y x t EI EI x x x ∂==∂∂∂==∂∂ 将边界条件带入梁模型的振型函数表达式()Y x ，式（23）
x Y
A C A C
==+⇒=-
[]()00
sinh cosh sin cos 00
x x dY A x B x C x D x dx B D B D
ββββββ==⎛⎫
=+-+= ⎪⎝⎭⇒+=⇒=-
[][][]22
2d cosh sinh cos sin 0
d cosh cos sinh sin 0x L
y A L B L C L D L x C L L D L L βββββββββ=⎛⎫=+--= ⎪⎝⎭⇒+++= [][][]33
3d sinh cosh sin cos 0
d cosh sin cosh cos 0
x L
y A L B L C L D L x C L L D L L βββββββββ=⎛⎫=++-= ⎪⎝⎭⇒-+-= (24)
通过(24)最后两个公式，可以得到
cosh cos sinh sin 0cosh sin cosh cos 0L L L L C L L L L D ββββββββ++⎡⎤⎧⎫⎧⎫=⎨⎬⎨⎬⎢⎥--⎣⎦⎩⎭⎩⎭
(25)
方程(25)的有解的条件时系数行列式为零，进一步化简为
cosh cos 1L L ββ= (26)
方程(26)为悬臂梁的频率方程，求解方程(26)可以得到悬臂梁的特征根值 方程(26)的解可以采用作图法 作图法：
05101520
-1
-0.5
0.5
1
牛顿迭代法进行求解。

令()cosh cos 1F L L L βββ=-，针对悬臂梁的阶次为自然数，即i N β=，梁的长度L 为已知。

因此，方程()0i F L β=的根为离散点。

表 悬臂梁的前5阶特征值
1L β 2L β 3L β 4L β 5L β 1.999 4.729
7.855
10.996
14.137
根据2
β=
，得到
(1,2,3)r r r ωβ==⋅⋅⋅ (27)
通过求解得到的特征根，可以得到系数C 与D 的比值r ξ.
r cosh cos sinh sin cosh sin cosh cos r r r r r r r r r
L L L L C D L L L L ββββξββββ-+⎛⎫
==-=- ⎪-+⎝⎭
(28)
因此，悬臂梁的第r 阶振动的固有频率r ω对应的振型函数为
()()sinh sin cosh cos r r r r r r Y x x x x x ββξββ=-+-
(29)
振型函数的绘制
令1L =，则位置坐标的取值范围为[]0,1x ∈，各阶的特征根为
表 悬臂梁的前5阶特征值
1β 2β 3β 4β 5β
1.999 4.729 7.855 10.996 14.137
绘制的各阶振动的振型如下图所示
0.20.4
0.60.81
-3-2
-1
1
2
3
x
Y (x )
各阶振型图
图 悬臂梁的各阶振动示意图
x
Y (x )
1阶振型图
x
Y (x )
2阶振型图
x Y (x
)3阶振型图
x
Y (x
)4阶振型图
x
Y
(
x )
第5阶振型图
附MATLAB 程序：
A.作图法
clear,clc
x=0:0.001:20;
y1=cos(x);
y2=-1./(cosh(x));
hold on
plot(x,y1,'-r')
plot(x,y2,'-b')
B. Newton 法程序
%二分法求解非线性方程
%调用eval 命令
clear,clc
syms f x; %定义变量
f=input('输入函数f(x)=')
A=input('输入求根区间[a,b]=')
e=input('输入精度要求e=')
%求解
while (A(2)-A(1))>e
c=(A(1)+A(2))/2;
x=A(1);
f1=eval(f);
x=c;
f2=eval(f);
if (f1*f2)>0
A(1)=c;
else A(2)=c;
end
end
c=(A(1)+A(2))/2;
fprintf('c=%0.8f\n',c)
C. 振型函数绘图
clear,clc
% 令bi=beta
b1=1.875;
b2=4.694;
b3=7.855;
b4=10.996;
b5=14.137;
b=[b1 b2 b3 b4 b5];
xi1=-(sinh(b1)+sin(b1))/(cosh(b1)+cos(b1)); xi2=-(sinh(b2)+sin(b2))/(cosh(b2)+cos(b2)); xi3=-(sinh(b3)+sin(b3))/(cosh(b3)+cos(b3)); xi4=-(sinh(b4)+sin(b4))/(cosh(b4)+cos(b4)); xi5=-(sinh(b5)+sin(b5))/(cosh(b5)+cos(b5)); xi=[xi1 xi2 xi3 xi4 xi5];
x=0:0.01:1;
n=length(x);
Y=zeros(n,5);
y0=zeros(n,1);
for i=1:5
for j=1:n
Y(j,i)=sinh(b(i)*x(j))-sin(b(i)*x(j))+xi(i)*(cosh(b(i)*x(j))-cos(b(i)*x(j)));
end
end
figure (1)
hold on
plot(x,Y(:,1),'-.g')
plot(x,Y(:,2),'-c')
plot(x,Y(:,3),'-.m')
plot(x,Y(:,4),'-b')
plot(x,Y(:,5),'-r')
plot(x,y0(:,1),':b')
box on
xlabel('x')
ylabel('Y(x)')
title('各阶振型图')
legend('第1阶振型','第2阶振型','第3阶振型','第3阶振型','第3阶振型','平衡位置')
figure (2)
hold on
plot(x,Y(:,1),'-b')
plot(x,y0(:,1),':r')
box on
xlabel('x')
ylabel('Y(x)')
title('1阶振型图')
legend('第1阶振型','平衡位置')
axis([0 1 -3 3])
figure (3)
hold on
plot(x,Y(:,2),'-b')
plot(x,y0(:,1),':r')
box on
xlabel('x')
ylabel('Y(x)')
title('2阶振型图')
legend('第2阶振型','平衡位置')
figure (4)
hold on
plot(x,Y(:,3),'-m')
plot(x,y0(:,1),':b')
xlabel('x')
ylabel('Y(x)')
box on
title('3阶振型图')
legend('第3阶振型','平衡位置') figure (5)
hold on
plot(x,Y(:,4),'-b')
plot(x,y0(:,1),':r')
xlabel('x')
ylabel('Y(x)')
box on
title('4阶振型图')
legend('第4阶振型','平衡位置') figure (6)
hold on
plot(x,Y(:,5),'-r')
plot(x,y0(:,1),':b')
xlabel('x')
ylabel('Y(x)')
box on
title('第5阶振型图')
legend('第5阶振型','平衡位置')。