Adjust the given code based on Newmark Beta Method based on input base acceleration file.

clear all;
clc;
close all;
%% Solution 1 & 2:
% Define parameters
m = 1; % Mass of each story (kg)
k = 1; % Stiffness of each story (N/m)
% Mass matrix
M = m * eye(3);
% Stiffness matrix
K = [2*k, -k, 0; -k, 2*k, -k; 0, -k, k];
% Solve the eigenvalue problem
[eigenvectors, eigenvalues] = eig(K, M);
% Extract eigenvalues
omega_squared = diag(eigenvalues);
% Convert to natural frequencies (rad/s)
natural_frequencies = sqrt(omega_squared);
% Convert to natural frequencies (Hz)
natural_frequencies_Hz = natural_frequencies / (2 * pi);
% Display results
disp('Natural Frequencies (Hz):');
disp(natural_frequencies_Hz);
% Adjust parameters to get the first natural frequency to 1 Hz
% This requires iterating and adjusting m and k until the first natural frequency is 1 Hz
desired_frequency = 1; % Hz
desired_omega = desired_frequency * 2 * pi; % rad/s
desired_omega_squared = desired_omega^2;
tolerance = 1e-3; % Tolerance for convergence
max_iterations = 1000; % Maximum number of iterations
iteration = 0;
while iteration < max_iterations
    iteration = iteration + 1;
    % Update stiffness matrix with current stiffness
    K = [2*k, -k, 0; -k, 2*k, -k; 0, -k, k];
    % Solve the eigenvalue problem again
    [eigenvectors, eigenvalues] = eig(K, M);
    % Extract the first eigenvalue
    first_omega_squared = eigenvalues(1, 1);
    % Adjust mass and stiffness
    k = k * (desired_omega_squared / first_omega_squared);
    M = m * eye(3);
    % Check convergence
    if abs(first_omega_squared - desired_omega_squared) < tolerance
        break;
    end
end
% Display final parameters
disp('Final Parameters:');
disp(['Mass per story (kg): ', num2str(m)]);
disp(['Stiffness per story (N/m): ', num2str(k)]);
% Display final natural frequencies
natural_frequencies = sqrt(diag(eigenvalues)) / (2 * pi);
disp('Final Natural Frequencies (Hz):');
disp(natural_frequencies);
% Normalize mode shapes for plotting
mode_shapes = eigenvectors ./ max(abs(eigenvectors), [], 1);
% Plot mode shapes
figure;
title('Mode Shapes of the 3-Story Shear Building');
%% Solution-3: Mode Shapes Plotting:
stories = [0; 1; 2; 3]; % Include ground level
for i = 1:3
    subplot(1, 3, i); % Change to a 1x3 grid
    plot([0; mode_shapes(:,i)], stories, '-o');
    title(['Mode Shape ', num2str(i)]);
    xlabel('Amplitude');
    ylabel('Story');
    grid on;
end
%% Solution-4: Implement Newmark-beta algorithm for dynamic analysis
% Parameters for Newmark-beta
beta = 1/4;      
gamma = 1/2;    
% Load earthquake data
data = load('CNP100.txt'); % Load earthquake data
dt = 0.01;                 % Sampling interval for 100 Hz
time = (0:length(data)-1)*dt;
% Initial conditions
u = zeros(3, 1); % Initial displacement
v = zeros(3, 1); % Initial velocity
a = M \((data(1) * ones(3, 1)) - K*u);    % Initial acceleration
% Time integration using Newmark-beta method
u_history = zeros(3, length(time));
a_history = zeros(3, length(time));
a_history(:, 1) = a;
for i = 1:length(time)-1
    % Predictors
    u_pred = u + (dt * v) + ((dt^2) * (0.5-beta) * a);
    v_pred = v + (dt * (1-gamma) * a);
    
    % Solve for acceleration
    a_new = M \ ((data(i+1) * ones(3, 1)) - (K * u_pred));
    
    % Correctors
    u = u_pred + (beta * (dt^2) * a_new);
    v = v_pred + (gamma * dt * a_new);
    
    % Store response
    u_history(:, i+1) = u;
    a_history(:, i+1) = a_new;
    a = a_new;
