newmark法求解非线性随机振动

大小不是白

大牛的有非线性系统结果比线性系统小，但我的比线性系统大，并且加大激励后非线性系统出现NaN，而线性系统结果正常。
求各位大神指教，已经纠结很久了。
matlab程序如下
N=30;  %%%%% 模态数量
L=20;
t_max=940;  %%%%%%%%%%%!!!!!!!!!!!!!!!!!!!!!!!  < tstep*numtstep
dt=0.5*1e-2;    %步长
E1=244.27.*(10.^9);
E2=122.56.*(10.^9);
M1=4429;
M2=3000;
V=0.3;
n=0;
h=0.6;
b=1;
A11=(E2./E1+n)./(n+1).*E1.*h./(1-(V.^2));
B11=n.*(E2./E1-1)./2./(n+1)./(n+2).*(h.^2).*E1./(1-(V.^2));
D11=((n.^3)+3.*(n.^2)+8.*n+3.*E2./E1.*(n.^2+n+2))./12./(n+1)./(n+2)./(n+3).*(h.^3).*E1./(1-(V.^2));
i1=(M2./M1+n)./(n+1).*M1.*h;  %%%%%%%  I1
A=D11-(B11.^2)./A11; %%%%%  
c=500;   %%%%%%  damping



M=eye(N);
C=c/i1*eye(N);
B=zeros(N);
for i=1:N
B(i,i)=(i*pi)^4*A/(i1*L^4);   %%%%%%%%%%%%%%%%%%%%4  
end
K=B;

%%%%%%%%% 非线性项
KK=zeros(N);
for i=1:N   
        KK(i,i)=(i*pi/L)^2*L;
end
RR=zeros(N);
for i=1:N   
        RR(i,i)=-(i*pi/L)^2*L;
end


%%%%%%%% random excitation
disq=zeros(N,1);
for i=1:N   
    kkkk=2*i-1;
    if kkkk<=N
    aa=2*sqrt(2)/pi/kkkk/i1;
    disq(kkkk,1)=aa;
    end
end    %%%%%%%  加权载荷（特征函数）


f=wt11e5(1,2)*disq;    %%%%%%  初始载荷值 t=0



u=zeros(N,1);%%%%[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]';
v=zeros(N,1);%%%[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]';

t(1)=0;               %初始时间
x(:,1)=u;             %初始位移
x1(:,1)=v;            %初始速度
a=M\f;
x2(:,1)=a;            %初始加速度

%%%%%%%%%%%newmark参数
gama=0.5;
beta=0.25;
a0=1/(beta*dt^2);
a1=gama/(beta*dt);
a2=1/(beta*dt);
a3=1/(2*beta)-1;
a4=gama/beta-1;
a5=dt*(gama/(2*beta)-1);
a6=dt*(1-gama);
a7=gama*dt;
%%%%%%%%%%%newmark参数


K1=K+a0*M+a1*C;      
i=2;
ll=t_max/dt;
for i=2:ll
      if mod(i,1000)==0
          outputa=i*dt;
      end
      t(i)=(i-1)*dt;
      Q=wt11e5(i,2)*disq;
      coef=A11/2/L/L/i1;      
      temp=coef*(x(:,i-1))'*KK*x(:,i-1);
      norn=1;iterx1=x(:,1);
      while (norn>0.001)
           q(:,i)=Q+M*(a0*x(:,i-1)+a2*x1(:,i-1)+a3*x2(:,i-1))+C*(a1*x(:,i-1)+a4*x1(:,i-1)+a5*x2(:,i-1));     
           iterx2=(K1-temp*RR)\q(:,i);
           diff=abs(iterx2-iterx1);
           norn=max(diff./abs(iterx1));
           iterx1=iterx2;
      end      
      x(:,i)=iterx1;  
      x2(:,i)=a0*(x(:,i)-x(:,i-1))-a2*x1(:,i-1)-a3*x2(:,i-1);  
       x1(:,i)=x1(:,i-1)+a6*x2(:,i-1)+a7*x2(:,i);
end  
plot(t,x(1,:))

心灯

没有文献信息，大家不了解背景，没法了解你的程序，所以注定你得不到想要的答复…

truleeee

??? Undefined function or method 'wt11e5' for input arguments of type 'double'.

Error in ==> newmak at 52
f=wt11e5(1,2)*disq;    %%%%%%  初始载荷值 t=0

运行报错了

大小不是白

不好意思，是我忽略了。@心灯 @truleeee 这是参考含分数阶梁的非线性随机振动做的功能梯度梁的随机振动，具体文献如下，wt11e5是激励谱

缱绻

你程序里也没调入txt文件啊

dollfish000

有没有Newmark解其它非线性问题的程序？

dollfish000

帮你出图了，那么久也没见来交流。

chayana

dollfish000 发表于 2017-11-26 14:59
帮你出图了，那么久也没见来交流。

怎么运行成功的？能分享一下程序吗？