Problem based on Modal period & Newmark-Beta Method

I've generated an MATLAB code based on these pointed problems, such that:

Q: Consider a 2D 3-story shear building model.

1. Form the mass and stiffness matrices for this building model. No damping in needed for this problem.

2. Using trial and errors, select meaningful story mass and stiffness parameters so that the first natural frequency of the building is 1 sec. Consider that the story mass and stiffness values are similar among the three stories; therefore, you need to set only one stiffness and one mass parameter. Use your engineering judgement and explain your logic. Find and plot the mode shapes of the building model.

3. Implement a Newmark-beta algorithm for time integration of the dynamic equilibrium equation. Use a unidirectional earthquake base excitation as input and solve for the acceleration response at the three DOFs. For the earthquake base excitation, use this input file with a sampling frequency of 100 Hz. Plot the time history of absolute acceleration response at each floor.

4. Resample the input signal at 20, 10, and 2 Hz and redo steps 4-5 with integration time step size of 0.05, 0.1, and 0.5 sec, respectively. Compare the results and explain your observations.

MATLAB CODE

-----------------------------------------------------------------------------------------------------------------

clear all;

clc;

close all;

%% INPUT Parameters:

m = 1; % Initial arbitrary mass of each story

k = 1; % Initial arbitrary stiffness of each story

%% Solution-1: Form Mass and Stiffness Matrices

M = m * eye(3); % Mass matrix

K = [ 2*k, -k, 0; -k, 2*k, -k; 0, -k, k]; % Stiffness matrix

%% Solution-2: Adjust Mass and Stiffness based on target frequency

target_frequency = 1; % Natural Frequency (Hz)

target_omega = 2 * pi * target_frequency; % Angular Frequency (rad/s)

% Adjust stiffness based on the target frequency

k = target_omega^2 * m / 2;

K = [ 2*k, -k, 0; -k, 2*k, -k; 0, -k, k]; % Updated Stiffness matrix

%% Solution-3: Find and Plot Mode Shapes

[eigVectors, eigValues] = eig(K, M); % Solve eigenvalue problem

omega = sqrt(diag(eigValues)); % Natural frequencies (rad/s)

figure;

for i = 1:3

subplot(3, 1, i);

plot([0 eigVectors(:,i)'], 'o-'); % Plot mode shapes

title(['Mode Shape ', num2str(i)]);

xlabel('Story');

ylabel('Displacement');

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

--------------------------------------------------------------------------------------------------------------------

I'm getting error in obtaining:

1. The first modal period as 1 second.
2. The correct Newmark Beta plots.