声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4306|回复: 8

[编程技巧] 请教newmark迭代求动力响应

[复制链接]
发表于 2009-12-1 15:00 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
最近在做车桥耦合的论文,计算用自己编的程序,因为把车和桥看成两个独立的子系统,在最后求动力响应的时候要用耦合条件,通过迭代的算法来实现,这部分程序写出来了,可是计算出每一时刻的位移都不收敛,造成死循环,请高手帮忙解决一下吧,我自己查了一个礼拜了,还是没头绪,不知道哪里错了!跪求解决办法,毕业在即,论文里就剩下这点东西了,拜托高手帮帮忙吧!


function coupling
format short e
%
%-------------------------选择时间步长Δt、参数β和γ,计算积分常数
beta=0.25;                    %β=0.25
gamma=0.5;                    %γ=0.5
dt=0.2;                       %间步长Δt=0.2(秒)
a0=1/(beta*dt^2);
a1=gamma/(beta*dt);
a2=1/(beta*dt);
a3=1/(2*beta)-1;
a4=gamma/beta-1;
a5=0.5*dt*(gamma/beta-2);
a6=(1-gamma)*dt;
a7=gamma*dt;
%-------------------------计算车辆和桥梁的刚度矩阵、质量矩阵、阻尼矩阵
[Kv,Mv,Cv,Kt,Ct,W]=vehicle;
[Kb,Mb,Cb,LOC]=bridge;
%[Kb,Mb,Cb,LOC]=bridgeplane;
%-------------------------计算车辆和桥梁的等效刚度矩阵Kb和Kv
Kb_=Kb+a0*Mb+a1*Cb;
Kv_=Kv+a0*Mv+a1*Cv;
%-------------------------取初始值
Ub=ones(size(Kb,1),1);               %桥梁初始位移
Ub1=ones(size(Kb,1),1);              %桥梁初始速度
Ub2=ones(size(Kb,1),1);              %桥梁初始加速度
Uv=ones(size(Kv,1),1);               %车辆初始位移
Uv1=ones(size(Kv,1),1);              %车辆初始速度
Uv2=ones(size(Kv,1),1);              %车辆初始加速度
ub(:,1)=Ub;
ub1(:,1)=Ub1;
ub2(:,1)=Ub2;
uv(:,1)=Uv;
uv1(:,1)=Uv1;
uv2(:,1)=Uv2;
L=4;%桥的长度
v=20*1000/3600;        %车速km/h
x=2;

