Earthquake Response Spectrum Analysis Using MATLAB

I found an old assignment from a long time ago in the trash bin…

This is a frequency-time hybrid algorithm based on the Newmark-β implicit integration method, utilizing the Piecewise Exact Method (PEM) to perform batch calculations of the elastic response spectrum for a single degree of freedom (SDOF) system under seismic excitation. The seismic acceleration time history is treated as a piecewise linear load, and the state transition matrix A and load interpolation matrix B are constructed using the analytical solution of Duhamel’s integral. The displacement, velocity, absolute acceleration, and pseudo-spectral quantities are solved recursively with unconditional stability and second-order accuracy for any damping ratio ζ and structural period T. First, the seismic records are read in units of gravitational acceleration g, converted to international standard units using a conversion factor of 9.8067 m/s², and then a parameter scan is performed for damping ratios ζ ∈ {0, 0.05, 0.10, 0.20} and structural periods T ∈ [0.1, 6] s (with a step size of 1 ms) to calculate and cache the peak responses for each condition. The program outputs the seismic time history curve and the absolute acceleration response spectrum, relative displacement response spectrum, relative velocity response spectrum, pseudo-acceleration response spectrum, and pseudo-velocity response spectrum corresponding to different damping levels.

The MATLAB plots are as follows:

Earthquake Response Spectrum Analysis Using MATLABEarthquake Response Spectrum Analysis Using MATLABEarthquake Response Spectrum Analysis Using MATLABEarthquake Response Spectrum Analysis Using MATLABEarthquake Response Spectrum Analysis Using MATLABEarthquake Response Spectrum Analysis Using MATLAB

The MATLAB code is as follows:

% Clear environment  
clear;  
clc;  
close all;  
  
% Read seismic records  
% Assume the seismic wave data is stored in the 'filename' file  
% The first column is time, the second column is acceleration (in g)  
filename = '01.txt'; % Please replace with the actual file name  
data = readmatrix(filename);  
time = data(:, 1); % Time column  
acceleration_g = data(:, 2); % Acceleration column (in g)  
acceleration = acceleration_g * 9.8067; % Convert acceleration to m/s^2  
  
% Initialize parameters  
dampingRatios = [0, 0.05, 0.1, 0.2]; % Damping ratio array  
structuralPeriods = 0.1:0.001:6; % Structural period range  
timeStep = time(2) - time(1); % Time interval  
N = length(acceleration); % Length of acceleration data  
  
% Initialize matrix to store maximum responses for different damping ratios  
maxAbsAccResponses = zeros(length(dampingRatios), length(structuralPeriods)); % Absolute acceleration  
maxRelDispResponses = zeros(length(dampingRatios), length(structuralPeriods)); % Relative displacement  
maxRelVelResponses = zeros(length(dampingRatios), length(structuralPeriods)); % Relative velocity  
maxPseudoAccResponses = zeros(length(dampingRatios), length(structuralPeriods)); % Pseudo-acceleration  
maxPseudoVelResponses = zeros(length(dampingRatios), length(structuralPeriods)); % Pseudo-velocity  
  
% Calculate seismic responses  
for k = 1:length(dampingRatios)  
    params.dampingRatio = dampingRatios(k); % Current damping ratio  
  
    % Initialize storage vectors  
    response.displacement = zeros(1, N); % Relative displacement  
    response.velocity = zeros(1, N); % Relative velocity  
    response.absoluteAcceleration = zeros(1, N); % Absolute acceleration  
    response.pseudoAcceleration = zeros(1, N); % Pseudo-acceleration  
    response.pseudoVelocity = zeros(1, N); % Pseudo-velocity  
  
    for idx = 1:length(structuralPeriods)  
        T = structuralPeriods(idx); % Current structural period  
  
        % Calculate natural frequency and damped frequency  
naturalFrequency = 2 * pi / T;  
dampedFrequency = naturalFrequency * sqrt(1 - params.dampingRatio^2);  
  
        % Pre-calculate common values  
