离合器预减振超载造成变速箱怠速异响matlab与python仿真分析

1. 概述

变速箱空挡噪声产生于车辆静止，发动机怠速运转，变速箱位于空挡，离合器在接合状态下，在离合器壳和变速箱壳体处可以听见“嘎啦、嘎啦”的响声。虽然变速箱空挡噪声对变速箱寿命影响不大，但是对用户的体验影响很大，如何避免变速箱出现空挡噪声是汽车NHV工程师面临的重要挑战之一。

变速箱产生空挡噪声的主要原因是发动机的角速度波动，产生角速度波动的主要原因：一方面是由于气缸内的燃气压缩产生不规则的脉冲扭矩；另一方面是由于活塞往复运动产生周期性的惯性力矩。周期性变化的发动机角加速度经离合器减振器传给变速箱输入轴，在花键齿间和齿轮齿间产生冲击，该冲击又经变速箱轴承传给变速箱壳体，变速箱壳体振动，发出噪声。

为减小变速箱空挡噪声，离合器减振器通常带有预减振级，刚度较小，以减小发动机传给变速箱的振动，但是如果预减振级的转角设计的过小，则其最大传递扭矩过小，在发动机怠速情况下，经发动机传给变速箱的扭矩可能超过其最大传递扭矩，从而离合器减振器的第二级弹簧参与工作，造成振动增大，发出噪声。因而离合器减振器匹配不合适，造成变速箱输入轴加速度过大，也是变速箱空挡噪声产生的重要原因。

为研究离合器减振器对变速箱空挡噪声的影响，本文建立了发动机、离合器、变速箱的简化动力学模型，对系统的振动情况进行了仿真计算，并对一变速箱空挡噪声案例进行对比分析，验证了仿真计算的可行性，为变速箱离合器匹配设计提供依据。

1. 变速箱空挡噪声的理论分析

齿轮噪声产生的临界条件是：相啮合的齿轮或花键的惯性扭矩大于其上的阻尼力矩。

当时，无负荷的啮合处便发出噪声。

这里，T为从动件上阻尼力矩（Nm）,J为从动件的转动惯量（包括附属于从动件的所有对转动惯量有影响的部分）（kgm2）, θ 为从动件的角加速度（rad/sec2）。

图1为变速箱空挡噪声产生的原理图，发动机经离合器将扭矩及运动传递给变速箱输入轴，然后又经过主动齿轮传给从动齿轮。当主动齿轮加速时，会推动从动齿轮加速，主从动齿轮总是在工作齿面（驱动侧）啮合，此时不会出现齿轮敲击声；当主动齿轮减速时，如果从动轮的阻尼力矩T小于从动齿轮的惯性力矩的绝对值，则从动齿轮不能及时减速，主从动齿轮会从工作齿面脱开，并在非工作齿面（背隙侧）接触，发生敲击声。

理想情况下，飞轮端转速波动经离合器减振器后传到变速箱输入轴端，转速波动会减小，变速箱主从动齿轮的转速波动较小，从动齿轮的角加速度θ较小，从动齿轮的惯性力矩的绝对值较小，不会发生空挡噪声。但是，如果离合器减振器设计不当，则可能使进入变速箱输入轴的转速波动加大，使从动齿轮的惯性力矩的绝对值大于其阻尼力矩，造成主从动齿轮从工作齿面脱开，在非工作齿面接触，产生敲击声。

另外，变速箱的油温对变速箱的齿轮噪声也有影响，当变速箱油温上升，变速箱润滑油的粘度降低，齿轮接触表面的阻尼减小，造成阻尼力矩减小，齿轮噪声加大，这是变速箱有时候在冷箱不出现空挡噪声，热箱后出现空挡噪声的原因。

图1 变速箱空挡噪声产生原理

1. 仿真计算分析

离合器减振器的性能对变速箱空挡噪声有重要的影响，为了分析其影响，可设计发动机、离合器和变速箱的动力学模型。