i=1;
t(1)=0;
I=1;
%t(i)*v<=L
%--------------------开始一个newmark时间步的迭代----------------
while I<=2;
    k=0;%迭代次数
    delta=0.03;
    %--------------计算中间变量
    Sb=Mb*(a0*ub(:,i)+a2*ub1(:,i)+a3*ub2(:,i))+Cb*(a1*ub(:,i)+a4*ub1(:,i)+a5*ub2(:,i));
    Sv=Mv*(a0*uv(:,i)+a2*uv1(:,i)+a3*uv2(:,i))+Cv*(a1*uv(:,i)+a4*uv1(:,i)+a5*uv2(:,i));
    %--------------计算路面不平度
    for j=1:2
        [r,r1]=roadsurface(x,16);
        rb(j)=r;
        rb1(j)=r1;
        x=x+v*dt;         
    end
   
    %-----------------------开始迭代
    while delta>0.02             %判断是否满足精度要求
          %-------------桥梁在t+Δt时刻的加速度、速度
          if(k==0)
              ub2(:,i+1)=a0*(ub(:,i)-ub(:,i))-a2*ub1(:,i)-a3*ub2(:,i)
              ub1(:,i+1)=ub1(:,i)+a6*ub2(:,i)+a7*ub2(:,i+1)
              
          else
              ub2(:,i+1)=a0*(ub(:,i+1)-ub(:,i))-a2*ub1(:,i)-a3*ub2(:,i);
              ub1(:,i+1)=ub1(:,i)+a6*ub2(:,i)+a7*ub2(:,i+1);
              
          end
   
          %--------------形成车辆方程右端荷载项
          %za=zeros(4,1);
          za(1)=ub(3*LOC(6,I),i)+rb(2);
          za(2)=ub(3*LOC(7,I),i)+rb(2);
          za(3)=ub(3*LOC(5,I),i)+rb(1);
          za(4)=ub(3*LOC(8,I),i)+rb(1);
          za=reshape(za,4,1);
          %za1=zeros(4,1);
          za1(1)=ub1(3*LOC(6,I),i)+rb1(2);
          za1(2)=ub1(3*LOC(7,I),i)+rb1(2);
          za1(3)=ub1(3*LOC(5,I),i)+rb1(1);
          za1(4)=ub1(3*LOC(8,I),i)+rb1(1);
          za1=reshape(za1,4,1);
         
          %----------------车辆等效荷载列阵
          Pv=Kt*za+Ct*za1;
          Pv_=Pv+Sv;               
         
          %-------------求解t+Δt时刻的车辆位移列阵
         
          uv(:,i+1)=inv(Kv_)*Pv_;
         
          %-------------车辆在t+Δt时刻的加速度、速度
          uv2(:,i+1)=a0*(uv(:,i+1)-uv(:,i))-a2*uv1(:,i)-a3*uv2(:,i);
          uv1(:,i+1)=uv1(:,i)+a6*uv2(:,i)+a7*uv2(:,i+1);
          %-------------计算桥梁方程右端荷载项
          Pb=zeros(size(Kb,1),1);
          for j=i:4
              Pb0(j)=W(j)+Kt(3+j,j)*(uv(j+3,i+1)-za(j,1))+Ct(3+j,j)*(uv1(j+3,i+1)-za1(j,1));
          end
          Pb(3*LOC(6,I),1)=Pb0(1);
          Pb(3*LOC(7,I),1)=Pb0(2);
          Pb(3*LOC(5,I),1)=Pb0(3);
          Pb(3*LOC(8,I),1)=Pb0(4);
          Pb_=Pb+Sb;               %桥梁等效荷载列阵
         
          %--------------求解t+Δt时刻的桥梁位移列阵
          ub(:,i+1)=inv(Kb_)*Pb_;
         
         
          %--------------桥梁在t+Δt时刻的加速度、速度
          ub2(:,i+1)=a0*(ub(:,i+1)-ub(:,i))-a2*ub1(:,i)-a3*ub2(:,i);
          ub1(:,i+1)=ub1(:,i)+a6*ub2(:,i)+a7*ub2(:,i+1);
         
          %---------------------判断桥梁位移是否收敛-----------------
          u1=ub(:,i);
          u2=ub(:,i+1);
         
          delta=norm((u2-u1)/(u2+eps))
          k=k+1
    end
   
    %--------------------进入下一时刻


    i=i+1;
    t(i)=t(i-1)+dt;
    x=v*dt;
    I=I+1;
end

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2009-12-1 19:22 | 显示全部楼层
把其余子程序一起放上来才好调试!
 楼主| 发表于 2009-12-1 20:35 | 显示全部楼层
[Kv,Mv,Cv,Kt,Ct,W]=vehicle;是求解车辆刚度矩阵、质量矩阵、阻尼矩阵和车轮的刚度矩阵、阻尼矩阵,W为每个车轮承担的重力
[Kb,Mb,Cb,LOC]=bridge;Kb,Mb,Cb分别为桥梁单元的刚度矩阵,质量矩阵和阻尼矩阵,LOC为单元的节点编号,按列排列
[r,r1]=roadsurface(x,16);求路面不平度,和此刻路面不平度的变化率
这些子程序都ok没有任何问题,就最后这个有问题,因为每个时刻的动力响应都是假设已知的,然后通过迭代计算,但是每次迭代时候,这次计算的位移和上次的位移,误差总是1,所以无限循环下去,麻烦哪位大侠能帮帮忙,万分感谢啊!
发表于 2009-12-1 21:01 | 显示全部楼层
个人水平/时间有限, 建议楼主想法简化下!
好像很复杂, 不太敢看!:@L
发表于 2009-12-2 08:44 | 显示全部楼层
不收敛 一般是参数选择不合适 当然首先要确定程序正确
 楼主| 发表于 2009-12-2 13:23 | 显示全部楼层
看来只能简化再简化了,别的好办法我也没了,除非用ansys做,不过那样论文没程序支持,发布到好的地方
发表于 2013-10-16 23:02 | 显示全部楼层

请问您问题解决了吗??
发表于 2018-4-16 19:39 | 显示全部楼层
楼主,我现在遇到同样的问题了,能否向您请教!,没有头绪
发表于 2023-2-3 15:30 来自手机 | 显示全部楼层
非线性迭代最好是把步长搞小一点,要不难收敛,才0.02,改成0.00002慢慢等结果
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-3-28 19:56 , Processed in 0.109190 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表