matlab计算单位冲击响应,冲击响应谱计算的matlab程序

《冲击响应谱计算的matlab程序》由会员分享，可在线阅读，更多相关《冲击响应谱计算的matlab程序(7页珍藏版)》请在人人文库网上搜索。

1、disp( )disp( srs.m ver 2.0 July 3, 2006)disp( by Tom Irvine Email: tomirvineaol.com)disp( )disp( This program calculates the shock response spectrum)disp( of an acceleration time history, which is pre-loaded into Matlab.)disp( The time history must have two columns: time(sec) & acceleration)disp( )%。

2、clear t;clear y;clear yy;clear n;clear fn;clear a1;clear a2clear b1;clear b2;clear jnum;clear THM;clear resp;clear x_pos;clear x_neg;%iunit=input( Enter acceleration unit: 1= G 2= m/sec2 );%disp( )disp( Select file input method );disp( 1=external ASCII file );disp( 2=file preloaded into Matlab );fil。

3、e_choice = input();%if(file_choice=1)filename, pathname = uigetfile(*.*);filename = fullfile(pathname, filename);% fid = fopen(filename,r);THM = fscanf(fid,%g %g,2 inf);THM=THM;elseTHM = input( Enter the matrix name: );end%t=double(THM(:,1);y=double(THM(:,2);%tmx=max(t);tmi=min(t);n = length(y);%out。

4、1 = sprintf(n %d samples n,n);disp(out1)%dt=(tmx-tmi)/(n-1);sr=1./dt;%out1 = sprintf( SR = %g samples/sec dt = %g sec n,sr,dt);disp(out1)%fn(1)=input( Enter the starting frequency (Hz) );if fn(1)sr/30.fn(1)=sr/30.;end%idamp=input( Enter damping format: 1= damping ratio 2= Q );%disp( )if(idamp=1)damp。

5、=input( Enter damping ratio (typically 0.05) );elseQ=input( Enter the amplification factor (typically Q=10) );damp=1./(2.*Q);end%disp( )disp( Select algorithm: )disp( 1=Kelly-Richman 2=Smallwood );ialgorithm=input( );%tmax=(tmx-tmi) + 1./fn(1);limit = round( tmax/dt );n=limit;yy=zeros(1,limit);for i。

6、=1:length(y)yy(i)=y(i);end %disp( )disp( Calculating response. )% SRS engine%for j=1:1000%omega=2.*pi*fn(j);omegad=omega*sqrt(1.-(damp2);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。

7、*dt*exp(-domegadt);b2(j)=b2(j)*( (omega/omegad)*(1.-2.*(damp2)*sind -2.*damp*cosd );b3(j)=0;%elseE=exp(-damp*omega*dt);K=omegad*dt;C=E*cos(K);S=E*sin(K);Sp=S/K;%a1(j)=2*C;a2(j)=-E2;b1(j)=1.-Sp;b2(j)=2.*(Sp-C);b3(j)=E2-Sp;endforward= b1(j), b2(j), b3(j) ; back = 1, -a1(j), -a2(j) ; % resp=filter(forw。

8、ard,back,yy);%x_pos(j)= max(resp);x_neg(j)= min(resp);% jnum=j; if fn(j) sr/8.breakendfn(j+1)=fn(1)*(2. (j*(1./12.); end% Output options%disp( )disp( Select output option );choice=input( 1=plot only 2=plot & output text file );disp( )%if choice = 2 %writefname, writepname = uiputfile(*,Save SRS data。

9、 as);writepfname = fullfile(writepname, writefname);writedata = fn x_pos (abs(x_neg) ;fid = fopen(writepfname,w);fprintf(fid, %g %g %gn,writedata);fclose(fid);% disp( Enter output filename );% SRS_filename = input( ,s);% fid = fopen(SRS_filename,w);% for j=1:jnum% fprintf(fid,%7.2f %10.3f %10.3f n,f。

10、n(j),x_pos(j),abs(x_neg(j);% end% fclose(fid);end% Plot SRS%disp( )disp( Plotting output. )% Find limits for plot%srs_max = max(x_pos);if max( abs(x_neg) ) srs_maxsrs_max = max( abs(x_neg );endsrs_min = min(x_pos);if min( abs(x_neg) ) = 0.1fmin=0.1;endif fn(1) = 1fmin=1;endif fn(1) = 10fmin=10;endif。

11、 fn(1) = 100fmin=100;endaxis(fmin,fmax,ymin,ymax);%disp( )disp( Plot pseudo velocity? );vchoice=input( 1=yes 2=no );if(vchoice=1)figure(2);% Convert to pseudo velocity%for j=1:jnum if iunit=1 x_pos(j)=386.*x_pos(j)/(2.*pi*fn(j);x_neg(j)=386.*x_neg(j)/(2.*pi*fn(j); elsex_pos(j)=x_pos(j)/(2.*pi*fn(j);。

12、x_neg(j)=x_neg(j)/(2.*pi*fn(j); endend %srs_max = max(x_pos);if max( abs(x_neg) ) srs_maxsrs_max = max( abs(x_neg );endsrs_min = min(x_pos);if min( abs(x_neg) ) srs_minsrs_min = min( abs(x_neg );end %plot(fn,x_pos,fn,abs(x_neg),-.);%if iunit=1ylabel(Velocity (in/sec);elseylabel(Velocity (m/sec); end。

13、xlabel(Natural Frequency (Hz);Q=1./(2.*damp);out5 = sprintf( Pseudo Velocity Shock Response Spectrum Q=%g ,Q);title(out5);grid;set(gca,MinorGridLineStyle,none,GridLineStyle,:,XScale,log,YScale,log);legend (positive,negative,2);%ymax= 10(round(log10(srs_max)+0.8);ymin= 10(round(log10(srs_min)-0.6);%axis(fmin,fmax,ymin,ymax);end。