end
%% Solution-5: Plot time history of absolute acceleration response
figure;
for i = 1:3
    subplot(3, 1, i);
    plot(time, a_history(i, :));
    title(['Absolute Acceleration Response Sampled at 100 Hz at DOF ', num2str(i)]);
    xlabel('Time (s)');
    ylabel('Acceleration (m/s^2)');
end
%% Solution-6: Resample input signal and repeat the analysis
frequencies = [20, 10, 2];
time_steps = [0.05, 0.1, 0.5];
for j = 1:length(frequencies)
    new_fs = frequencies(j);
    new_dt = 1/new_fs;
    new_data = resample(data, new_fs, 100);
    new_time = (0:length(new_data)-1)*new_dt;
    
    % Repeat the Newmark-beta integration with new data
    u = zeros(3, 1); % Initial displacement
    v = zeros(3, 1); % Initial velocity
    a = M \ ((new_data(1) * ones(3, 1)) - K*u); % Initial acceleration
    
    u_history = zeros(3, length(new_time));
    a_history = zeros(3, length(new_time));
    a_history(:, 1) = a;
    
    for i = 1:length(new_time)-1
        % Predictors
        u_pred = u + (new_dt * v) + ((new_dt^2) * (0.5-beta) * a);
        v_pred = v + (new_dt * (1-gamma) * a);
        
        % Solve for acceleration
        a_new = M \ ((new_data(i+1) * ones(3, 1)) - K * u_pred);
        
        % Correctors
        u = u_pred + (beta * (new_dt^2) * a_new);
        v = v_pred + (gamma * new_dt * a_new);
        
        % Store response
        u_history(:, i+1) = u;
        a_history(:, i+1) = a_new;
        a = a_new;
    end
    
    % Plot time history of absolute acceleration response
    figure;
    for i = 1:3
        subplot(3, 1, i);
        plot(new_time, a_history(i, :));
        title(['Absolute Acceleration Response at DOF ', num2str(i), ' (Resampled at ', num2str(new_fs), ' Hz)']);
        xlabel('Time (s)');
        ylabel('Acceleration (m/s^2)');
    end
end

Answer 1

[akshatsood]

Hi @Parvesh Deepan,

I understand that you want to Adjust the given code based on Newmark Beta Method based on input base acceleration file. As per my understandning, you need to perform the following changes to bring this into place.

Initial Conditions - Initialize displacement and velocity vectors:

u = zeros(3, 1); % Initial displacement
v = zeros(3, 1); % Initial velocity

Newmark-Beta Parameters - Ensure "beta" and "gamma" are set for the average acceleration method:

beta = 1/4;      
gamma = 1/2;

Time Integration Loop - Implement the Newmark-Beta method:

for i = 1:length(time)-1
    % Predictors
    u_pred = u + (dt * v) + ((dt^2) * (0.5-beta) * a);
    v_pred = v + (dt * (1-gamma) * a);

    % Solve for acceleration
    a_new = M \ ((data(i+1) * ones(3, 1)) - (K * u_pred));

    % Correctors
    u = u_pred + (beta * (dt^2) * a_new);
    v = v_pred + (gamma * dt * a_new);

    % Store response
    u_history(:, i+1) = u;
    a_history(:, i+1) = a_new;
    a = a_new;
end

Plotting Results - Ensure plots are correctly labeled and display the correct data:

figure;
for i = 1:3
    subplot(3, 1, i);
    plot(time, a_history(i, :));
    title(['Absolute Acceleration Response at DOF ', num2str(i)]);
    xlabel('Time (s)');
    ylabel('Acceleration (m/s^2)');
end

Resampling - Use the "resample" function correctly for new sampling frequencies:

new_data = resample(data, new_fs, 100);

These snippets and explanations should help guide you in making the necessary adjustments to your code.

I hope this helps.