satlxl 发表于 2007-5-16 17:03

我的程序得到的怎么不是正余弦的振动呢?

%用newmark法求解非惯性场中的振动问题
clc
clear
%定义柔性梁材料参数,l为柔性梁的长度,a为横截面积,ii为转动惯量
l=1.2;aa=0.0006;   ii=1.8e-007;
p=7866; ee=2.07*1e+11;r0=0.05;
n=2   %%%n是取的阶数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dt=0.1;
t_0=0;
t_max=4;
jisuanbushu=(t_max-t_0)/dt+1;
q1=zeros(n,1);
q2=zeros(n,1);
u=zeros(2*n,jisuanbushu);
v=zeros(2*n,jisuanbushu);
a=zeros(2*n,jisuanbushu);
%计算模态矩阵ph1为纵向变形的,ph2为横向变形的         
syms x t real
for i=1:n
    ph1(i)=sin((2*i-1)*pi*x/(2*l));
    if i==1
      ph2(i)=cos(1.875*x/l)-cosh(1.875*x/l)-(cos(1.875)+cosh(1.875))*(sin(1.875*x/l)-sinh(1.875*x/l))/(sin(1.875)+sinh(1.875));
    elseif i==2
      ph2(i)=cos(4.694*x/l)-cosh(4.694*x/l)-(cos(4.694)+cosh(4.694))*(sin(4.694*x/l)-sinh(4.694*x/l))/(sin(4.694)+sinh(4.694));
    else
      ph2(i)=cos((i-0.5)*pi*x/l)-cosh((i-0.5)*pi*x/l)-(cos((i-0.5)*pi)+cosh((i-0.5)*pi))*(sin((i-0.5)*pi*x/l)-sinh((i-0.5)*pi*x/l))/(sin((i-0.5)*pi)+sinh((i-0.5)*pi));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k1=int(ee*aa*(diff(ph1)).'*(diff(ph1)),x,0,l);
k1=eval(k1);
lxl1=(diff(ph2,2)).'*(diff(ph2,2));
lxl2=(1-x)*(diff(ph2)).'*(diff(ph2));
lxl3=ph1.'*ph2;
lxl4=ph1.'*ph1;
lxl5=ph2'*ph2;
for i=1:n
    for j=1:n      
      k2(i,j)=ee*ii*quadl(inline(lxl1(i,j)),0,l);
      d0(i,j)=p*aa*quadl(inline(lxl2(i,j)),0,l);
      d1(i,j)=p*aa*quadl(inline(x*lxl2(i,j)),0,l);
      rr(i,j)=p*aa*quadl(inline(lxl3(i,j)),0,l);
      m1(i,j)=p*aa*quadl(inline(lxl4(i,j)),0,l);
      m2(i,j)=p*aa*quadl(inline(lxl5(i,j)),0,l);
    end
end
u01=int(p*aa*ph1,0,l);
u02=int(p*aa*ph2,0,l);
u11=int(p*aa*x*ph1,0,l);
u12=int(p*aa*x*ph2,0,l);
gq1q2=-rr;
gq2q1=rr';
%%%%%%%%%%%%%%%%%%%%%%%%%
t=t_0-dt;
for i=1:1:jisuanbushu-1
    t=t+dt   
    omega=49.6*(1-exp(-t));
    domega=49.6*exp(-t);
    if t>=4
      break
    end
    %计算各系数矩阵的值
    mq1c=-rr*q2;
    mq2c=(r0*u02+u12+q1'*rr)';
    kq1q1=k1-(omega)^2*m1;
    kq2q2=k2-(omega)*m2+(omega)^2*(r0*d0+d1);
    qq1=(omega)^2*(r0*u01'+u11');
   
    M=;
    C=2*omega*;
    K=;   
    Q=;
    Q=eval(Q);
   

    %利用newmark法求解梁的振动位移
    b=0.25;
    r=0.5;

    a0=1/(b*(dt)^2);
    a1=r/(b*dt);
    a2=1/(b*dt);
    a3=1/(2*b)-1;
    a4=r/b-1;
    a5=0.5*dt*((r/b)-2);
    a6=dt*(1-r);
    a7=r*dt;

    K=K+a0*M+a1*C;
    Q=Q+M*(a6*u(:,i)+a2*v(:,i)+a3*a(:,i))+C*(a1*u(:,i)+a4*v(:,i)+a5*a(:,i));
    u(:,i+1)=inv(K)*Q;
    a(:,i+1)=a0*(u(:,i+1)-u(:,i))-a2*v(:,i)-a3*a(:,i);
    v(:,i+1)=v(:,i)+a6*a(:,i)+a7*a(:,i+1);
    for j=1:n
      q1(j,1)=u(j,i+1);
      q1(j,1)=u(n+j,i+1);
    end
end
%画出梁末端位移响应图**************************************************

T=t_0:dt:t_max;
figure(1)            %%纵向振动位移
plot(T,u(1,:))
title('柔性叶片末端纵向振动位移图')
xlabel('时间(t)')
ylabel('位移(m)')
figure(2)         %%横向振动位移
plot(T,u(n+1,:))
title('柔性叶片末端横向振动位移图')
xlabel('时间(t)')
ylabel('位移(m)')

得到的图形是一条直线,最后特别大。如果积分步长特别小,得到的就对

无水1324 发表于 2007-5-19 21:57

本帖最后由 VibInfo 于 2016-5-11 16:03 编辑

原帖由 satlxl 于 2007-5-16 17:03 发表
%用newmark法求解非惯性场中的振动问题
clc
clear
%定义柔性梁材料参数,l为柔性梁的长度,a为横截面积,ii为转动惯量
l=1.2;aa=0.0006;   ii=1.8e-007;
p=7866; ee=2.07*1e+11;r0=0.05;
n=2   %%%   ...

步长小,就对这个怎么解释?

satlxl 发表于 2007-5-20 19:58

谢谢你的意见,我试了一下,步长取0.05就发散了,原来得到的直线的原因是后面突然间数字变的很大,
前面的就约等于零了,怎么回事呢??
页: [1]
查看完整版本: 我的程序得到的怎么不是正余弦的振动呢?