请教一个纽马克程序问题
下面是用纽马克法求振动响应的matlab代码, 初始条件u,v,a都为零时，结果良好可是当v=;时，得出的u值特别大，与初始条件u,v,a都为零时相差甚远，这是什么原因？
请教论坛里的朋友了！！！
clc;clear;
j1=76318.474; j2=5255.968; j3=1063.431;
j4=7858.041; j5=821.753; j6=2290.285; j7=11342.453;
k1=6.464e8; k2=6.486e8; k3=5.929e8;
k4=1.981e8; k5=1.980e8; k6=11.595e8; M0=2000;
J=[j1 0 0 0 0 0 0; 0 j2 0 0 0 0 0; 0 0 j3 0 0 0 0; 0 0 0 j4 0 0 0; ...
0 0 0 0 j5 0 0; 0 0 0 0 0 j6 0; 0 0 0 0 0 0 j7]; %% J=diag(); %%by ChaChing
K=[k1 -k1 0 0 0 0 0; -k1 k1+k2 -k2 0 0 0 0; 0 -k2 k2+k3 -k3 0 0 0
0 0 -k3 k3+k4 -k4 0 0; 0 0 0 -k4 k4+k5 -k5 0; 0 0 0 0 -k5 k5+k6 -k6; 0 0 0 0 0 -k6 k6];
degree=length(J); C=zeros(degree,degree);
u(:,1)=zeros(degree,1); v(:,1)=; a(:,1)=zeros(degree,1);
n=10;%时间段数
tend=0.05; dt=tend/n; gaJa=0.5; beta=0.25;
alpha0=1/beta/dt^2; alpha1=gaJa/beta/dt; alpha2=1/beta/dt; alpha3=1/2/beta-1;
alpha4=gaJa/beta-1; alpha5=dt/2*(gaJa/beta-2); alpha6=dt*(1-gaJa); alpha7=gaJa*dt;
Kn=K+alpha0*J+alpha1*C; %%%%形成等效刚度阵
%%%%不同时刻等效荷载Pn,位移u,速度v,加速度a的迭代计算
for i=1:1:n
P(:,i+1)=[(M0/tend)*dt*i;0;0;0;0;0;(-M0/tend)*dt*i];
Pn(:,i+1)=P(:,i+1)+J*(alpha0*u(:,i)+alpha2*v(:,i)+alpha3*a(:,i)) ...
+C*(alpha1*u(:,i)+alpha4*v(:,i)+alpha5*a(:,i));
u(:,i+1)=Kn\Pn(:,i+1);
a(:,i+1)=alpha0*(u(:,i+1)-u(:,i))-alpha2*v(:,i)-alpha3*a(:,i);
v(:,i+1)=v(:,i)+alpha6*a(:,i)+alpha7*a(:,i+1);
end
N=400; %总时间段数
for j=n+1:1:N
P2=;
Pn(:,j+1)=P2+J*(alpha0*u(:,j)+alpha2*v(:,j)+alpha3*a(:,j)) ...
+C*(alpha1*u(:,j)+alpha4*v(:,j)+alpha5*a(:,j));
u(:,j+1)=Kn\Pn(:,j+1);
a(:,j+1)=alpha0*(u(:,j+1)-u(:,j))-alpha2*v(:,j)-alpha3*a(:,j);
v(:,j+1)=v(:,j)+alpha6*a(:,j)+alpha7*a(:,j+1);
end
%%%%输出
u1=u'
[ 本帖最后由 ChaChing 于 2009-4-9 13:19 编辑 ] LZ这个不是编程问题了!?
应属专业问题, 待专家并有兴趣者路过
页:
[1]