冲击响应谱计算的matlab程序

《冲击响应谱计算的matlab程序》由会员分享，可在线阅读，更多相关《冲击响应谱计算的matlab程序（7页珍藏版）》请在装配图网上搜索。

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

2、%clear t;clear y;clear yy;clear n;clear fn;clear al;clear a2clear bl;clear b2;clear jnum;clear THM;clear resp;clear x_pos;clear x_neg;%iunit二input( Enter acceleration unit:1= G 2= m./sec 2);%disp()dispC Select file input method );dispCl=external ASCII file );dispC 2=file preloaded into Matlab );file

3、_choice = input ()；%if (file_choice=l)filename, pathname二 uigetf *);filename 二 fullfile(pathname, filename):%fid 二 fopen (filename, r);THM = fscanf(fid/%g %g , 2 inf);THM二THM;elseTHM 二 inp ut ( Enter the matrix name:);end%trouble(THM(:, 1);y=double(THM(:, 2);%tmx=max(t);tmi=min(t);n 二 length(y);%out

4、l 二 sprintf (n %d samples n，, n):disp(outl)%dt二(tmx-tmi)/ (nl);sr= 1 /dt;%outl 二 sprintf C SR 二 %g samples/secdt 二 %g sec n , sr, dt);disp(outl)%fn 二 input ( Enter the starting frequency (Hz) );if fn(l) sr/30.fn(l)=sr/30.;end %idamp=input (J Enter damping format:1= damping ratio 2= Q ，)；%disp()if(id

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

6、1, limit);for i=l:length(y)yy(i)=y(i);end disp( )% SRS engine for j二1:1000%omega二2 *pi*fn(j);omegad=omega*sqrt(1 -(damp 2); cosd二cos(omegad*dt);sind二sin(omegad*dt);domegadt二damp*omega*dt;if (ialgorithm=l)al (j)二2. *exp(-domegadt) *cosd;a2 (j)二-exp (-2 *domegadt);bl (j)二2. domegadt;b2 (j)=omega*dt*ex

7、p(-domegadt);b2 (j)=b2 (j) * ( (omega/omegad) * (1 一2. * (damp 2) *sind一2.*dampcosd );b3(j)=0;%elseE=exp(-damp*omega*dt);K=omegad*dt;C二E*cos (K);S二E*sin(K);Sp=S/K;%al(j)=2*C; a2(j)=-E2; bl (j)=l. -Sp; b2(j)=2.*(Sp-C); b3(j)=E2-SP;end forward= bl (j),b2(j),b3(j);back二1, -al(j), -a2(j);%resp二filter(fo

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

9、a as);writepfname = fullfile(writepname, writefname): writedata 二 fn x_pos (abs(x_neg);fid 二 fopen(writepfname,);fprintf (fid, %g %g %gn，, writedata); fclose (fid);% disp( Enter output filename );% SRS_filename = input ( , s);% fid = f open(SRS_f i1ename,，w)；% for j=l:jnum%fprintf (fid, %7. 2f %10.

10、3f %10. 3fn , fn (j), x_pos (j), abs (x_neg (j);% end% fclose(fid);end% Plot SRS%disp()disp( Plot ting output )% Find limits for plot%srs_max 二 max(x_pos);if max ( abs(x_neg) ) srs_maxsrsmax 二 max( abs(x_neg );endsrs_min 二 min(x_pos);if min ( abs(x_neg) ) = 0. 1fmin二 0. 1;endif fn(l) = 1fmin二1;endif

11、 fn(l) = 10fmin二10;endif fn(l) = 100fmin二100;endaxis (Lfmin, fmax, ymin, ymax);%disp( )dispC Plot pseudo velocity?，)；vchoice二input( l=yes 2=no );if (vchoice二二1)figure(2);% Convert to pseudo velocity%for j=l:jnumif iunit=lx_pos (j)=386. *x_pos (j) / (2. *pi*fn (j);x_neg (j)=386. *x_neg (j) / (2. *pi*

12、fn(j);elseX_pos (j) =x_pos (j)/(2. *pi*fn(j);x_neg(j) =x_neg(j)/(2. *pi*fn(j);endend%srs_max 二 max(x_pos);if max ( abs(x_neg) ) srs_maxsrsmax 二 max( abs(x_neg );endsrs_min 二 min(x_pos);if min ( abs(x_neg) ) srs_minsrsmin 二 min( abs(x_neg );end%plot (fn, x_pos, fn, abs (x_neg)、- )；%if iunit二二1ylabel

13、(? Velocity (in/sec);elseylabel ( Velocity (m/sec);endxlabel( 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 J :，XScale , log ，,，YScale,，log);legend ( positive , negative , 2);%ymax二 10 (round(loglO(srs_max)+0. 8);ymin二 10 (round (loglO(srs_min)-0. 6);%axis (fmin, fmax, ymin, ymax); end