e_t = exp(-params.dampingRatio * naturalFrequency * timeStep);  
s = sin(dampedFrequency * timeStep);  
c = cos(dampedFrequency * timeStep);  
d_f = (2 * params.dampingRatio^2 - 1) / (naturalFrequency^2 * timeStep);  
d_3t = params.dampingRatio / (naturalFrequency^3 * timeStep);  
  
        % Assemble A matrix  
        A = [  
            e_t * (s * params.dampingRatio / sqrt(1 - params.dampingRatio^2) + c), e_t * s / dampedFrequency;  
            -naturalFrequency * e_t * s / sqrt(1 - params.dampingRatio^2), e_t * (-s * params.dampingRatio / sqrt(1 - params.dampingRatio^2) + c)  
        ];  
  
        % Assemble B matrix  
        B = [  
            e_t * ((d_f + params.dampingRatio / naturalFrequency) * s / dampedFrequency + (2 * d_3t + 1 / naturalFrequency^2) * c) - 2 * d_3t, ...  
            -e_t * (d_f * s / dampedFrequency + 2 * d_3t * c) - 1 / naturalFrequency^2 + 2 * d_3t;  
            e_t * ((d_f + params.dampingRatio / naturalFrequency) * (c - params.dampingRatio / sqrt(1 - params.dampingRatio^2) * s) - ...  
            (2 * d_3t + 1 / naturalFrequency^2) * (dampedFrequency * s + params.dampingRatio * naturalFrequency * c)) + 1 / (naturalFrequency^2 * timeStep), ...  
            e_t * (1 / (naturalFrequency^2 * timeStep) * c + s * params.dampingRatio / (naturalFrequency * dampedFrequency * timeStep)) - 1 / (naturalFrequency^2 * timeStep)  
        ];  
  
        % Calculate responses based on seismic records  
        for i = 1:(N-1)  
            response.displacement(i+1) = A(1,1) * response.displacement(i) + A(1,2) * response.velocity(i) + B(1,1) * acceleration(i) + B(1,2) * acceleration(i+1);  
            response.velocity(i+1) = A(2,1) * response.displacement(i) + A(2,2) * response.velocity(i) + B(2,1) * acceleration(i) + B(2,2) * acceleration(i+1);  
            response.absoluteAcceleration(i+1) = -2 * params.dampingRatio * naturalFrequency * response.velocity(i+1) - naturalFrequency^2 * response.displacement(i+1);  
              
            % Calculate pseudo-acceleration and pseudo-velocity  
            response.pseudoAcceleration(i+1) = response.absoluteAcceleration(i+1) + acceleration(i+1); % Pseudo-acceleration is the sum of absolute acceleration and ground acceleration  
            response.pseudoVelocity(i+1) = response.velocity(i+1) + (acceleration(i+1) - acceleration(i)) * timeStep; % Pseudo-velocity is relative velocity plus the velocity change caused by ground acceleration  
        end  
  
        % Store maximum responses  
        maxAbsAccResponses(k, idx) = max(abs(response.absoluteAcceleration));  
        maxRelDispResponses(k, idx) = max(abs(response.displacement));  
        maxRelVelResponses(k, idx) = max(abs(response.velocity));  
        maxPseudoAccResponses(k, idx) = max(abs(response.pseudoAcceleration));  
        maxPseudoVelResponses(k, idx) = max(abs(response.pseudoVelocity));  
  
        % Reset storage vectors for the next period's calculation  
        response.displacement = zeros(1, N);  
        response.velocity = zeros(1, N);  
        response.absoluteAcceleration = zeros(1, N);  
        response.pseudoAcceleration = zeros(1, N);  
        response.pseudoVelocity = zeros(1, N);  
    end  
end  
  
% Plot charts  
% Plot seismic record chart  
figure;  
plot(time, acceleration);  
title('EARTHQUAKE RECORD');  
xlabel('Time (s)');  
ylabel('Acceleration (m/s^2)');  
grid on;  
  
outputFileName1 = 'EARTHQUAKE RECORD.png';  
exportgraphics(gca, outputFileName1, 'Resolution', 500);  
  
% Plot absolute acceleration response spectrum  
figure;  
hold on;  
colors = {'b', 'r', 'g', 'm'}; % Define color array  
for k = 1:length(dampingRatios)  
    plot(structuralPeriods, maxAbsAccResponses(k, :), ['-', colors{k}]);  
    legendEntries{k} = ['\zeta = ', num2str(dampingRatios(k))]; % Create legend entry  
end  
hold off;  
title('Absolute Acceleration Response Spectrum');  
xlabel('Structural Period (s)');  
ylabel('Absolute Acceleration (m/s^2)');  
legend(legendEntries, 'Location', 'Best');  
grid on;  
  
