使用MATLAB计算动力总成悬置系统固有频率和主振型

在进行汽车相关设计计算时，经常需要计算动力总成悬置系统的固有频率以及主振型。本文试图编写MATLAB脚本来完成动力总成悬置系统固有频率和主振型的计算任务。

一、实现步骤

1.1 矩阵计算

1. 计算质量矩阵M（1）
2. 计算各悬置刚度矩阵k（2）
3. 计算角度转置矩阵T（4）
4. 计算位置转移矩阵B（3）
5. 累加计算刚度矩阵K（5）

M=\begin{bmatrix}
m&0&0&0&0&0\\
0&m&0&0&0&0\\
0&0&m&0&0&0\\
0&0&0&I_{xx}&-I_{xy}&-I_{xz}\\
0&0&0&-I_{xy}&I_{yy}&-I_{yz}\\
0&0&0&-I_{xz}&-I_{yz}&I_{zz}\\
\end{bmatrix}\tag{1}

k_i=\begin{bmatrix}
k_{ui}&0&0\\
0&k_{vi}&0\\
0&0&k_{wi}\\
\end{bmatrix}\tag{2}

B_i=\begin{bmatrix}
1&0&0&0&Z_i&-Y_i\\
0&1&0&-Z_i&0&X_i\\
0&0&1&Y_i&-X_i&0\\
\end{bmatrix}\tag{3}

T_i=\begin{bmatrix}
cos\alpha_{ui}&cos\beta_{ui}&cos\gamma_{ui}\\
cos\alpha_{vi}&cos\beta_{vi}&cos\gamma_{vi}\\
cos\alpha_{wi}&cos\beta_{wi}&cos\gamma_{wi}\\
\end{bmatrix}\tag{4}

K=\sum_{i=1}^nB_i^\mathrm{T}T_i^\mathrm{T}K_iB_iT_i\tag{5}

（3）根据实际情况，如果不是相对于中心位置的可能是以下表达式（参考文献[3]中内容，程序中未给出）

B_i=\begin{bmatrix}
1&0&0&0&Z_i-Z_0&Y_0-Y_i\\
0&1&0&Z_0-Z_i&0&Y_i-Y_0\\
0&0&1&Y_i-Y_0&X_0-X_i&0\\
\end{bmatrix}\tag{3-2}

（4）还有一种表达形式，只有两个角度（参考文献[1]中内容，程序中未给出）

T_i=\begin{bmatrix}
cos\theta_{1i}&0&sin\theta_{1i}\\
sin\theta_{1i}\theta_{2i}&cos\theta_{2i}&-cos\theta_{1i}cos\theta_{2i}\\
-sin\theta_{1i}&sin\theta_{2i}&cos\theta_{1i}cos\theta_{2i}\\
\end{bmatrix}\tag{4-2}

1.2 频率计算和主振型

1. 由M和K求得矩阵，取特征值，利用ω=2πf求得固有频率f
2. [v,d]=eig来求解主振型

二、代码与说明

clc;
clear;
%创建参数变量：
% 悬置数目
num=4;
% 质量(kg)
m=277.33;
%创建参数数组
% 惯性参数(xx,yy,zz,xy,yz,zx)
aj=[16.87,7.92,16.55,-2.47,2.95,-0.95];
% 悬挂位置(x1,y1,z1;x2,y2,z2;x3,y3,z3...)
al=[-9.3,209.5,-136.5;20.5,-436.3,-156.5;-359.2,12.5,164.5;584.9,37.1,239.5]'*10^-3;
% 初始动刚度值
ak=[120,35,100;155,90,150;100,70,110;80,42,115]'*10^3;
% 初始角度值(如果多于三个就再加入，后续程序里会乘以cos)
aa{1}=[0,90,90;90,0,90;90,90,0];
aa{2}=[0,90,90;90,0,90;90,90,0];
aa{3}=[0,90,90;90,0,90;90,90,0];
aa{4}=[0,90,90;90,0,90;90,90,0];
%计算质量矩阵
M=[diag([m,m,m]),zeros(3);zeros(3),[aj(1),-aj(4),-aj(6);-aj(4),aj(2),-aj(5);-aj(6),-aj(5),aj(3)]];
%计算刚度矩阵
% 计算需要的各个矩阵
for i=1:num
	%计算ki
	K{i}=diag([ak(1,i),ak(2,i),ak(3,i)]);
	%计算Bi
    B{i}=[eye(3),[0,al(3,i),-al(2,i);-al(3,i),0,al(1,i);al(2,i),-al(1,i),0]];
    %计算Ti
    T{i}=cosd(aa{i});
end
% 计算刚度矩阵K
Ks=0;
for i=1:num
	Ks=Ks+B{i}'*T{i}'*K{i}*T{i}*B{i};
end

%计算频率和主阵型
A=(inv(M))*Ks;
[V,E1]=eig(A);
f=diag(sqrt(E1)/(2*pi));
[f,Sort]=sort(f);
f
V1=V(:,Sort);
v1

三、正确性验证

使用数据为：

表1 动力总成的质量和惯性参数(质量单位为kg，惯性参数单位为kg·m²)

以上数据实验时可采用三线摆测量。

表2 悬置的初始位置(单位：mm)

表3 初始刚度值(单位：N/mm)

在MATLAB 2020a中运行以上代码，输出以下结果：

f =

   5.133572625279529
   6.154126964783480
   6.715844929048983
   8.495706155632371
   9.802800097989474
  15.708570201117141


V1 =

  -0.002575595333373   0.039466990231208   0.604485422397946  -0.042057586371625   0.092059939681550  -0.001923323097360
  -0.984141131745956   0.022430044492999  -0.005879380718305   0.005788641434482   0.002382522157190  -0.000381130252290
  -0.031125123189836  -0.314088774322631   0.080306281140817   0.217013382019956   0.029371516378800   0.011868177386655
  -0.172570999898838  -0.903582353870596   0.091580682327157  -0.951268659010678  -0.213972175366066   0.203700986952554
  -0.024014007619249  -0.245061865768309   0.069605796709097  -0.017779503634877  -0.102534614271487  -0.955297284911824
   0.011550032591457  -0.150900971999743  -0.784143940797264  -0.214183731216263   0.966622556205327  -0.213935287617928

利用某软件对以上数据进行计算，进行对比

可以发现计算结果近似。

参考文献：

[1]曾令贤.用MATLAB计算发动机悬置系统的固有频率和主振型[J].汽车科技,2005(04):27-29.
[2]王星. 电动汽车动力总成悬置系统优化设计[D].合肥工业大学,2015.
[3]薛华,刘志强,刘岩,苏迎.基于Matlab的动力总成悬置系统解耦优化[J].噪声与振动控制,2015,35(02):65-68.

本人为初学者，不能保证计算的准确性，敬请谅解！