Newmark-Beta in dynamic analysis for solving nonlinear stiffness matrix problem

Hello everyone i am solving multi degree freedom of system by using Newmarks beta method in which stiffness coefficient are nonlineare . When i am ruing the code, the result is look like scary . geting displacement value in kilometerers.I am attaching code . Please help me out.

Thanks

__________________________________________________________________________________________________________________

clear all

clc

% Define parameters

n = 3; % Number of degrees of freedom

m = eye(n); % Mass matrix (identity matrix for simplicity)

c = 0.1 * eye(n); % Damping matrix (proportional damping)

k = diag([11, 8, 9]); % Nonlinear stiffness matrix

% Time parameters

dt = 0.01; % Time step

t_final = 10; % Final time

t = 0:dt:t_final;

% Sine loading parameters

omega = 2*pi*0.5; % Frequency of the sine loading

F_sine = @(t) 200*10^3*sin(omega*t); % Sine loading function

% Initial conditions

x0 = zeros(n, 1); % Initial displacement

v0 = zeros(n, 1); % Initial velocity

% Newmark-Beta parameters

beta = 0.25; % Newmark-Beta parameter

gamma = 0.5; % Newmark-Beta parameter

% Initialize arrays

x = zeros(n, length(t)); % Displacement array

v = zeros(n, length(t)); % Velocity array

a = zeros(n, length(t)); % Acceleration array

% Initial conditions

x(:,1) = x0;

v(:,1) = v0;

a(:,1) = m \ (-c*v(:,1) - NonlinearStiffness(x(:,1), k));

% Newmark-Beta iteration

for i = 2:length(t)

% Predictor step

x_pred = x(:,i-1) + dt*v(:,i-1) + (dt^2/2) * ((1-2*beta)*a(:,i-1));

v_pred = v(:,i-1) + dt*((1-gamma)*a(:,i-1));

% Corrector step

a_pred = m \ (-c*v_pred - NonlinearStiffness(x_pred, k) + F_sine(t(i)));

x(:,i) = x(:,i-1) + dt*v_pred + (dt^2/2) * (beta*a_pred + (1-beta)*a(:,i-1));

v(:,i) = v_pred + dt*((1-gamma)*a_pred + gamma*a(:,i-1));

% Update acceleration for next iteration

a(:,i) = a_pred;

end

% Plot results

plot(t, x(1,:));

xlabel('Time');

ylabel('Displacement DOF 1 ')

% Nonlinear stiffness function

function F = NonlinearStiffness(x, k)

F = k * x.^3; % Element-wise cubic nonlinear stiffness

end

________________________________________________________________________

First figure is the displcement time seris and ssecond is the obtain displcamenet.

Answer 1

[VBBV]

clear all
clc
% Define parameters
n = 3; % Number of degrees of freedom
m = eye(n); % Mass matrix (identity matrix for simplicity)
c = 0.1 * eye(n); % Damping matrix (proportional damping)
k = diag([11, 8, 9]); % Nonlinear stiffness matrix
% Time parameters
dt = 0.001; % Time step
t_final = 10; % Final time
t = 0:dt:t_final;
% Sine loading parameters
omega = 2*pi*0.5; % Frequency of the sine loading
F_sine = @(t) 200*10^3*sin(omega*t); % Sine loading function
% Initial conditions
x0 = zeros(n, 1); % Initial displacement
v0 = zeros(n, 1); % Initial velocity
% Newmark-Beta parameters
beta = 0.25; % Newmark-Beta parameter
gamma = 0.5; % Newmark-Beta parameter
% Initialize arrays
x = zeros(n, length(t)); % Displacement array
v = zeros(n, length(t)); % Velocity array
a = zeros(n, length(t)); % Acceleration array
% Initial conditions
x(:,1) = x0;
v(:,1) = v0;
a(:,1) = m \ (-c*v(:,1) - NonlinearStiffness(x(:,1), k));
% Newmark-Beta iteration
for i = 2:length(t)
    % Predictor step
    x_pred = x(:,i-1) + dt*v(:,i-1) + (dt^2/2) * ((1-2*beta)*a(:,i-1) + 2*beta*a(:,i-1));
    v_pred = v(:,i-1) + dt*((1-gamma)*a(:,i-1));
    
    % Corrector step
    a_pred = m \ (-c*v_pred - NonlinearStiffness(x_pred, k) + F_sine(t(i)));
    x(:,i) = x(:,i-1) + dt*v_pred + (dt^2/2) * (beta*a_pred + (1-beta)*a(:,i-1));
    v(:,i) = v_pred + dt*((1-gamma)*a_pred + dt*gamma*a(:,i));
    
    % Update acceleration for next iteration
    a(:,i) = a_pred;
    
end
% Plot results
  
    plot(t, x(1:3,:));
    xlabel('Time');
    ylabel('Displacement DOF 1 ')
    xlim([0 5]); grid 

% Nonlinear stiffness function
function F = NonlinearStiffness(x, k)
    F = k * x.^3; % Element-wise cubic nonlinear stiffness
end

[VBBV]

There are some changes that need to be done for dispalcement in predictor step and velocity in corrector step according to the Newmark Beta algortihm. See this link for further information

https://en.wikipedia.org/wiki/Newmark-beta_method

[Ram]

Thanks

[Parvesh Deepan]

Can you please solve this problem based on Newmark Beta method (problem is attached herewith).