tianya7071 发表于 2007-8-23 16:37

帮我看看啊!下面这个程序怎么用Newmark法来求解

下面的程序想用Newmark方法来求解,不知道怎么弄,请大家指点~
主程序
function Rotors_System_Func
clear
clc
global BN Nb1 Nb2 w1 w2 Ro1 Ri1 Ro2 Ri2
n_one_T=1;   
n_T=1;
%e217QT(85*15*28)
Ro1=48.828;
Ri1=68.672;
Nb1=10;
%6D2272822(110*140*19)
Ro2=66.515;
Ri2=58.5;
Nb2=34;
n_n=0;
w2_min=100;
w2_max=200;
w2_step=100;
w1_rpm=11000;
q_initial(1:8,1)=1e-11;
BN=Ri1/(Ri1+Ro1)*Nb1;
tic
w2_rpm=10400;
%for w2_rpm=w2_min:w2_step:w2_max
    n_n=n_n+1
    w2=w2_rpm/60;
    w1=w1_rpm/60;
    T_vc=2*pi/w2;
    dt=T_vc/n_one_T;
    time=n_T*T_vc;
    n=round(time/dt);
    t_span(1:n)=linspace(0,time,n);
    %options=odeset('RelTol',1e-5);
    =ode45('Rotors_System_Sub_Func',t_span,q_initial);
%end
toc
%subplot(2,1,1);plot(w2_rpm,q(:,1),'k.')
%subplot(2,1,2);
plot(q(8000:100:10000,1),q(8000:100:10000,2))
这是子程序:
function dq=Rotors_System_Sub_Func(t,q)
global BN Nb1 Nb2 w1 w2 Ro1 Ri1 Ro2 Ri2
r01=2e-2;
r02=5e-2;
m1=7.86;
m2=11.93;
W1=m1*9.8;
W2=m2*9.8;
Fx11=0;
Fy11=0;
%Deformation=zeros(10,1);
for i=1:Nb1
    Deformation=zeros(Nb1,1);
    sita(i)=2*pi/Nb1*(i-1)+BN/Nb1*w1*2*pi*t;
    Deformation(i,1)=q(1,1)*cos(sita(i))+q(2,1)*sin(sita(i))-r01;
    if Deformation(i)<=0
      Deformation(i)=0;
    end
    Kb1=1.4656e6;
    Kb2=1.7518e6;
    fx11=Kb1*Deformation(i)^1.5*cos(sita(i));
    fy11=Kb1*Deformation(i)^1.5*sin(sita(i));
   
    Fx11=Fx11+fx11;
    Fy11=Fy11+fy11;
end
Fx12=Fx11;
Fy12=Fy11;
Fx21=0;
Fy21=0;
for j=1:Nb2
    sita(j)=2*pi/Nb2*(j-1)+(Ro2/(Ro2+Ri2)*w1+Ri2/(Ro2+Ri2)*w2)*2*pi*t;
    Deformation1(j,1)=q(3,1)*cos(sita(j))+q(4,1)*sin(sita(j))-r02;
    if Deformation1(j)<=0
      Deformation1(j)=0;
    end
   
    fx21=Kb2*Deformation1(j)^1.1*cos(sita(i));
    fy21=Kb2*Deformation1(j)^1.15*sin(sita(i));
   
Fx21=Fx21+fx21;
Fy21=Fy21+fy21;
end
Fx22=0;
Fy22=0;
for k=1:Nb1
    sita(k)=2*pi/Nb1*(k-1)+BN/Nb1*w2*2*pi*t;
    Deformation2(k,1)=q(3,1)*cos(sita(i))+q(4,1)*sin(sita(i))-r01;
    if Deformation2(k)<=0
      Deformation2(k)=0;
    end
   
    fx22=Kb1*Deformation2(k)^1.5*cos(sita(k));
    fy22=Kb1*Deformation2(k)^1.5*sin(sita(k));
   
    Fx22=Fx22+fx22;
    Fy22=Fy22+fy22;
