maker系列--威尔逊theta法
maker系列--威尔逊theta法 <BR>function maker4(m,c,k,x0,v0,time,delt) <BR>%线性加速度法计算三自由度简谐受迫振动:改编于尚涛等的书 <BR>%将beta改为1/4，即为纽马克beta法 <BR>%m-为质量，c-阻尼，k-刚度，x0-初位移，v0-除速度；time-仿真时间，delt-时间步长,num-自由度数 <BR>% MX''+CX'+KX=F <BR>% X(t+delt)=^(-1){F(t+delt)-C-K} <BR>% X(t+delt)'=X(t)'+delt*/2 <BR>% X(t+delt)=X(t)+delt*X(t)'+(delt^2)*X(t)''/3+(delt^2)*X(t+delt)''/6 <BR>%http://www.dytrol.com chinamaker@dytrol.com <BR>%2003.11.27 <BR><BR>n=time/delt; <BR>disp=zeros(3,n);%设定存储位移矩阵 <BR>m_inv=inv(m+delt*c/2+delt^2*k/6); <BR>=eig(inv(m)*k); <BR>diag(sqrt(fre));%固有频率 <BR>i=1; <BR>beta=1/4; <BR>for t=0:delt:time <BR> f=';%外扰力 <BR> if t==0 <BR> xdd0=inv(m)*(f-k*x0-c*v0);%初始加速度 <BR> else <BR> xdd=m_inv*(f-c*(v0+delt/2*xdd0)-k*(x0+delt*v0+(1/2-beta)*delt^2*xdd0));%计算加速度 <BR> 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);%计算位移 <BR> xd=v0+delt*(xdd0+xdd)/2;%计算速度 <BR> xdd0=xdd; <BR> v0=xd; <BR> x0=x; <BR> disp(:,i)=x0; <BR> i=i+1; <BR> end <BR>end <BR>t=1:n; <BR>figure('numbertitle','off','name','自由度1的位移','pos',); <BR>plot(t,disp(1,:)),grid,xlabel('时间(s/10)'),title('自由度1的时程曲线'); <BR>figure('numbertitle','off','name','自由度2的位移','pos',); <BR>plot(t,disp(2,:)),grid,xlabel('时间(s/10)'),title('自由度2的时程曲线'); <BR>figure('numbertitle','off','name','自由度3的位移','pos',); <BR>plot(t,disp(3,:)),grid,xlabel('时间(s/10)'),title('自由度3的时程曲线'); <BR>%end <BR><BR>实例 <BR>>> m=2*; <BR>c=1.5*; <BR>k=50*; <BR>x0='; <BR>v0='; <BR>delt=0.1; <BR>time=100; <BR>>> maker4(m,c,k,x0,v0,time,delt) <BR> 不是很明白
页:
[1]