基于欧拉伯努利梁模型的本征模分析
61
|
Myron
|
|
0
|
数学,物理学,视频总结

1591 字
|
12 分钟 
标题：Modal Analysis (Eigenvalue Analysis/Free Vibration Analysis) of beam: Theory and MATLAB Coding
链接：https://youtu.be/yhJC16mHNBM
原作者：Dr. Subhajit Mondal

以一个Euler-Bernoulli Beam模型为例，解释本征模分析，包括理论和代码。

本征值问题

下面是本征值问题的最初形态与求解方式。

\begin{aligned}
&\vb*A \va*x = \lambda \va*x \\
&\vb*A \va*x – \lambda \va*x = 0 \\
&\qty(\vb*A – \lambda \vb*I) \va*x = 0 \\
&\mqty|\vb*A – \lambda \vb*I| = 0 \\
\end{aligned}

一般先求本征值$\lambda$，再求解本征向量$\va*x$。比如$\vb*A = \mqty(2 & 1 \\ 1 & 2)$，可转化为$\mqty|2 – \lambda & 1 \\ 1 & 2 – \lambda| = 0$，可以解出本征值与本征向量为如下形式：

\begin{aligned}
\lambda_1 = 1 &\hspace{2em} \va*x_1=\mqty(-1 \\ 1) \\
\lambda_2 = 3 &\hspace{2em} \va*x_2=\mqty(1 \\ 1) \\
\end{aligned}

欧拉伯努利杆物理模型

描述这个杆的物理模型的方程如下，其中$\qty[\,]$表示矩阵。

\begin{aligned}
&\qty[M]\qty{\ddot{u}} + \qty[C]\qty{\dot{u}}+ \qty[K] \qty{\dot{u}}=F\qty{t} \\
\end{aligned}

其中$\qty[M]$是质量矩阵，$\qty[C]$是阻尼项，$\qty[K]$是劲度系数，F是外力。如果去掉外力，考虑初始给一点扰动，然后让他自由振动。其方程退化为如下形式

\begin{aligned}
&\qty[M]\qty{\ddot{u}} + \qty[C]\qty{\dot{u}}+ \qty[K] \qty{\dot{u}}=0 \\
\end{aligned}

此时因为含有阻尼项，求解该方程时会得到复数的特征值和特征向量。$\lambda=\lambda_r + \lambda_i i$，为了得到实数特征值和特征向量，也可以把阻尼项去掉，方程变为如下形式：

\begin{equation}\label{eq:1224_1}
\qty[M]\qty{\ddot{u}} + \qty[K] \qty{\dot{u}}= 0
\end{equation}

此时可以得到方程的解为下面形式：

\begin{equation}\label{eq:1224_2}
u = \qty{\phi} \sin{\omega t}
\end{equation}

其中$\phi$是一个向量，表示特征向量对应的模。将解\eqref{eq:1224_2}代入到原方程\eqref{eq:1224_1}可得：

\begin{equation}
\begin{aligned}
-\qty[M]\qty{\phi}\omega^2\sin{\omega t} + \qty[K]\qty{\phi}\sin{\omega t} &= 0 \\
\qty(-\qty[M]\omega^2 + \qty[K])\qty{\phi}\sin{\omega t} &= 0 \\
\qty(-\qty[M]\omega^2 + \qty[K])\qty{\phi} &= 0 \\
\text{det}\qty(-\qty[M]\omega^2 + \qty[K]) = \mqty|\qty[K] – \qty[M]\omega^2| &= 0 \\
\end{aligned}
\end{equation}

其中$\mqty|\qty[K] – \qty[M]\omega^2| = 0$项类似本征值问题的$\mqty|\vb*A – \lambda \vb*I| = 0$。

mass and stiffness 矩阵

Euler-Bernoulli Beam Finite Element – Deriving the Mass and Stiffness Matrices有矩阵完整推导过程，其中M矩阵的表达式如下所示：

\begin{equation}\label{eq:mass_matrix}
\begin{aligned}
M = \frac{\rho A l}{420}\mqty[
156 & 22l & 54 & -13l \\
22l & 4l^2 & 13l & -3l^2 \\
54 & 13l & 156 & -22l \\
-13l & -3l^2 & -22l & 4l^2 \\
]
\end{aligned}
\end{equation}

