冲击响应谱计算的matlab程序

disp(' ')

disp(' srs.m ver 2.0 July 3, 2006')

disp(' by Tom Irvine Email: tomirvine@aol.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(' ') %

clear t; clear y; clear yy; clear n; clear fn; clear a1; clear a2 clear b1; 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(' ')

disp(' Select file input method '); disp(' 1=external ASCII file ');

disp(' 2=file preloaded into Matlab '); file_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'; else

THM = input(' Enter the matrix name: '); end %

t=double(THM(:,1));

y=double(THM(:,2)); %

tmx=max(t); tmi=min(t); n = length(y); %

out1 = 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=input(' Enter damping ratio (typically 0.05) '); else

Q=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=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.-(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)*(1.-2.*(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)=1.-Sp; 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;

if fn(j) > sr/8.

break end

fn(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 as'); 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(' Enter output filename '); % SRS_filename = input(' ','s'); %

% fid = fopen(SRS_filename,'w'); % for j=1:jnum

% fprintf(fid,'%7.2f .3f .3f \\n',fn(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_max srs_max = max( abs(x_neg )); end

srs_min = min(x_pos);

if min( abs(x_neg) ) < srs_min srs_min = min( abs(x_neg )); end %

figure(1);

plot(fn,x_pos,fn,abs(x_neg),'-.'); %

if iunit==1

ylabel('Peak Accel (G)'); else

ylabel('Peak Accel (m/sec^2)'); end

xlabel('Natural Frequency (Hz)'); Q=1./(2.*damp);

out5 = sprintf(' Acceleration 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)); %

fmax=max(fn); fmin=fmax/10.; %

fmax= 10^(round(log10(fmax)+0.5)); %

if fn(1) >= 0.1 fmin=0.1; end

if fn(1) >= 1 fmin=1; end

if fn(1) >= 10 fmin=10; end

if fn(1) >= 100 fmin=100; end

axis([fmin,fmax,ymin,ymax]); %

disp(' ')

本文来源：https://www.bwwdw.com/article/f806.html