我现在用的就是wilson法求解振动方程，结果不收敛。正不知道怎么办呢？请各位给我解答好吗？谢谢！
我的源程序如下：
%wilson-sita法，采用增量形式
sita=1.4;ts=0.01;bt=sita*ts;
u=repmat(0,n-1,1);
v=repmat(0,n-1,1);
U0=repmat(0,2*(n-1),q);%在wilson法求解中各时间步长下的位移；
U1=repmat(0,2*(n-1),q);%在wilson法求解中各时间步长下的速度；
U2=repmat(0,2*(n-1),q);%在wilson法求解中各时间步长下的加速度；
pk=(6*M/(bt^2)+3*C/bt+K);%有效刚度矩阵pk；
tU0=[u;v];tU1=[u;v];tU2=[u;v];%初始位移、速度和加速度变量；
for i=1:q
TU=pinv(pk)*(F(1:2*(n-1),i+1)-F(1:2*(n-1),i)+M*(6*tU1/bt+3*tU2)+C*(3*tU1+ts*tU2/2));
dU2=6*TU/bt^2-6*tU1/bt-3*tU2;%加速度增量；
U2(1:2*(n-1),i)=tU2+dU2/sita;
dU1=ts*tU2+ts*dU2/2;%速度增量；
U1(1:2*(n-1),i)=tU1+dU1;
dU0=ts*tU1+ts^2*tU2/2+ts^2*dU2/6;%位移增量；
U0(1:2*(n-1),i)=tU0+dU0;
tU0=U0(1:2*(n-1),i);
tU1=U1(1:2*(n-1),i);
tU2=U2(1:2*(n-1),i);
end
maker系列--威尔逊theta法
function maker4(m,c,k,x0,v0,time,delt)
%线性加速度法计算三自由度简谐受迫振动:改编于尚涛等的书
%将beta改为1/4，即为纽马克beta法
%m-为质量，c-阻尼，k-刚度，x0-初位移，v0-除速度；time-仿真时间，delt-时间步长,num-自由度数
% MX''+CX'+KX=F
% X(t+delt)=[M+delt*C/2+delt^2*K/6]^(-1){F(t+delt)-C[X(t)+delt*X(t)'']-K[X(t)+delt*X(t)'+delt^2*X(t)''/3]}
% X(t+delt)'=X(t)'+delt*[X(t)''+X(t+delt)]/2
% X(t+delt)=X(t)+delt*X(t)'+(delt^2)*X(t)''/3+(delt^2)*X(t+delt)''/6
%http://www.dytrol.com [email protected]
%2003.11.27
n=time/delt;
disp=zeros(3,n);%设定存储位移矩阵
m_inv=inv(m+delt*c/2+delt^2*k/6);
[mod,fre]=eig(inv(m)*k);
diag(sqrt(fre));%固有频率
i=1;
beta=1/4;
for t=0:delt:time
f=[2.0*sin(3.5*t) -2.0*cos(2*t) 1.0*sin(3*t)]';%外扰力
if t==0
xdd0=inv(m)*(f-k*x0-c*v0);%初始加速度
else
xdd=m_inv*(f-c*(v0+delt/2*xdd0)-k*(x0+delt*v0+(1/2-beta)*delt^2*xdd0));%计算加速度
x=m_inv*(m*(x0+delt*v0+delt^2/3*xdd0)+c*(delt/2*x0+delt^2/3*v0+delt^3/12*xdd0)+delt^2/6*f);%计算位移
xd=v0+delt*(xdd0+xdd)/2;%计算速度
xdd0=xdd;
v0=xd;
x0=x;
disp(:,i)=x0;
i=i+1;
end
end
t=1:n;
figure('numbertitle','off','name','自由度1的位移','pos',[450 180 400 420]);
plot(t,disp(1,:)),grid,xlabel('时间(s/10)'),title('自由度1的时程曲线');
figure('numbertitle','off','name','自由度2的位移','pos',[450 180 400 420]);
plot(t,disp(2,:)),grid,xlabel('时间(s/10)'),title('自由度2的时程曲线');
figure('numbertitle','off','name','自由度3的位移','pos',[450 180 400 420]);
plot(t,disp(3,:)),grid,xlabel('时间(s/10)'),title('自由度3的时程曲线');
%end
实例
>> m=2*[1 0 0;0 1 0;0 0 1];
c=1.5*[2 -1 0;-1 2 -1;0 -1 2];
k=50*[2 -1 0;-1 2 -2;0 -1 2];
x0=[1 1 1]';
v0=[1 1 1]';
delt=0.1;
time=100;
>> maker4(m,c,k,x0,v0,time,delt)
来源于：动力学与控制技术论坛 http://www.dytrol.com
评论0