K矩阵的表达式如下：

\begin{equation}\label{eq:stiffness_matrix}
\begin{aligned}
K = \mqty[
\frac{12EI}{L^3} & \frac{6EI}{L^2} & -\frac{12EI}{L^3} & \frac{6EI}{L^2} \\
\frac{6EI}{L^2} & \frac{4EI}{L} & -\frac{6EI}{L^2} & \frac{2EI}{L} \\
-\frac{12EI}{L^3} & -\frac{6EI}{L^2} & \frac{12EI}{L^3} & -\frac{6EI}{L^2} \\
\frac{6EI}{L^2} & \frac{2EI}{L} & -\frac{6EI}{L^2} & \frac{4EI}{L} \\
]
\end{aligned}
\end{equation}

杆上每一个有限元节点都有两个自由度，一个是角度$\theta$，一个是距离$y$。根据有限元划分了多少个节点，整个系统就有多少个自由度。最终矩阵的行列数$n$等于最终的自由度数量。

有限元数值方法

我们划分成5个区间，总共6个节点，共计12个自由度，最终合成的矩阵如下：

\begin{equation}\label{eq:stiffness_matrix_big}\small
K_g = \mqty[
\frac{12EI}{L^3} & \frac{6EI}{L^2} & -\frac{12EI}{L^3} & \frac{6EI}{L^2} &
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\frac{6EI}{L^2} & \frac{4EI}{L} & -\frac{6EI}{L^2} & \frac{2EI}{L} &
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
-\frac{12EI}{L^3} & -\frac{6EI}{L^2} &
\frac{24EI}{L^3}&
0 &
-\frac{12EI}{L^3} & \frac{6EI}{L^2} &
0 & 0 & 0 & 0 & 0 & 0 \\
\frac{6EI}{L^2} & \frac{2EI}{L} &
0 &
\frac{8EI}{L} &
-\frac{6EI}{L^2} & \frac{2EI}{L} &
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 &
-\frac{12EI}{L^3} & -\frac{6EI}{L^2} &
\frac{24EI}{L^3}&
0 &
-\frac{12EI}{L^3} & \frac{6EI}{L^2} &
0 & 0 & 0 & 0 \\
0 & 0 &
\frac{6EI}{L^2} & \frac{2EI}{L} &
0 &
\frac{8EI}{L} &
-\frac{6EI}{L^2} & \frac{2EI}{L} &
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 &
-\frac{12EI}{L^3} & -\frac{6EI}{L^2} &
\frac{24EI}{L^3} &
0 &
-\frac{12EI}{L^3} & \frac{6EI}{L^2} &
0 & 0 \\
0 & 0 & 0 & 0 &
\frac{6EI}{L^2} & \frac{2EI}{L} &
0 &
\frac{8EI}{L} &
-\frac{6EI}{L^2} & \frac{2EI}{L} &
0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 &
-\frac{12EI}{L^3} & -\frac{6EI}{L^2} &
\frac{24EI}{L^3} &
0 &
-\frac{12EI}{L^3} & \frac{6EI}{L^2} \\
0 & 0 & 0 & 0 & 0 & 0 &
\frac{6EI}{L^2} & \frac{2EI}{L} &
0 &
\frac{8EI}{L} &
-\frac{6EI}{L^2} & \frac{2EI}{L} \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &
-\frac{12EI}{L^3} & -\frac{6EI}{L^2} & \frac{12EI}{L^3} & -\frac{6EI}{L^2} \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &
\frac{6EI}{L^2} & \frac{2EI}{L} & -\frac{6EI}{L^2} & \frac{4EI}{L} \\
]_\text{12×12}
\end{equation}

可以看到矩阵元素都在对角线附近，称之为Banded matrices。同时矩阵对称，对角线上的元素都为正，可以得到矩阵应该为正定或者半正定矩阵。此时如果不给边界条件的话，其本征值与本征向量应该都为零。也就是杆的运动应该类似刚体平移或者旋转，而没有形变引起的振动。

左侧固定边界，$y=0\, \theta=0$。可以将$K_{G11}$与$K_{G22}$的数值设置得特别大，来模拟其为固定边界。

同理可以得到质量的12×12的矩阵。

