冲击响应谱计算的matlab程序.docx
- 文档编号:10437611
- 上传时间:2023-05-25
- 格式:DOCX
- 页数:10
- 大小:15.78KB
冲击响应谱计算的matlab程序.docx
《冲击响应谱计算的matlab程序.docx》由会员分享,可在线阅读,更多相关《冲击响应谱计算的matlab程序.docx(10页珍藏版)》请在冰点文库上搜索。
冲击响应谱计算的matlab程序
disp('')
disp('verJuly3,2006')
disp('byTomIrvineEmail')
disp('')
disp('Thisprogramcalculatestheshockresponsespectrum')
disp('ofanaccelerationtimehistory,whichispre-loadedintoMatlab.')
disp('Thetimehistorymusthavetwocolumns:
time(sec)&acceleration')
disp('')
%
cleart;
cleary;
clearyy;
clearn;
clearfn;
cleara1;
cleara2
clearb1;
clearb2;
clearjnum;
clearTHM;
clearresp;
clearx_pos;
clearx_neg;
%
iunit=input('Enteraccelerationunit:
1=G2=m/sec^2');
%
disp('')
disp('Selectfileinputmethod');
disp('1=externalASCIIfile');
disp('2=filepreloadedintoMatlab');
file_choice=input('');
%
if(file_choice==1)
[filename,pathname]=uigetfile('*.*');
filename=fullfile(pathname,filename);
%
fid=fopen(filename,'r');
THM=fscanf(fid,'%g%g',[2inf]);
THM=THM';
else
THM=input('Enterthematrixname:
');
end
%
t=double(THM(:
1));
y=double(THM(:
2));
%
tmx=max(t);
tmi=min(t);
n=length(y);
%
out1=sprintf('\n%dsamples\n',n);
disp(out1)
%
dt=(tmx-tmi)/(n-1);
sr=1./dt;
%
out1=sprintf('SR=%gsamples/secdt=%gsec\n',sr,dt);
disp(out1)
%
fn
(1)=input('Enterthestartingfrequency(Hz)');
iffn
(1)>sr/30.
fn
(1)=sr/30.;
end
%
idamp=input('Enterdampingformat:
1=dampingratio2=Q');
%
disp('')
if(idamp==1)
damp=input('Enterdampingratio(typically');
else
Q=input('Entertheamplificationfactor(typicallyQ=10)');
damp=1./(2.*Q);
end
%
disp('')
disp('Selectalgorithm:
')
disp('1=Kelly-Richman2=Smallwood');
ialgorithm=input('');
%
tmax=(tmx-tmi)+1./fn
(1);
limit=round(tmax/dt);
n=limit;
yy=zeros(1,limit);
fori=1:
length(y)
yy(i)=y(i);
end
%
disp('')
disp('Calculatingresponse.....')
%
%SRSengine
%
forj=1:
1000
%
omega=2.*pi*fn(j);
omegad=omega*sqrt(damp^2));
cosd=cos(omegad*dt);
sind=sin(omegad*dt);
domegadt=damp*omega*dt;
%
if(ialgorithm==1)
a1(j)=2.*exp(-domegadt)*cosd;
a2(j)=-exp(-2.*domegadt);
b1(j)=2.*domegadt;
b2(j)=omega*dt*exp(-domegadt);
b2(j)=b2(j)*((omega/omegad)*.*(damp^2))*sind-2.*damp*cosd);
b3(j)=0;
%
else
E=exp(-damp*omega*dt);
K=omegad*dt;
C=E*cos(K);
S=E*sin(K);
Sp=S/K;
%
a1(j)=2*C;
a2(j)=-E^2;
b1(j)=;
b2(j)=2.*(Sp-C);
b3(j)=E^2-Sp;
end
forward=[b1(j),b2(j),b3(j)];
back=[1,-a1(j),-a2(j)];
%
resp=filter(forward,back,yy);
%
x_pos(j)=max(resp);
x_neg(j)=min(resp);
%
jnum=j;
iffn(j)>sr/8.
break
end
fn(j+1)=fn
(1)*(2.^(j*(1./12.)));
end
%
%Outputoptions
%
disp('')
disp('Selectoutputoption');
choice=input('1=plotonly2=plot&outputtextfile');
disp('')
%
ifchoice==2
%%
[writefname,writepname]=uiputfile('*','SaveSRSdataas');
writepfname=fullfile(writepname,writefname);
writedata=[fn'x_pos'(abs(x_neg))'];
fid=fopen(writepfname,'w');
fprintf(fid,'%g%g%g\n',writedata');
fclose(fid);
%%
%disp('Enteroutputfilename');
%SRS_filename=input('','s');
%
%fid=fopen(SRS_filename,'w');
%forj=1:
jnum
%fprintf(fid,'%%%\n',fn(j),x_pos(j),abs(x_neg(j)));
%end
%fclose(fid);
end
%
%PlotSRS
%
disp('')
disp('Plottingoutput.....')
%
%Findlimitsforplot
%
srs_max=max(x_pos);
ifmax(abs(x_neg))>srs_max
srs_max=max(abs(x_neg));
end
srs_min=min(x_pos);
ifmin(abs(x_neg)) srs_min=min(abs(x_neg)); end % figure (1); plot(fn,x_pos,fn,abs(x_neg),'-.'); % ifiunit==1 ylabel('PeakAccel(G)'); else ylabel('PeakAccel(m/sec^2)'); end xlabel('NaturalFrequency(Hz)'); Q=1./(2.*damp); out5=sprintf('AccelerationShockResponseSpectrumQ=%g',Q); title(out5); grid; set(gca,'MinorGridLineStyle','none','GridLineStyle',': ','XScale','log','YScale','log'); legend('positive','negative',2); % ymax=10^(round(log10(srs_max)+); ymin=10^(round(log10(srs_min)); % fmax=max(fn); fmin=fmax/10.; % fmax=10^(round(log10(fmax)+); % iffn (1)>= fmin=; end iffn (1)>=1 fmin=1; end iffn (1)>=10 fmin=10; end iffn (1)>=100 fmin=100; end axis([fmin,fmax,ymin,ymax]); % disp('') disp('Plotpseudovelocity'); vchoice=input('1=yes2=no'); if(vchoice==1) figure (2); % %Converttopseudovelocity % forj=1: jnum ifiunit==1 x_pos(j)=386.*x_pos(j)/(2.*pi*fn(j)); x_neg(j)=386.*x_neg(j)/(2.*pi*fn(j)); else x_pos(j)=x_pos(j)/(2.*pi*fn(j)); x_neg(j)=x_neg(j)/(2.*pi*fn(j)); end end % srs_max=max(x_pos); ifmax(abs(x_neg))>srs_max srs_max=max(abs(x_neg)); end srs_min=min(x_pos); ifmin(abs(x_neg)) srs_min=min(abs(x_neg)); end % plot(fn,x_pos,fn,abs(x_neg),'-.'); % ifiunit==1 ylabel('Velocity(in/sec)'); else ylabel('Velocity(m/sec)'); end xlabel('NaturalFrequency(Hz)'); Q=1./(2.*damp); out5=sprintf('PseudoVelocityShockResponseSpectrumQ=%g',Q); title(out5); grid; set(gca,'MinorGridLineStyle','none','GridLineStyle',': ','XScale','log','YScale','log'); legend('positive','negative',2); % ymax=10^(round(log10(srs_max)+); ymin=10^(round(log10(srs_min)); % axis([fmin,fmax,ymin,ymax]); end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 冲击 响应 计算 matlab 程序