线性系统为什么会有这样的相图？

siyanger

做了一个三自由度的线性系统，用ode45解出来却有这样的相图和响应。又是刚度矩阵值的问题吗？
w=0.2;
T=2*pi/w;
%tspan=[1000,1290];%1.0e-003 *
tspan=[0:T/400:50*T]
X0=[ -0.0089;-0.7717; -0.0237];
Y0=[X0;0;0;0];
[t,y]=ode45('cdyanalysis2',tspan,Y0);


function fungear=cdylanalysis2(t,y)
w=0.2;
E=0.2;
FVector1 = 0.5+E*w^2*cos(w*t);
FVector3 = 0.48+E*w^2*cos(w*t);
            
K =[1.0000    0.005   0;     1.42      0.09    2.05;
    0         0.03    1.08];

K11=K(1,1); K12=K(1,2); K21=K(2,1);
K22=K(2,2); K23=K(2,3);K32=K(3,2); K33=K(3,3);

C  =[ 0.0121    0    0;
         0.0121    0    0.012;
          0    0    0.012];
      
C11=C(1,1);C21=C(2,1);C23=C(2,3);C33=C(3,3);

fungear=[y(4);
          y(5);
          y(6);
           -C11*y(4)-K11*y(1)+ K12*y(2)+FVector1;
           C21*y(4)+C23*y(6)+K21*y(1)-K22*y(2)+K23*y(3);
           -C33*y(6)+K32*y(2)-K33*y(3)+FVector3];

siyanger

如下图：

siyanger

将时间各往后取，响应图还是那样，好象永远得不到稳态解似的。
相图还是一样的。

siyanger

用下面的程序画出的图就更不对了。是程序有问题还是系统有问题呢？
%　 poincare 图      
zz=400;
x1=y(:,1);
dx1=y(:,4);
xx1=x1(20:zz:n2);
dxx1=dx1(20:zz:n2);         
            figure('unit','normalized','color',[1,1,1]);
            H_1 = get(gcf);
            set(gcf,'Name','poincare 图','numbertitle','off');
            set(gcf,'defaultuicontrolunits','normalized');
            set(gcf,'defaultuicontrolfontname','仿宋体');
            set(gcf,'defaultuicontrolhorizontalAlignment','left');
            plot_1 = plot(xx1, dxx1,'k*');
            axis([-2 2  -2.5 2.5])
            grid off
% 功率谱           
Y=fft(x1);
Y(1)=[];
n=length(Y);
Power=abs(Y(1:n2/2)).^2/n2^2;
freq=400/T*(1:n2/2)./n;
figure
plot(freq,Power,'k')
axis([0 2  0 0.05])
title('功率谱');
xlabel('频率');ylabel('位移')

toes

有意思。
没有仔细看你的系统模型，试着算了一下，在我试算的时间范围内确实没有得到平稳的结果。后来发现阻尼好像比较小，你可以试试，如果阻尼系数增大到原来的10倍，就比较容易得到平稳结果了。
不知道阻尼是否取的合适？阻尼比是多少呢？

个人意见，仅供参考。

[ 本帖最后由 toes 于 2006-9-8 19:09 编辑 ]

siyanger

谢谢toes！我将上述阻尼比增大，得到一个发散的结果。阻尼比原来取的为０．０３。但因为一直得不到收敛的结果，我将刚度和阻尼矩阵都改了，没按原来的系统结构参数改。
我将阻尼改为：
   
C  =[ 0.12    0    0;
         0.12     0    0.12;
         0　　   0    0.12];      
怎么还是不行。

toes

估计是你没有完全消除初值的影响。
由于初值的影响，系统的响应在开始一段时间会出现以固有频率为频率的振荡。为了消除初值的影响，需要将前面一段响应去掉，仅分析后面稳定的响应。
在你原来的基础上，我将阻尼系数全部乘以10。计算400个周期，取后一半数据分析，得到的结果如下图：

400t

如果计算800个周期，取后一半数据分析，得到的结果如下图：