考虑到发动机怠速、变速箱位于空挡时，发动机转动惯量比离合器、变速箱大的多，假设发动机运动不会受到离合器及变速箱的影响，可以将发动机的飞轮转速作为动力学模型的输入。另外，为了简化模型，我们不需要求哪一对齿轮出现敲击，将所有与输入轴连接的件等效地换算到一轴上，即求各零件（离合器从动片、各拟合齿轮、中间轴、轴承等）到一轴上的当量转动惯量。这样，只需要考虑离合器减振器的刚度和阻尼，即建立单自由度振动模型。离合器减震器是非线性的，有两级刚度，假设在怠速情况下，离合器减振器要么在第一级刚度下工作，要么在第二级刚度下工作，这样，可进一步将动力学模型简化为线性单自由度模型，其本质为离合器隔振模型，如图2所示。u为发动机飞轮带离合器压盘的转速，x为变速箱输入轴的转速，M为变速箱输入轴及其相连零件的当量转动惯量，c和k为离合器扭转减振器的阻尼及刚度。

图2 离合器隔振模型

其微分方程为：

设系统的输出y(t)=x(t),系统的输入为u(t),对微分方程进行拉氏变换可得到系统的传递函数模型TF：

根据传递函数模型，可由飞轮的转速激励得到变速箱输入轴的转速响应。这里，我们不仅希望得到变速箱输入轴的转速响应，还希望得到变速箱输入轴的加速度响应，以估算从动齿轮的惯性力矩。为此，设计一个新的微分方程：

并转换为现在最常用的状态空间模型SS：

dX=AX+Bv

y=CX+Dv

X为状态变量，

即x1为变速箱输入轴的角位移，x2 为变速箱输入轴的转速。

    

   

    

    D=1/m

y为变速箱输入轴的角加速度。

根据线性微分方程的性质，当vt=u(t)时，

系统的状态为：

系统的输出为：

 

用matlab编程进行仿真计算：

1. 用A、B、C、D矩阵及函数ss()建立状态空间模型

sys=ss(A, B, C, D);

1. 用函数lsim()进行仿真计算

[y,t,X]=lsim(sys,u,t,X0);

X0=[0,0]表示系统初始状态为0

u为试验测得的飞轮转速波动数据（测量值减均值），t为对应的时间

X为系统状态变量，X（：，1）为输入轴转速，X（：，2）为输入轴角加速度

y为系统输出，y为输入轴角加速度导数

1. 利用c*X(:,2)+k*X(:,1)计算变速箱输入轴的转速，利用c*y+k*X(:,2)计算变速箱输入轴的角加速度

 

五、附录

1、  单自由度无阻尼振动程序

k1=7.5/(pi/180);

m1=0.22;

c1=0;

t=0:0.001:1;

f0=705/60*3;

%f_dx0=9*(pi/180)*cos(2*pi*f0*t);

%f_x0=9*(pi/180)/(2*pi*f0)*sin(2*pi*f0*t);

%f1=c1*f_dx0+k1*f_x0;

wn=sqrt(k1/m1);

s=2*pi*f0/wn;

B=9*(pi/180)/(2*pi*f0);

x0=0;

dx0=0;

x=x0*cos(wn*t)+dx0/wn*sin(wn*t)+B/(1-s^2)*(sin(2*pi*f0*t)-s*sin(wn*t));

figure

subplot(3,1,1)

%plot(t,f1)

subplot(3,1,2)

plot(t,x*180/pi)

2、  单自由度有阻尼振动程序

k1=7.5/(pi/180);

m1=0.22;

wn=sqrt(k1/m1);

s_=0.5;

wd=sqrt(1-s_^2)*wn;

c1=2*s_*wn*m1;

t=0:0.001:1;

f0=705/60*3;

f_dx0=9*(pi/180)*cos(2*pi*f0*t);%%%%%%%%%%%%%%%%%

f_x0=9*(pi/180)/(2*pi*f0)*sin(2*pi*f0*t);

f1=0*c1*f_dx0+k1*f_x0;%%%%%%没考虑本阻尼造成的地基力的输入!!!!!!!!!!!!!!

 

 

