多自由度微分方程的求解

luoluo

原帖由 咕噜噜 于 2007-6-29 21:04 发表

                               
登录/注册后可看大图

首先是meshkc没有定义
其次你编这个程序最好写为m函数文件的形式
f没有定义

我只是拷了一部分代码，那个meshkc是以一个函数的格式写的，f也定义了

luoluo

原帖由 无水1324 于 2007-6-29 20:31 发表
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;. ...

我不太明白你的意思，该怎么指定呢，不要笑哦，我比较笨。:@(

咕噜噜

你这种写法通常是在m函数文件中使用的，m函数中并不直接定义，调用的时候才定义

luoluo

我的函数是这样的
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];

无水1324

就是你就算得时刻t，取该时刻对应的km和cm

luoluo

求啮合刚度和啮合阻尼的函数：
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

appleseed05

齿轮的东西我不太懂，但是对于二阶微分方程处理成一阶的一般是这样（如下图）

然后就可以RK了

luoluo

原帖由 无水1324 于 2007-6-29 21:26 发表
就是你就算得时刻t，取该时刻对应的km和cm

我明白你的意思了，我再改改看吧，谢谢了！

luoluo

原帖由 appleseed05 于 2007-6-29 21:49 发表
齿轮的东西我不太懂，但是对于二阶微分方程处理成一阶的一般是这样（如下图）

然后就可以RK了

谢谢，只不过我的刚度和阻尼是时变的，如果是常数的话我想我也能解出来吧

VibrationMaster

1.楼主说的是变系数, 并非是非线性问题.
2.不管是线性还是非线性,ode45采用的数值解法总是可以的.
3.楼主所说的离散形式,当然在离散时刻上是已知的,.当采用ode45时一般是等步长的,如果这个步长与阻尼和刚度的采样间隔一致,可以采用插值将阻尼和刚度的采样间隔之间插值,以便与ode45的步长一致. 插值方法很多,找本数值方法看一看.

无水1324

原帖由 VibrationMaster 于 2007-6-30 07:12 发表
1.楼主说的是变系数, 并非是非线性问题.
2.不管是线性还是非线性,ode45采用的数值解法总是可以的.
3.楼主所说的离散形式,当然在离散时刻上是已知的,.当采用ode45时一般是等步长的,如果这个步长与阻尼和刚度的 ...

“ 采用ode45时一般是等步长的”ode45时一般是并不是等步长的。这点要注意

无水1324

但是你提出的方法不失为一个好方法

luoluo

其实km和cm的取值我是等间隔取的，每个间隔之间的dt是知道的，在那个meshkc(n)的函数里面可以体现出来

luoluo

无水1324还在吗？我的qq是32677381，我们用qq聊一下吧，我这个人太笨了，尽管明白你说的是什么意思，可是还是不知道该怎么修改，可以加一下你的qq吗？

无水1324

就是在每个计算时刻选取一个km，cm就是了，要增加一个循环