end
P=1470;
Cx11=1.500;Cy11=1.500;
Cx12=2.500;Cy12=2.500;
Cx21=2.000;Cy21=2.000;
Cx22=7.000;Cy22=7.000;
dq(1:4,1)=q(5:8,1);
dq(5:8,1)=1000*[
          -1/m1*((Cx11+Cx12)*q(5,1)+Fx11+Fx12-(W1+Fx21+Cx21*q(7,1)));...
          -1/m1*((Cy11+Cy12)*q(6,1)+Fy11+Fy12-(Fy21+Cy21*q(8,1)));
          -1/m2*((Cx21+Cx22)*q(7,1)+Fx21+Fx22-(P+W2));
          -1/m2*((Cy21+Cy22)*q(8,1)+Fy21+Fy22)]

咕噜噜 发表于 2007-8-23 16:57

回复 #1 tianya7071 的帖子

:@L :@L 程序好长啊
你能不能把你的问题直接说出来,你原来用什么方法,现在要用Newmark法?

无水1324 发表于 2007-8-23 19:23

回复 #1 tianya7071 的帖子

Newmark算法编辑好之后,是不是可以将ode45换成Newmark函数了

tianya7071 发表于 2007-8-24 08:39

由于是用有限元列的振动方程,质量阵,刚度阵比较大用ode45求解比较麻烦,在网上看到有用Newmark求响应的,用ode45求响应时需要用的主程序和子程序两个m文件,用Newmark求时是不是只要一个m函数文件就可以了 ,程序中的循环比较多不知道怎么求解
请大家讨论~~

无水1324 发表于 2007-8-24 09:13

回复 #4 tianya7071 的帖子

具体的用Newmark计算没有做过,看你这个Rotors_System_Sub_Func函数确实比叫复杂的,里面的参数受很多限制,在用Newmark的时候也必须一样,循环是不能避免的。

无水1324 发表于 2007-8-24 09:21

回复 #4 tianya7071 的帖子

你这个转子碰摩系统的方程求解吧。应该用ode45求解比较快。我以前遇到过几个人都是用ode45求解该类方程的。

用Newmark你还是找个例子学习一下,因为具体的我也没有做过。

tianya7071 发表于 2007-8-24 09:44

我做的是非线性振动微分方程,是由于力的非线性引起的问题,系统自由度比较多,而且非线性力和响应有关,感到用响应求非线性力,然后再用非线性力求响应,这样反复的循环不知道怎么弄才好

无水1324 发表于 2007-8-25 13:19

回复 #7 tianya7071 的帖子

非线性力和响应?这个没有看懂,应该是与激励有关系吧

咕噜噜 发表于 2007-8-25 14:52

回复 #7 tianya7071 的帖子

非线性力和响应有关?什么意思?应该是非线性力是位移响应的非线性函数吧?

tianya7071 发表于 2007-8-26 21:07

方程是这样的
Mx 十 Cx + Kx = g + f(x,z) 其中 ,M,C ,K ER"",x ER",g ER"是与时间有关的力矢量,
f(x,z)是非线性力矢,M后面乘的是加速度向量,C后面乘的是速度向量。
f(x,z)用for i=1:Nb1
    Deformation=zeros(Nb1,1);
    sita(i)=2*pi/Nb1*(i-1)+BN/Nb1*w1*2*pi*t;
    Deformation(i,1)=xcos(sita(i))+z*sin(sita(i))-r01;
    if Deformation(i)<=0
      Deformation(i)=0;
    end
    Kb1=1.4656e6;
    Kb2=1.7518e6;
    fx11=Kb1*Deformation(i)^1.5*cos(sita(i));
    fy11=Kb1*Deformation(i)^1.5*sin(sita(i));
   
    Fx11=Fx11+fx11;
    Fy11=Fy11+fy11;
end
Fx12=Fx11;
Fy12=Fy11
来迭代
请教如何用Newmark法来求解
欢迎大家讨论!

zhaowu 发表于 2007-8-31 20:52

看看http://forum.vibunion.com/forum/thread-50906-1-1.html

tianya7071 发表于 2007-9-4 16:52

谢谢楼主~
页: [1]
查看完整版本: 帮我看看啊!下面这个程序怎么用Newmark法来求解