s=2*pi*f0/wn;

B0=9*(pi/180)/(2*pi*f0);

B=B0/sqrt((1-s^2)^2+(2*s_*s)^2);

 

%%%%%%%%%%%%%%%%%%

if s>1

    B=-B;

end 

%%%%%%%%%%%%%%%%%%%%%%

 

phai=atan(2*s_*s/(1-s^2));

 

x0=0.0;

dx0=0.0;

x=exp(-s_*wn*t).*(x0*cos(wd*t)+(dx0+s_*wn*x0)/wd*sin(wd*t))...

+B*exp(-s_*wn*t).*(sin(phai)*cos(wd*t)+(s_*wn*sin(phai)-2*pi*f0*cos(phai))/wd*sin(wd*t))...

+B*sin(2*pi*f0*t-phai);

 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%https://wenku.baidu.com/view/009915e19b89680203d82503.html

%%%第6页公式10如果能改下就更好了。从9推导10时，默认了一个条件就是X与X0同号，

%%%即为正值。如果阻尼omiga>1(即激励频率大于无阻尼固有频率时),应为负值。

%%%增加一个条件就可以使本公式既可以满足有阻尼振动，也可以满足无阻尼振动。从数学和编程上来说，更合理。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure

subplot(3,1,1)

plot(t,f1)

subplot(3,1,2)

plot(t,x*180/pi)

 

grid on

3、  FFT程序

%%%%%%傅里叶变换/逆变换/短时傅里叶变换%%%%%%

 

data=load('testdata.mat');%

 

t=data.x(:,1);%

y=data.x(:,2);

 

Fs=1000;

t1=min(t):1/Fs:max(t);

y1=interp1(t,y,t1);

%y1=sin(2*pi*10*t1)

subplot(3,1,1);

plot(t,y,'r*',t1,y1)

xlabel('Time(s)');%在第一个窗口画波形

 

sigLength=length(y1);

%

Nshow=sigLength/10;%显示的点数（频率段，因频谱对称只需一半内）

%

n=0:sigLength-1;

Y=fft(y1-mean(y1),length(y1)); %在sigLength这个有限区间内做快速傅立叶变换！！！减均值

%！！！！！减均值只考虑波动的频谱情况！！！

%mag=abs(Y);

f=Fs*n/sigLength;

 

subplot(3,1,2);plot(f(1:Nshow),abs(Y(1:Nshow)));

xlabel('Frequency(Hz)');

%在第一个窗口画率谱,

grid on;

 

 

4、  SS模型仿真程序

k1=7.55/(pi/180);

m1=0.3511;

M=[m1];

K=[k1];

n=size(M,1);

 

wn=sqrt(k1/m1);

s_=0.25;

wd=sqrt(1-s_^2)*wn;

c1=2*s_*wn*m1;

C=[c1];

G=[C,M;M,zeros(n)];

H=[K,zeros(n);zeros(n),-M];

ssA=-G\H;

ssB=G\[eye(n);zeros(n)];%G dX +H X =E u => dX="-G\H" X +"G\E" u

%ssB=[zeros(n),inv(M);inv(M),-inv(M)^2]*[eye(n);zeros(n)];

ssC=[-M\K,-M\C];

ssD=inv(M);%设输出Y= d^2X=[-M\K,-M\C]*[X,dX]+[M]\U

%[t,f_x0,f_dx0]=wav_trapz;

data=load('testdata.mat');%

t=data.x(:,1);%

y=data.x(:,2);

Fs=1000;

t1=min(t):1/Fs:5;%max(t);%到5s?????

y1=interp1(t,y,t1);

y1=(y1-mean(y1))/60;%计算相对转速并换算到RPS

t=t1;U=y1;

%t=t-min(t);%将仿真开始时间调整为0

%f1=C(1)*f_dx0+K(1)*f_x0;

%U=[f1];

sys=ss(ssA,ssB,ssC,ssD);

X0=[0.0,0];%X0=zeros(2*n,1);%X0=[x1_0,...,dx1_0,...]

