马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
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 15:00
提升卡
变色卡
显身卡
楼主
发表于 2009-12-1 21:01