800t

siyanger

非常感谢！我按照你说的方法做出来和你一样的相图了，下面我想做一个多自由度的间隙非线性系统。
发现能在论坛上和大家交流真是收获很大，比找同学和老师讨论要方便多啦！

toes

原帖由 siyanger 于 2006-9-9 10:00 发表
非常感谢！我按照你说的方法做出来和你一样的相图了，下面我想做一个多自由度的间隙非线性系统。
发现能在论坛上和大家交流真是收获很大，比找同学和老师讨论要方便多啦！

多自由度的间隙非线性系统，如果做出来了给大家鉴赏一下。

glise

原帖由 siyanger 于 2006-9-8 22:00 发表
谢谢toes！我将上述阻尼比增大，得到一个发散的结果。阻尼比原来取的为０．０３。但因为一直得不到收敛的结果，我将刚度和阻尼矩阵都改了，没按原来的系统结构参数改。
我将阻尼改为：
   
C  =[ 0.12    0    ...

这个阻尼比给得好像偏大了，实际系统很少有这么大的

toes

阻尼可能是偏大的，所以最好算一下阻尼比，我们并不清楚结构的实际情况，阻尼最好还是楼主结合实际给一下。
我算了一下，阻尼小的情况下系统很难得到稳定的解，最大我算到了10000个周期，再算下去我的电脑负担有点重，呵呵。
如果不改阻尼大家有什么好的建议求解该题？谢谢。

siyanger

原帖由 toes 于 2006-9-9 16:18 发表


多自由度的间隙非线性系统，如果做出来了给大家鉴赏一下。

好啊，我想做些参数对沌混啊，分叉什么的影响方面分析，但前提是得做出来。

siyanger

原帖由 toes 于 2006-9-9 20:24 发表
阻尼可能是偏大的，所以最好算一下阻尼比，我们并不清楚结构的实际情况，阻尼最好还是楼主结合实际给一下。
我算了一下，阻尼小的情况下系统很难得到稳定的解，最大我算到了10000个周期，再算下去我的电脑负担有 ...

这个阻尼是比一般的系统要大。是齿轮的啮合阻尼。阻尼系数为：C=2*Ksi*sqrt(K/(1/m1+1/m2))，Ksi为阻尼比，0.03-0.17。以上的公式参考　李润方，王建军的　《齿轮系统动力学》Ｐ１６１。
一般啮合刚度为１０^8左右，而这个系统中质量为0.2kg左右吧，这样在数学模型中Ｃ前面还有系数，最后算来阻尼矩阵中最大的元素为０.５８，最小的为０.１。（化为无量纲后）
以前用这个阻尼系数用近似解析法算的时候，结果还可以的。但用数值法解还是未知数呢！

siyanger

原帖由 toes 于 2006-9-9 20:24 发表
阻尼可能是偏大的，所以最好算一下阻尼比，我们并不清楚结构的实际情况，阻尼最好还是楼主结合实际给一下。
我算了一下，阻尼小的情况下系统很难得到稳定的解，最大我算到了10000个周期，再算下去我的电脑负担有 ...

这个阻尼是比一般的系统要大。是齿轮的啮合阻尼。阻尼系数为：C=2*Ksi*sqrt(K/(1/m1+1/m2))，Ksi为阻尼比，0.03-0.17。以上的公式参考　李润方，王建军的　《齿轮系统动力学》Ｐ１６１。
一般啮合刚度为１０^8左右，而这个系统中质量为0.2kg左右吧，这样在数学模型中Ｃ前面还有系数，最后算来阻尼矩阵中最大的元素为０.５８，最小的为０.１。（化为无量纲后）
以前用这个阻尼系数用近似解析法算的时候，结果还可以的。但用数值法解还是未知数呢！

wendyrong

我是新手，为什么我把第一个程序导入后出现这句呢：
Error: File: TRY3.m Line: 9 Column: 1
Function definitions are not permitted at the prompt or in scripts.