[Y,t,X]=lsim(sys,U,t,X0);

X(:,3)=Y;

Cs=[K(1)C(1) 0;0 K(1) C(1)];%Y=C*X

Ys=X*Cs';

 

subplot(3,1,1)

plot(t,U*60)

xlabel('飞轮转速f_dx0(rpm)');

grid on

subplot(3,1,2)

plot(t,Ys(:,1)*60)

xlabel('离合器从动盘转速dx(rpm)');

subplot(3,1,3)

plot(t,Ys(:,2))

xlabel('离合器从动盘加速度a(rps^2)');

 

 

5、  单自由度振动理论PPT

 

https://wenku.baidu.com/view/009915e19b89680203d82503.html

第6页公式10如果能改下就更好了。从9推导10时，默认了一个条件就是X与X0同号，即为正值。如果阻尼omiga>1(即激励频率大于无阻尼固有频率时),应为负值。增加一个条件就可以使本公式既可以满足有阻尼振动，也可以满足无阻尼振动。从数学和编程上来说，更合理。

 

6、改用python进行仿真计算

 

import matplotlib.pyplot as plt
from scipy.signal import lti, lsim
import numpy as np
import pandas as pd
#plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
#plt.rcParams['axes.unicode_minus']=False #用来正常显示负号

k1=450/(np.pi/180);#7.55;450
m1=0.2331;#0.066251;#0.2331;
M=np.array([[m1]]);
K=np.array([[k1]]);
wn=np.sqrt(k1/m1);
s_=0.25;
n=M.shape[0]
wd=np.sqrt(1-s_**2)*wn;
c1=2*s_*wn*m1;
C=np.array([[c1]]);
G=np.r_[np.c_[C,M],np.c_[M,np.zeros((n,n))]] #G=[C,M;M,O]
H=np.r_[np.c_[K,np.zeros((n,n))],np.c_[np.zeros((n,n)),-M]] #H=[K,O;O,-M]
ssA=-np.linalg.inv(G)@H; 
ssB=np.linalg.inv(G)@np.row_stack([np.eye(n),np.zeros((n,n))]); #G dX +H X =E u => dX="-G\H" X +"G\E" u 
#ssB=[zeros(n),inv(M);inv(M),-inv(M)^2]*[eye(n);zeros(n)];
#ssC=[-M\K,-M\C];
ssC=np.column_stack([-np.linalg.inv(M)@K,-np.linalg.inv(M)@C])
ssD=np.linalg.inv(M); #设输出Y= d^2X=[-M\K,-M\C]*[X,dX]+[M]\U
#t,f_x0,f_dx0]=wav_trapz;
data = pd.read_csv('testdata.txt', header=None, )
#data.head(10)
Fs=data.values[0]
y=data.values[1:]
y1=(y-np.mean(y))/60;#计算相对转速并换算到RPS
t=np.linspace(0,1/Fs[0]*(np.size(y)-1), num=np.size(y))  #Fs[0]才是采样频率的值！！！
U=y1
sys=lti(ssA,ssB,ssC,ssD);
X0=[0.0,0]#X0=zeros(2*n,1);%X0=[x1_0,...,dx1_0,...]
tout,Y,X=lsim(sys,U,t,X0);
#X(:,3)=Y;
X=np.c_[X,Y]
#Cs=[K(1) C(1) 0;0 K(1) C(1)];%Y=C*X
Cs=np.array([[k1,c1,0],[0,k1,c1]])
Ys=X@Cs.T

plt.figure()
plt.subplot(3,1,1)
plt.plot(t,U*60)
plt.xlabel('speed of flywheel/f_dx0(rpm)'); 
#grid on
plt.grid(alpha=0.3)
plt.subplot(3,1,2)
plt.plot(t,Ys[:,0]*60)
plt.xlabel('speed of clutch/dx(rpm)');
plt.subplot(3,1,3)
plt.plot(t,Ys[:,1])
plt.xlabel('acc. of clutch/a(rps^2)');