outputFileName2 = 'Absolute Acceleration Response Spectrum.png';  
exportgraphics(gca, outputFileName2, 'Resolution', 500);  
  
% Plot relative displacement response spectrum  
figure;  
hold on;  
for k = 1:length(dampingRatios)  
    plot(structuralPeriods, maxRelDispResponses(k, :), ['-', colors{k}]);  
end  
hold off;  
title('Relative Displacement Response Spectrum');  
xlabel('Structural Period (s)');  
ylabel('Relative Displacement (m)');  
legend(legendEntries, 'Location', 'Best');  
grid on;  
  
outputFileName3 = 'Relative Displacement Response Spectrum.png';  
exportgraphics(gca, outputFileName3, 'Resolution', 500);  
  
% Plot relative velocity response spectrum  
figure;  
hold on;  
for k = 1:length(dampingRatios)  
    plot(structuralPeriods, maxRelVelResponses(k, :), ['-', colors{k}]);  
end  
hold off;  
title('Relative Velocity Response Spectrum');  
xlabel('Structural Period (s)');  
ylabel('Relative Velocity (m/s)');  
legend(legendEntries, 'Location', 'Best');  
grid on;  
  
outputFileName4 = 'Relative Velocity Response Spectrum.png';  
exportgraphics(gca, outputFileName4, 'Resolution', 500);  
  
% Plot pseudo-acceleration response spectrum  
figure;  
hold on;  
for k = 1:length(dampingRatios)  
    plot(structuralPeriods, maxPseudoAccResponses(k, :), ['-', colors{k}]);  
end  
hold off;  
title('Pseudo-Acceleration Response Spectrum');  
xlabel('Structural Period (s)');  
ylabel('Pseudo-Acceleration (m/s^2)');  
legend(legendEntries, 'Location', 'Best');  
grid on;  
  
outputFileName5 = 'Pseudo-Acceleration Response Spectrum.png';  
exportgraphics(gca, outputFileName5, 'Resolution', 500);  
  
% Plot pseudo-velocity response spectrum  
figure;  
hold on;  
for k = 1:length(dampingRatios)  
    plot(structuralPeriods, maxPseudoVelResponses(k, :), ['-', colors{k}]);  
end  
hold off;  
title('Pseudo-Velocity Response Spectrum');  
xlabel('Structural Period (s)');  
ylabel('Pseudo-Velocity (m/s)');  
legend(legendEntries, 'Location', 'Best');  
grid on;  
  
outputFileName6 = 'Pseudo-Velocity Response Spectrum.png';  
exportgraphics(gca, outputFileName6, 'Resolution', 500);  
  
% Set output file names and paths  
% outputFileName1 = 'EARTHQUAKE RECORD.png';  
% outputFileName2 = 'Absolute Acceleration Response Spectrum.png';     
% outputFileName3 = 'Relative Displacement Response Spectrum.png';  
% outputFileName4 = 'Relative Velocity Response Spectrum.png';  
% outputFileName5 = 'Pseudo-Acceleration Response Spectrum.png';  
% outputFileName6 = 'Pseudo-Velocity Response Spectrum.png';  
  
% Use exportgraphics function to export images, specifying resolution of 300 DPI  
% exportgraphics(gca, outputFileName1, 'Resolution', 500);  
% exportgraphics(gca, outputFileName2, 'Resolution', 500);  
% exportgraphics(gca, outputFileName3, 'Resolution', 500);  
% exportgraphics(gca, outputFileName4, 'Resolution', 500);  
% exportgraphics(gca, outputFileName5, 'Resolution', 500);  
% exportgraphics(gca, outputFileName6, 'Resolution', 500);

Related posts

  1. MATLAB Signal Processing | Fast DCT Algorithm Using FFT (Complete Code Included)
  2. Analysis of Spring-Mass-Damper System Behavior Using MATLAB
  3. Reliability Verification of PCBA-SMT-AOI Visual Inspection Equipment
  4. Constant Amplitude Multi-Beamforming Algorithm for Millimeter Wave Communication Based on MATLAB [Includes MATLAB Source Code]
  5. Deploying Wireless Sensor Networks Using MATLAB with Sample Code
  6. New HTTP Smuggling Attack Technique Allows Attackers to Inject Malicious Requests

Leave a Comment