matlab代码

下面是构造矩阵求解本征值的部分。

nelm=5;
ndof= 2*nelm+2;
M(ndof,ndof)=0.0;
K(ndof,ndof)=0.0;

roh=7800;
b=39E-3;
d=5.933E-3;
A= b*d;
E=2E11;
I=6.7772E-10;
l=1.0/nelm;
Ke(1,1)= (E*I/l^3)*12;
Ke(1,2)= (E*I/l^3)*6*l;
Ke(1,3)= (E*I/l^3)*(-12);
Ke(1,4)= (E*I/l^3)*6*l; 
Ke(2,2)= (E*I/l^3)*4*l^2;
Ke(2,3)= (E*I/l^3)*(-6*l);
Ke(2,4)= (E*I/l^3)*(2*l^2); 
Ke(3,3)= (E*I/l^3)*12;
Ke(3,4)= (E*I/l^3)*(-6*l); 
Ke(4,4)= (E*I/l^3)*(4*l^2);
Me= (roh*A*l/420)*[156 22*l 54 -13*l;
                   22*l 4*l^2 13*l -3*l^2;
                   54 13*l 156 -22*l;
                   -13*l -3*l^2 -22*l 4*l^2] ;
for i=1:4
    for j=1:4
        Ke(j,i)= Ke(i,j);
    end 
end
p=0;
q=0;
for ne=1:nelm 
   for ii=1:4
       for jj=1:4    
           M(p+ii,q+jj)=M(p+ii,q+jj)+ Me(ii,jj);
           K(p+ii,q+jj)=K(p+ii,q+jj)+ Ke(ii,jj);
        end
   end
   p=p+2;
   q=q+2;
end  
  
   K(1,1)=1E10;
   K(2,2)=1E10;
[evec,eval]=eig(K,M);

下面是把求解出来的本征值，本征向量读取并且作图的部分：

for i=1:ndof;
    eval(i,i)=sqrt(eval(i,i));
    fhz(i,i)= eval(i,i)/(2*pi);
end

for i=1:nelm
    k1=i*2-1;
   V(i,1)=evec(k1,1);
   V(i,2)=evec(k1,2);
   
end
plot(V(:,1),'r--o', 'linewidth',2, 'markersize',5, 'markeredgecolor','r','markerfacecolor','r')
hold on
plot(V(:,2),'g--o', 'linewidth',2, 'markersize',5, 'markeredgecolor','g','markerfacecolor','g')

legend({'Mode 1', 'Mode 2'},'Location','northeast')

title('Eigenvalue analysis of a Beam', 'FontName','TimesNewRoman','fontweight','bold','fontangle',...
'italic','fontsize',16,'color','b')

xlabel('Node Number', 'FontName','TimesNew弹性柱子的振动Roman','fontweight','bold','fontangle',...
'italic','fontsize',16,'color','k')

ylabel('Magnitude', 'FontName','TimesNewRoman','fontweight','bold','fontangle',...
'italic','fontsize',16,'color','k')
hold off

运行结果

本征向量

频率最低的两个模
将精度调高之后，可显示频率更高的模

本征值

这里显示前三个本征值（频率）在不同精度下的直（从小到达排序）：

本征值精度5精度10精度20精度50精度100
$\lambda_1$30.470330.469330.460830.231231.5947
$\lambda_2$191.0483190.9590190.9516190.9232191.0978
$\lambda_3$536.5935534.8089534.6811534.6649534.7132
文章标题：基于欧拉伯努利梁模型的本征模分析
文章作者：Myron
转载链接：https://phyiscs.com/eigenvalue-analysis-of-a-beam.html
beameigenvalueEuler-Bernoullimatlabmodal
上一篇
下一篇

推荐文章

重整化群简介
Wigner 定理
神经网络在磁约束聚变物理上应用汇总
文献阅读：《EuroPED-NN: Uncertainty aware surrogate model》
《托卡马克台基等离子体密度预测模型》
《基于神经网络加速模型的自洽核心-台基输运模拟》
《基于神经网络模型的KSTAR台基区剥离-气球模不稳定性研究》
时空的经典理论(一) — 时空的几何
BOUT++5.1.1对Python的依赖
Tokamak平衡计算（一）——磁面坐标系