Prepare an MATLAB code based on Newmark Beta Algorithm problem

19 次查看（过去 30 天）

显示更早的评论

YQ 2024-6-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2131656-prepare-an-matlab-code-based-on-newmark-beta-algorithm-problem

评论： Umar 2024-6-25

The problem is attached herewith along with needed text file.

4 个评论
显示 2更早的评论隐藏 2更早的评论

YQ 2024-6-25

I've written this code.

%%%%%%%%%%%%%%%%%%%%%%%%% Assignment Problem %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

clc;

close all;

%% INPUT Parameters:

m = 1; % mass of each story (arbitrary value, to be adjusted)

k = 1; % stiffness of each story (arbitrary value, to be adjusted)

%% Solution-1:

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

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

%% Solution-2:

target_frequency = 1; % Natural Frequency (Hz)

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

k = target_omega^2 * m / 2; % Mass and stiffness Adjustment Based on the target frequency

K = [ 2*k, -k, 0;-k, 2*k, -k;0, -k, k]; % Stiffness matrix Updation with the new k value

%% Solution-3:

[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:

% Parameters for Newmark-beta

beta = 1/4;

gamma = 1/2;

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);

v = zeros(3, 1);

a = M \ (data(1) * ones(3, 1)); % initial acceleration

% Time integration

u_history = zeros(3, length(time));

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;

end

%% Solution-5: Plot time history of absolute acceleration response --

figure;

for i = 1:3

subplot(3, 1, i);

plot(time, u_history(i, :));

title(['Absolute Acceleration Response at DOF ', num2str(i)]);

xlabel('Time (s)');

ylabel('Acceleration (m/s^2)');

end

%% Solution-6:

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);

v = zeros(3, 1);

a = M \ (new_data(1) * ones(3, 1)); % initial acceleration

u_history = zeros(3, length(new_time));

for i = 1:length(new_time)-1

u_pred = u + new_dt * v + new_dt^2 * (0.5-beta) * a;

v_pred = v + new_dt * (1-gamma) * a;

a_new = M \ (new_data(i+1) * ones(3, 1) - K * u_pred);

u = u_pred + beta * new_dt^2 * a_new;

v = v_pred + gamma * new_dt * a_new;

u_history(:, i+1) = u;

end

% Plot time history of absolute acceleration response

figure;

for i = 1:3

subplot(3, 1, i);

plot(new_time, u_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

YQ 2024-6-25

But the time history responses are incorrect.

Umar 2024-6-25

Hi YQ, The line a = M \ (data(1) * ones(3, 1)); calculates the initial acceleration using the earthquake data at the first time step. However, this calculation should be based on the initial displacement and velocity of the system, not the earthquake data.

By initializing the initial acceleration as zero and calculating it based on the initial displacement and velocity, you can ensure that the time history responses are calculated correctly.

Hope this will help resolve your issue.

请先登录，再进行评论。

请先登录，再回答此问题。