请教newmark迭代求动力响应
最近在做车桥耦合的论文，计算用自己编的程序，因为把车和桥看成两个独立的子系统，在最后求动力响应的时候要用耦合条件，通过迭代的算法来实现，这部分程序写出来了，可是计算出每一时刻的位移都不收敛，造成死循环，请高手帮忙解决一下吧，我自己查了一个礼拜了，还是没头绪，不知道哪里错了！跪求解决办法，毕业在即，论文里就剩下这点东西了，拜托高手帮帮忙吧！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;
%-------------------------计算车辆和桥梁的刚度矩阵、质量矩阵、阻尼矩阵
=vehicle;
=bridge;
%=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
=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 把其余子程序一起放上来才好调试！ =vehicle;是求解车辆刚度矩阵、质量矩阵、阻尼矩阵和车轮的刚度矩阵、阻尼矩阵，W为每个车轮承担的重力
=bridge;Kb，Mb，Cb分别为桥梁单元的刚度矩阵，质量矩阵和阻尼矩阵，LOC为单元的节点编号，按列排列
=roadsurface(x,16);求路面不平度，和此刻路面不平度的变化率
这些子程序都ok没有任何问题，就最后这个有问题，因为每个时刻的动力响应都是假设已知的，然后通过迭代计算，但是每次迭代时候，这次计算的位移和上次的位移，误差总是1，所以无限循环下去，麻烦哪位大侠能帮帮忙，万分感谢啊！ 个人水平/时间有限, 建议楼主想法简化下!
好像很复杂, 不太敢看!:@L 不收敛 一般是参数选择不合适 当然首先要确定程序正确 看来只能简化再简化了，别的好办法我也没了，除非用ansys做，不过那样论文没程序支持，发布到好的地方 2bu2 发表于 2009-12-2 13:23 static/image/common/back.gif
看来只能简化再简化了，别的好办法我也没了，除非用ansys做，不过那样论文没程序支持，发布到好的地方
请问您问题解决了吗？？ 楼主，我现在遇到同样的问题了，能否向您请教！，没有头绪 非线性迭代最好是把步长搞小一点，要不难收敛，才0.02，改成0.00002慢慢等结果
页:
[1]