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

tianya7071

下面的程序想用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);
    [t,q]=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)]

咕噜噜

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

无水1324

Newmark算法编辑好之后，是不是可以将ｏｄｅ４５换成Newmark函数了

tianya7071

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

无水1324

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

无水1324

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

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

tianya7071

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

无水1324

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

咕噜噜

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

tianya7071

方程是这样的
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

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

tianya7071

谢谢楼主~