[讨论]参数激励非线性系统分析

无水1324

原帖由 xiaojin831108 于 2007-3-22 17:41 发表
裂纹或松动的转子系统列出系统的微分方程后能用matlab的ode45解出数值解不？
微分方程是强非线性的，能不能用matlab的.m文件编程求解啊？
还有里面的油膜力计算式表达好复杂啊？有没有简化的或有谁已经计算出 ...

强非线的基本上都可以用ode45求解，油膜力复杂，只有有显式表达式就可以了，这方面我没有做过！看了一些文章，觉得应该做出没有问题！

无水1324

原帖由 21172485 于 2007-3-23 14:28 发表
在陈予恕的《非线性振动》中提及：对于强非线性系统，谐波平衡法是一个有效的方法，在很多情况下可以获得满意的结果
这个系统随着激励w的不同分很多情况，有共振情况和非共振情况或者w=0的自由振动情况

谐波平衡法 是一种比较好的方法，可以求出低次近似解析解，但是在求高次时就会出现“矛盾问题”（详见胡海岩的文集），而且对所求解进行稳定性分析，必须借助于线性化或者其它的方法。

xiaojin831108

%function [dx,dy]=sdzz1(t,x,y)
clear all;close all;
t=0:20;
x=[0,0,0,0,0,0];
y=[0,0,0,0,0,0,0,0];
dx=zeros(6,1)
dy=zeros(8,1)
%参数设置
m1=4
m2=32.1
m3=50
R=25
L=12
c=0.11
miu=0.018
c1=1050
c2=2100
omiga=50  %?
k=2.5*10^7
ks=1000  %??  
ks1=7.5*10^7
ks2=2.5*10^9
cs=20  %????
cs1=350
cs2=500
b=0.02  %??
tau=omiga*t
P=m2/2
delta=((miu*omiga*R*L)/P)*(R/c)^2*((L/2*R)^2)

pi=3.14159
%无量纲化
epus1=c1/(m1*omiga)
epus3=c1/(m1*omiga)
epus2=c2/(m2*omiga)
epus4=cs/(m3*omiga)
eta1=k/(m1*omiga^2)
eta3=k/(m1*omiga^2)
eta2=k/(m2*omiga^2)
eta4=ks/(m3*omiga^2)
M1=(m1*c*omiga^2)/(delta*P)
M3=(m1*c*omiga^2)/(delta*P)
M4=(m3*c*omiga^2)/(delta*P)
g=9.8;
G=g/(omiga^2)
%微分方程组程序
dx(1)=x(4)
dx(2)=x(5)
dx(3)=x(6)
dy(1)=y(5)
dy(2)=y(6)
dy(3)=y(7)
dy(4)=y(8)
%油膜力x,y方向的分量表达式
[fx1,fy1]=youmoli(x(1),y(1),dx(1),dy(1))
[fx2,fy2]=youmoli(x(3),y(3)-y(4),dx(3),dy(3)-dy(4))
%继续
dx(4)=-epus1*x(4)-eta1*(x(1)-x(2))+1/M1*fx1
dy(5)=-epus1*y(5)-eta1*(y(1)-y(2))+1/M1*fy1
dx(5)=-epus2*x(5)-eta2*(2*x(2)-x(1)-x(3))+b*cos(omiga*t)
dy(6)=-epus2*y(6)-eta2*(y(2)-y(1)-y(3))+b*sin(omiga*t)-G
dx(6)=-epus3*x(6)-eta3*(x(3)-x(2))+1/M3*fx2
dy(7)=-epus3*y(7)-eta3*(y(3)-y(2))+1/M3*fy2-G
dy(8)=-epus4*y(8)-eta4*y(4)-1/M4*fy2-G

运行后老是显示：
???  In an assignment  A(I) = B, the number of elements in B and
I must be the same.

Error in ==> sdzz1 at 60
dx(5)=-epus2*x(5)-eta2*(2*x(2)-x(1)-x(3))+b*cos(omiga*t)
请高手们帮我指点下迷津好不？？谢谢谢谢！！！！！

无水1324

最好给出youmoli函数，调试一下。
   dx(5)=-epus2*x(5)-eta2*(2*x(2)-x(1)-x(3))+b*cos(omiga*t)；
b*cos(omiga*t)在的t是一个向量，需要改为点乘

xiaojin831108

function [fx,fy]=youmoli(x,y,dx,dy)
%x=2;y=3;dx=3;dy=4;
alpha=atan((y+2*dx)/(x-2*dy))-pi/2*sign((y+2*dx)/(x-2*dy))-pi/2*sign(y+2*dx);
S=(x*cos(alpha)+y*sin(alpha))/(1-(x*cos(alpha)+y*sin(alpha)^2));
L=(2/(sqrt(1-x^2-y^2)))/(pi/2+atan(((y*cos(alpha)-x*sin(alpha))/sqrt(1-x^2-y^2))))
V=(2+(y*cos(alpha)-x*sin(alpha))*L)/(1-x^2-y^2)
fx=(sqrt((x-2*dy)+(y+2*dx)^2)/(1-x^2-y^2))*(3*x*V-sin(alpha)*L-2*cos(alpha)*S)
fy=(sqrt((x-2*dy)+(y+2*dx)^2)/(1-x^2-y^2))*(3*y*V-cos(alpha)*L-2*sin(alpha)*S)
这是youmoli函数，现在我又不知道用ode45能不能解耦合的强非线性的微分方程组了
这个是有关x,y的微分方程

xiaojin831108

比如已经有了松动转子系统的非微分方程组，我编了matlab的.m文件，然后在调用matlab的ode45来解这个方程组，感觉初值对其影响 很大，有时候要运行好久也没结果，像死循环，所有有哪个高手做过这个的，能不能告诉我这个初值怎么定义的啊？定的是多少

无水1324

ode45是可以求解耦合的强非线性的微分方程组的，
非线性方程对初值都很敏感，需要多次试，一般取0或者1

无水1324

你的问题就是出在b*cos(omiga*t)中；
t的取值是t=0:20;它是一个向量，而dx(5)是一个值，那么你要求时间在（0，20）中的解就必须一点一点求，增加循环：

t=0:20；
for i=1:length(t)
    ..................
dx(5)=-epus2*x(5)-eta2*(2*x(2)-x(1)-x(3))+b*cos(omiga*t(i));
    ..................

end

xiaojin831108

谢谢！方程可能没什么问题了，就是那个初值好难选啊，一直试的都不合适，提示精度有问题呵呵，不知道是不是选的不合适还是别的问题？

无水1324

你再多多试几个初值看一下。

xiaojin831108

楼主所说的初值一般是0或1吗？是不是初值只能试啊？有个问题是我试的时候每试一次程序都运行好久，不知道能不能出结果？一般松动转子系统运行时间一般是多少啊？

无水1324

只有试，
运行时间长是有可能的。
我原来经常遇到类似的问题，要运行几个小时，
选择初值时，要避免方程中出现无限大数或者分母中出现趋于0的数，

xiaojin831108

matlab有这样的提示：
Failure at t=2.620111e-008.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.293956e-023) at time t.
是不是就是楼主所说的方程中出现了无限大数或者分母中出现趋于0的数，以前楼主做过这个转子系统吗？我看有的说这个要是合适的话最多运行时间是580多秒.

无水1324

对
这种情况我也遇到过，计算到某一点的时候就进行不下去了，
我主要是做齿轮动力学方面的。

xiaojin831108

谢谢楼主！！！！
那这个通过试的方法太难了，因为有14个初值要定义，就用0或 1试啊？一般要是选的合适的话要运行多久啊？是不是运行时间太长就直接停下来算了啊？