四自由度时变系统数值方法求解讨论

咕噜噜

一个四自由度的振动方程：
[M]{x''}+[C(t)]{x'}+[K(t)]{x}={F}
其中刚度和阻尼是时变的，且为离散值，但是不知道刚度和阻尼是否可以写成时间t的函数，下面程序为luoluo自编程序，km,cm为不同时刻的刚度和阻尼值
function f=myMulti(t,x)
%初始数据
m1=0.96;I1=4.3659e-4;m2=2.88;I2=8.3602e-3;  %质量和转动惯量
z1=19;z2=48; %主动轮和从动轮的齿数
k1y=6.56e7;k2y=6.56e7;c1y=1.8e5;c2y=1.8e5;   %支承刚度和阻尼
%km=1.3e6;cm=3.62;  %啮合刚度和阻尼
n=300;
[km,cm]=meshkc(n); % 一个啮合周期内的km和cm值
f1=30;   %输入轴的转动频率
T1=1/f1; %输入轴旋转一周的时间
Tz=T1/z1; %一个啮合周期
nsum = n*z1; %总的采样点数
r1=0.02834;r2=0.07160;  %大、小齿轮基圆半径
L=0.016;    %齿宽
Cr=1.6456;  %重合度
M1=11.9;M2=48.8;    %输入输出扭矩

%将原微分方程组化为一阶方程组，并写成myMulti.m
f=[x(2);(-c1y*x(2)-k1y*x(1)-km*(x(1)+r1*x(5)-x(3)+r2*x(7))-cm*(x(2)+r1*x(6)-x(4)+r2*x(8)))/m1;...
    x(4);(-c2y*x(4)-k2y*x(3)+km*(x(1)+r1*x(5)-x(3)+r2*x(7))+cm*(x(2)+r1*x(6)-x(4)+r2*x(8)))/m2;...
    x(6);(-r1*(km*(x(1)+r1*x(5)-x(3)+r2*x(7))+cm*(x(2)+r1*x(6)-x(4)+r2*x(8)))-M1)/I1;...
    x(8);(r2*(km*(x(1)+r1*x(5)-x(3)+r2*x(7))+cm*(x(2)+r1*x(6)-x(4)+r2*x(8)))-M2)/I2];

求啮合刚度和啮合阻尼的函数：
function[km,cm] = meshkc(n)
%meshkc 计算齿轮副的啮合刚度及阻尼
%
%输入
%n 在一个啮合周期内的计算点数（缺省时计算300点）
%
%输出
%km 随时间变化的啮合刚度值（n*2）
%cm 随时间变化的啮合阻尼值(n*2)

%初始参数
z=[19,48];xi=[2.26,0.0];sigma=0.0;alpha=22.6*pi/180.0;
alpha0=20.0*pi/180.0;rb1=0.02834;
omega=[2*pi*30,74.6128];epsilon=1.6456;
pn=419.9012;
Tz=1.7544e-3;E=2.068e11;mu=0.3;

if nargin<1, n=300;
end
for i=1:2
if (z(i)>135),z (i)=135;
end
if (z(i)<25)
A(i)=11.05+0.166*(abs(z(i)-25))^0.2;
B(i)=0.32+0.05*(abs(z(i)-25))^0.35;
else
A(i)=11.05-0.166*(abs(z(i)-25))^0.2;
B(i)=0.32-0.05*(abs(z(i)-25))^0.35;
end
C(i)=(abs(z(i)-17))^0.45/38.4+1.22;
end
t01=((omega(1)+omega(2))*tan(alpha)-sqrt((z(2)+2+2*xi(2)-2*sigma)^2*...
omega(2)^2/(z(1)^2*(cos(alpha0))^2)-omega(1)^2))/omega(1)/omega(2);

nn=fix(n/epsilon);
t=0.0;
dt=Tz/(n-1);
for m=1:n
   km(m,1)=t;
   cm(m,1)=t;

   for i=1:2
    %轮齿的弯曲刚度(kw(i)(i=1,2));
    lamda(i)=0.5*(z(i)+2+2*xi(i))-0.5*sqrt(1+omega(i)^2*(t+t01)^2)*...
             z(i)*cos(alpha0);
    kw(i)=E*(1+xi(i))^B(i)*(1+lamda(i))^C(i)/A(i);

    %轮齿的接触刚度(kc(i)(i=1,2));
    rho(1)=omega(1)*(t+t01)*rb1;
    rho(2)=((omega(1)+omega(2))/omega(2)*tan(alpha)-omega(1)*(t+t01))*rb1;
    if(i==1),j=2;else,j=1;
    end
    kc(i)=pi*E/((1-mu^2)*(log(pi*E*rho(i)*(rho(1)+rho(2)))-...
          log(2*(1-mu^2)*rho(j)*pn)+0.814));

    %单齿刚度(k(i)(i=1,2));
    k(i)=kw(i)*kc(i)/(kw(i)+kc(i));
  end
  if(m<=nn) %单齿啮合区内啮合刚度
    km(m,2)=k(1)*k(2)/(k(1)+k(2));
  else    %双齿啮合区内啮合刚度
    km(m,2)=k(1)*k(2)/(k(1)+k(2))+km(m-nn,2);
  end
  t=t+dt;
  gxi=0.07;
  r1=55.25/2000;r2=182.75/2000;J1=3.698332e-4;J2=15.257864e-3;
  cm(m,2)=2*gxi*sqrt(km(m,2)*r1^2*r2^2*J1*J2/(r1^2*J1*r2^2*J2));
end
此为典型的变系数线性系统的求解，要求求出其数值解，目前提议可以用ode45来求解，但是由于阻尼是离散值，所以出现了一定的问题，也有其他方法但是不是很可行。大家可各述高见！！

咕噜噜

不知道luoluo最近上来没有，你觉得把阻尼和刚度循环如何

octopussheng

如果能够得到阻尼和时间的关系，如写成一一对应的向量形式，则建议使用下面的方法处理。

下面是一个求解含一个时变项的常微分方程求解方法：
function  xdot = fun1 ( t , x , beta , omega , A , w0 , theta )
% FUN1.M： 时变微分方程例子
% 时变项是 A * sin ( w0 * t – theta )
xdot ( 2 ) =  beta * x (2 ) - omega ^2 * x (1) + A * sin ( w0 * t – theta ) ;  
xdot ( 1 ) = x ( 2 ) ;
xdot = xdot ( : );
% 使得xdot是一列

下面是求解的函数：
beta = .1 ;
omega = 2 ;
A = .1 ;
w0 = 1.3 ;
theta = pi /4 ;
X0 = [ 0 1]’;
t0 = 0 ;
tf = 20 ;
options = [];
[ t , y] = ode23 ( @fun1, [t0 , tf ], X0, options, beta, omega, A, w0, theta );
plot ( t , y );

可以参照一下试试！
我用这个方法解决了一个含有随机数据的微分方程的求解问题！

hohoo

咕主任 这个程序中的阻尼和刚度计算公式是怎么得到的？

无水1324

楼主辛苦了。
还有这个帖子主要是集中讨论多自由度系统，特别是时变参数的多自由度系统的数值解法。至于其中参数的来由、出处，只好由楼主说明。

  希望大家在多自由度的系统求解上发表你的看法！

咕噜噜

这个是实验所得到的，不用考虑，只把他当做离散值吧

咕噜噜

时变刚度和阻尼系统的参数求解用ode45或者ode23都不是问题，很简单，这里主要是考虑了楼主的系统中阻尼和刚度是离散值

octopussheng

对啊，就是把阻尼刚度看成随时间的离散值，如果不能得到一个随时间的离散序列，那也没办法求解了，呵呵

咕噜噜

主要是看这个随时间的离散序列是如何的，luoluo这几天没上来，也不知道她的具体是如何的阻尼和刚度离散值，但是我想应该是可以形成一个时间的离散序列的吧

octopussheng

肯定是这样的吧，不然他这个根本没法求解了，呵呵，方程参数都没法写进去

hnczf739892267

学习了，学习了