Graph not showing up when i run code
1 次查看(过去 30 天)
显示 更早的评论
no graphs show up when i run this code
clear all;
clc;
close all;
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODAL ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
clc
k1 = 0.5e9;
k2 = 1e9;
M1 = 65040;
M2 = 62080;
%Assumed modal damping:
z1= 0.001;
z2= 0.1;
%
%Define effective mass and stiffness matrices:
M=[M1,0;0,M2];
K=[k1+k2,-k2;-k2,k2];
%
[V,D]=eig(inv(M)*K);
Nat_freq=sqrt(D); %Natural frequencies in rad/s
Nat_freq_Hz=Nat_freq/(2*pi) %Natural frequencies in Hz
%
%Vnorm - normalised modal matrix
Normal_modes=[V(:,1)/V(1,1) V(:,2)/V(2 ,2)];
Vnorm_columns=[V(:,1)/V(1,1) V(:,2)/V(1,2)];
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GENERATE FORCE (Square Pulses) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=1024;
%fs=400 %Sample rate (samples/s)
fs=1000;
f_Nyquist=fs/2 %Nyquist frequency
df=fs/N %Frequency resolution
dt=1/fs %Time increment of samples
%Tlength=dt*(N-1) %Simulation duration in s
Tlength= 40.96; %Duration of simulation(s)[40.96]
t=0:dt:Tlength; %time array with increment dt
Ndata=length(t) %Data length over the duration of simulation
%
%======= Generation of random excitation ========
%%%CHANGE TO YOUR DATA
F0=25e3 %Amplitude of sinusoidal excitation
Force=0.5*F0*(square(t.*(pi),5)+1); %Zero mean with amplitude F0
Force=Force';
%======= End of random signal generation =======
%
figure(1)
plot(t,Force,'b') %Plot 'Force' time history
grid
xlabel('Time(s)')
ylabel('Force(N)')
title('Excitation force time history')
legend('Random signal')
%
%Show frequency domain features of the excitation
fft_force=fft(Force,N); %FFT of 'Force'
fft_force=fft_force'; %Change from row vector to column vector
freq=0:df:df*(N-1); %Generate a frequency array with increment df
%
figure(2)
plot(freq,db(abs(fft_force)),'r') %Plot mag of 'fft_force' on db scale
grid on
xlabel('Frequency(Hz)')
ylabel('Magnitude of fft force')
title('Magnitude of FFT of the excitation force')
%
%Right-Hand-Side(RHS) force vector of EoM:
%NB: One row for one EoM
F=[zeros(1,Ndata);Force']; %F size is 2xNdata
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUPERPOSITION METHOD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Defines the normal mode matrix from the Modal Analysis results
%============================ Begin [4] =============================
disp('Normal modes from modal analysis are used below')
X1=Vnorm_columns(:,1);
X2=Vnorm_columns(:,2);
disp('Modal matrix X:')
X=[X1 X2] %Form modal matrix for subsequent solution
%Modal mass and stiffness matrices
disp('Modal mass matrix Mg:')
Mg=X'*M*X;
disp('Modal stiffness matrix Kg:')
Kg=X'*K*X;
%
%Modal force vector for Ndata time steps (Ndata columns)
Fg=X'*F; %A 4xNdata matrix with one row for one mode of de-coupled EoMs
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NUMERICAL METHOD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Part 6 solves the de-coupled SDOF EoM for each mode
%============================ Begin [6] =============================
Tspan=[0 Tlength];
q0=[0;0];
%
%NB: The user function file 'Example1_AnyInput.m' should be located in the
%same folder.
%
%Solving the 1st mode
mq=Mg(1,1);
kq=Kg(1,1);
cq = 2*z1*Nat_freq(1,1)*Mg(1,1);
fq=Fg(1,:);
[T1 Q1]=ode45(@(t,y) Example1_AnyInput(t,y,mq,kq,cq,fq,dt),Tspan,q0);
%
%Solving the 2nd mode
mq=Mg(2,2);
kq=Kg(2,2);
cq = 2*z2*Nat_freq(2,2)*Mg(2,2);
fq=Fg(2,:);
[T2 Q2]=ode45(@(t,y) Example1_AnyInput(t,y,mq,kq,cq,fq,dt),Tspan,q0);
%
%Below change the response output from 'ode45' with variable-stepsize into
%fixed stepsize with time array 't' using interpolation function 'PCHIP'
Q1N=interp1(T1,Q1,t,'PCHIP');
Q2N=interp1(T2,Q2,t,'PCHIP');
%plots the modal co-ordinate time history for each mode
figure(31)
plot(t,Q1N(:,1))
grid on
xlabel('Time(s)')
ylabel('Displacement q1')
title('1st modal co-ordinate time history')
%
figure(32)
plot(t,Q2N(:,1))
grid on
xlabel('Time(s)')
ylabel('Displacement q2')
title('2nd modal co-ordinate time history')
% Applies modal superposition transformation: x=X*q
x_rows=X*[Q1N(:,1)';Q2N(:,1)']; % Real Displacement
x_dot =X*[Q1N(:,2)';Q2N(:,2)']; % Velocity
%NB: dimensions in above multiplication X(2x2)*Q(2xNdata)-> x(2xNdata)
x=x_rows'; %Transpose to put each displacement history as a column
figure(41)
plot(t,x(:,1))
grid on
xlabel('Time(s)')
ylabel('Displacement x1(m)')
title('m1 displacement response')
%
figure(42)
plot(t,x(:,2))
grid on
xlabel('Time(s)')
ylabel('Displacement x2(m)')
title('m2 displacement response')
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FREQUENCY ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculate the FFT of modal co-ordinates
fft_q1=fft(Q1N(:,1),N); %FFT with N data length
fft_q2=fft(Q2N(:,1),N);
%
%Below plot magnitude of FFT of modal co-ordinates
%
%Plot on linear scale
figure(5)
plot(freq,(abs(fft_q1)),'r')
grid on
hold on %to allow adding more curves to the same plot window
plot(freq,(abs(fft_q2)),'b')
hold off %to end adding more curves to this plot window
xlabel('Frequency(Hz)')
ylabel('Mag of FFT of Modal Disp')
title('Mag of FFT of Modal response plotted in [0,fs]')
legend('FFT for mode 1','FFT for mode 2')
xlim([0 f_Nyquist]) %plot in frequency range [0,f_Nyquist]
%
%Plot on dB scale
figure(6)
plot(freq,db(abs(fft_q1)),'r')
grid on
hold on
plot(freq,db(abs(fft_q2)),'b')
hold off
xlabel('Frequency(Hz)')
ylabel('Mag of FFT of Modal Disp (dB)')
title('Mag of FFT of Modal response plotted in [0,fs/2]')
legend('FFT for mode 1','FFT for mode 2')
xlim([0 f_Nyquist]) %plot in frequency range [0,f_Nyquist]
% calculates and plots the FFT of physical co-ordinates
%============================== Begin [11] ==============================
%Calculate the FFT of the displacement of each mass
fft_d1=fft(x(:,1),N);
fft_d2=fft(x(:,2),N);
%
%Below plot magnitude of FFT of the displacements of the masses
%
%Plot on linear scale
figure(7)
plot(freq,(abs(fft_d1)),'r')
grid on
hold on
plot(freq,(abs(fft_d2)),'b')
hold off
xlim([0 f_Nyquist]) %plot in frequency range [0,f_Nyquist]
xlabel('Frequency(Hz)')
ylabel('Mag of FFT of mass displacement')
title('Mag of FFT of mass displacement plotted in [0,fs/2]')
legend('FFT of m1 disp','FFT of m2 disp')
%
%Plot on dB scale
figure(8)
plot(freq,db(abs(fft_d1)),'r')
grid on
hold on
plot(freq,db(abs(fft_d2)),'b')
hold off
xlabel('Frequency(Hz)')
ylabel('Mag of FFT of Modal Disp (dB)')
title('Mag of FFT of mass displacement (dB) plotted in [0,fs/2]')
legend('FFT of disp 1','FFT of disp 2')
xlim([0 f_Nyquist]) %plot in frequency range [0,f_Nyquist]
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TRANSMISSIBILITY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_dot = x_dot'; %convert the velocity vector from row vector to a column vector
c1 = 2*z1*sqrt(k1*M1); %damping coefficient
Fground = x(:,1).*k1 + x_dot(:,1).*c1; %Force acting on the ground
plot(t,Fground) % plot the force acting on the ground in the time domain
xlabel('time (sec)')
ylabel('Ground Force')
freq=0:df:((df*N)/2)-1; %Generate a frequency array with increment df
L = length(freq);
% Fast Fourier Transform of Input force and Ouput force
FIN = fft(Force,N);
FOUT = fft(Fground,N);
TR = abs(FOUT)./abs(FIN); %Transmissibility
figure(9)
plot(freq,TR(1:L,1)) %plot the Transmissibility in the frequency domain
xlabel('frequency (Hz)')
ylabel('Force Transmissibility')
grid on
hold on
figure(10)
plot(freq,db(TR(1:L,1))) %plot the Transmissibility in the frequency domain dB scale
xlabel('frequency (Hz)')
ylabel('Force Transmissibility (dB)')
grid on
hold on
回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!

