# Newmark's Method for Nonlinear Systems

178 Ansichten (letzte 30 Tage)
Markus Landwehr am 23 Mai 2017
Bearbeitet: zhang am 1 Apr. 2024 um 8:13
Hello everyone,
based on the book "Dynamics of Structures" by Chopra I would like to simulate nonlinear vibrations in Matlab with the Newmark´s method for nonlinear systems. I attached the book chapter where the algorithm (modified Newton-Raphson and Newmark´s-method) are explained.
My current implementation of these algorithm:
function [u, udot, u2dot] = newmark_int_nlin(t,p,u0,udot0,m,k,c,gamma,beta,Werkstoffmodell,Solver)
%Newmark's Method for nonlinear systems
%--------------------------------------------------------------------------
% Integrates a nonlinear 1-DOF system with mass "m", spring stiffness "k", damping
% coeffiecient "c" and nonlinear force "fs", when subjected to an external load P(t).
% Returns the displacement, velocity and acceleration of the system with
% respect to an inertial frame of reference.
%
% SYNTAX
% [u, udot, u2dot] = newmark_int_nlin(t,p,u0,udot0,m,k,c,gamma,beta,Werkstoffmodell,Solver)
%
% INPUT
% [t] : Time Vector [n,1]
% [p] : Externally Applied Load [n,1]
% [u0]: Initial Position [1,1]
% [udot0]: Initial Velocity [1,1]
% [m]: System Mass [1,1]
% [k]: System Stiffness [1,1]
% [c]: System Damping [1,1]
% [gamma]: Newmark coefficient [1,1]
% [beta]: Newmark coefficient [1,1]
% [Werkstoffmodell]: material model parameters
% [Solver]: solver parameters
%
% OUTPUT
% [u]: Displacemente Response [n,1]
% [udot]: Velocity [n,1]
% [u2dot]: Acceleration [n,1]
%
% N = number of time steps
%
% The options include changing the value of the "gamma" and "beta"
% coefficient which appear in the formulation of the method. By default
% these values are set to gamma = 1/2 and alpha = 1/4.
dt = t(2) - t(1); %timestep
a = m/(beta*dt) + gamma*c/beta; %newmark coefficient
b = 0.5*m/beta + dt*(0.5*gamma/beta - 1)*c; %newmark coefficient
TOL = 1e-6; %Tolerance
j_max = 100; %max iterations
dp = diff(p); %external force
u = zeros(length(t),1); udot = u; u2dot = u;
u(1,1) = u0; %initial condition
udot(1,1) = udot0;%initial condition
u2dot(1,1) = 1/m*(p(1)-k*u0-c*udot0-Fresfun(u(1,1),t(1),Werkstoffmodell,Solver)); %initial condition
% u2dot(1) = 1/m*(p(1)-k*u0-c*udot0);
for i = 1:(length(t)-1) %for each timestep
dp_dach = dp(i) + a*udot(i,1) + b*u2dot(i,1);
% ki = (Fresfun(u(i+1,1),t(i),Werkstoffmodell,Solver)-Fresfun(u(i,1),t(i),Werkstoffmodell,Solver))/(u(i+1,1)-u(i,1));%tangent stiffness
ki = k;
ki_dach = ki + gamma/(beta*dt)*c + m/(beta*dt^2); %tangent stiffness
% modified Newton Raphson iterative procedure
% initialize data
j = 2;
u(i+1,1) = u(i,1);
fs(i,1) = Fresfun(u(i,1),t(i),Werkstoffmodell,Solver);
dR(1) = 0;
dR(2) = dp_dach;
kT = ki_dach;
% calculation for each iteration
while j < j_max
du(j) = (dR(j))/kT;
u(i+1,j) = u(i+1,j-1)+du(j);
fs(i,j) = Fresfun(u(i,j),t(i),Werkstoffmodell,Solver); %compute fs(j)
df(i,j) = fs(i,j)-fs(i,j-1)+(kT-k)*du(j);
dR(j+1) = dR(j)-df(i,j);
du_sum = sum(u,2);
if du(j)/du_sum(j) < TOL % repetition for next iteration
break;
end
j = j+1;
end
du(i) = du_sum(j);
dudot_i = gamma/(beta*dt)*du(i) - gamma/beta*udot(i) + dt*(1-0.5*gamma/beta)*u2dot(i);
du2dot_i = 1/(beta*dt^2)*du(i) - 1/(beta*dt)*udot(i) - 0.5/beta*u2dot(i);
u(i+1) = du(i) + u(i);
udot(i+1) = dudot_i + udot(i);
u2dot(i+1) = du2dot_i + u2dot(i);
end
In this state, the code is not working properly. I am having issues with the "modified newton raphson algorithm" and I think there are mistakes in my code.
I would appreciate if you could check my code and give me feedback.
Best regards Markus
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Antworten (4)

Christopher Wong am 6 Jul. 2018
Bearbeitet: Christopher Wong am 24 Sep. 2018
Hi, I've also been attempting these problems from Chopra (2014). I'm having some trouble following your syntax and would need more descriptive %documentation. I have not attempted the constant stiffness (modified) Newton-Raphson method yet, but here is a code I made that works for generic Newton-Raphson. It reproduced the results from Chopra, Example 5.5 exactly.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Newmark's Method for SDF Nonlinear Systems %%%
By: Christopher Wong %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Method as per A.K. Chopra %%%
Example 5.5, (2014) %%%
%%%Using Newton-Raphson iterations
% Clear workspace for each new run
clear
% Establish time step parameters
T_n=1; % natural period
dt=0.1; % time step size
t=0:dt:1.0; % total length of time
n=length(t)-1; % number of time steps
TOL=1e-3; % convergence criteria
% Determine which special case to use: constant avg. vs. linear accel
if dt/T_n<=0.551 % Use linear accel. method - closest to theoretical
gamma=0.5;
beta=1/6;
else % Use constant avg. accel. method - unconditionally stable (Example 5.5)
gamma=0.5;
beta=0.25;
end
%%%Establish system properties
xi=0.05; % percentage of critical damping
omega_n=2*pi/T_n; % natural angular frequency
m=0.4559; % mass
k=18; % stiffness
c=2*xi*m*omega_n; % damping constant
u_y=2; % yield deformation
%%%Input excitation function
p(t<0.6)=50*sin(pi*t(t<0.6)/0.6);
p(t>=0.6)=0;
%%%Establish initial conditions @ i=1
u(1)=0; % displacement
v(1)=0; % velocity
f_s(1)=0; % restoring force
k_T(1)=k; % tangent stiffness
a(1)=(p(1)-c*v(1)-f_s(1))/m; % acceleration
%%%Calculate Newmark constants
A1=m/(beta*dt.^2)+gamma*c/(beta*dt);
A2=m/(beta*dt)+c*(gamma/beta-1);
A3=m*(1/(2*beta)-1)+dt*c*(gamma/(2*beta)-1);
k_hat=k+A1;
%%%Calculations for each time step, i=0,1,2,...,n
%%%Inititialize time step
for i=1:n
p_hat(i+1)=p(i+1)+A1*u(i)+A2*v(i)+A3*a(i); % Chopra eqn. 5.4.12
u(i+1)=p_hat(i+1)/k_hat; % linear displacement
k_T(i+1)=k_T(i); % tangent stifness at beginning of time step
f_s(i+1)=u(i+1)*k_T(i+1); % linear restoring force
%%%Determine if iteration is linear or nonlinear
if abs(f_s(i+1))<=abs(k*u_y) % If linear
u(i+1)=u(i+1); % keep linear value
f_s(i+1)=f_s(i+1); % keep linear value
else % If nonlinear
u(i+1)=u(i); % restore value from previous nonlinear iteration
f_s(i+1)=f_s(i); % resore value from previous nonlinear iteration
end
R_hat(i+1)=p_hat(i+1)-f_s(i+1)-A1*u(i+1); % Compute initial residual
%%%Begin Netwon-Raphson iteration
while abs(R_hat(i+1))>TOL % Terminate if converged
k_T_hat(i+1)=k_T(i+1)+A1;
du=R_hat(i+1)/k_T_hat(i+1);
u(i+1)=u(i+1)+du;
f_s(i+1)=f_s(i)+k*(u(i+1)-u(i));
%%%Determine if restoring force is yielding
if abs(f_s(i+1))>=abs(k*u_y) % If yielding
f_s(i+1)=k*u_y;
k_T(i+1)=0;
else % If elastic
k_T(i+1)=18;
end
R_hat(i+1)=p_hat(i+1)-f_s(i+1)-A1*u(i+1); % Compute new residual
end
%%%Calculations for new velocity and acceleration
v(i+1)=gamma*(u(i+1)-u(i))/(beta*dt)+v(i)*(1-gamma/beta)+dt*a(i)*(1-gamma/(2*beta));
a(i+1)=(u(i+1)-u(i))/(beta*dt.^2)-v(i)/(beta*dt)-a(i)*(1/(2*beta)-1);
end
%%%Generate Solution Table
solution=[t;p;R_hat;k_T;k_T_hat;u;f_s;v;a];
solution=solution';
colNames={'t','p','R_hat','k_T','k_T_hat','u','f_s','v','a'};
solution=array2table(solution,'VariableNames',colNames);
%%%Generate plots
figure(3)
plot(t,u)
xlabel('\itt\rm, s')
ylabel('\itu\rm, cm')
title('Displacement vs. Time Plot')
figure(4)
plot(t,p)
xlabel('\itt\rm, s')
ylabel('\itp\rm, kN')
title('Excitation Function Plot')
figure(5)
plot(u,f_s)
xlabel('\itu\rm, cm')
ylabel('\itf_s\rm, kn')
title('Elastoplastic Loop Plot')
##### 1 Kommentar-1 ältere Kommentare anzeigen-1 ältere Kommentare ausblenden
Shivang Bhatt am 28 Aug. 2022
Bearbeitet: Shivang Bhatt am 28 Aug. 2022
Hey Christopher.
I've a little doubt. I'm using your code to run example 5.5 of A.K.Chopra sir's book and I am not getting the same answers. Can you help me? My guess is the wrong inputs I am putting in the code! You can contact me on 21se801@bvmengineering.ac.in
Thank you!

Melden Sie sich an, um zu kommentieren.

victorroda am 29 Nov. 2018
I have worked out this code, but I find that it does not work properly when the spring is unloading within the yield region.
Additional lines must be implemented in the yield condition to satisfy that if fS·dU<0, then the stiffness is no longer 0, but the stiffness of the linear region.
Regards,
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

civil s am 1 Apr. 2020
hi
my code doesnt converge for negative post yield stiffness.can anyone help please? it is ok for zero and positive post yield stiffness but it does not work for negative one.
pls let me know if anyone can help so i will send the code.
thanks
##### 3 Kommentare1 älteren Kommentar anzeigen1 älteren Kommentar ausblenden
civil s am 1 Apr. 2020
no i do not have any problem with zero and positive post yield stiffness. my problem is the negative one. and i am sure there is a solution for it. but i dont know what the solution is :D
Ya Su am 12 Apr. 2021
Have you found the solution for negative stiffness?

Melden Sie sich an, um zu kommentieren.

zhang am 31 Mär. 2024 um 15:48
Bearbeitet: zhang am 1 Apr. 2024 um 8:13
Here is a step by step implementation of the algorithm in Dr. Chopra's book.
It produces the almost identical result to Example 5.5 and 5.6 except for some rounding error in numerical computing.
The code of example 5-6 differs from example 5-5 only at step 3.3, check the comment.
%% pre-processing
clear;
clc;
% system properties
dampingRatio = 0.05;
mass = 0.2533;
stiffness = 10;
angularFreq = sqrt(stiffness / mass);
damping = 2 * dampingRatio * mass * angularFreq;
yieldDeform = 0.75;
% % timestep and instants
timeStep = 0.1;
time = (0:timeStep:1.0)';
numSteps = numel(time) - 1;
% excitation function
forceFunc = @(t) (t < 0.6) .* 10 .* sin(pi * t / 0.6);
force = forceFunc(time);
% specify method
method = 1;
switch method
case 2 % linear acceleration
gamma = 0.5;
beta = 1/6;
otherwise % average constant acceleration
gamma = 0.5;
beta = 0.25;
end
% criteria of convergence
tolerance = 1e-3;
maxIteration = 10;
% Preallocate arrays for speed
adjustedForce = zeros(numSteps + 1, 1);
displacement = zeros(numSteps + 1, 1);
tangentStiffness = zeros(numSteps + 1, 1);
restoreForce = zeros(numSteps + 1, 1);
residual = zeros(numSteps + 1, maxIteration);
deltaDisplacement = zeros(numSteps + 1, maxIteration);
adjustedTangentStiffness = zeros(numSteps + 1, 1);
velocity = zeros(numSteps + 1, 1);
acceleration = zeros(numSteps + 1, 1);
%% nonlinear solver
% Step 1.0 - initial calculations
% - 1.1
displacement(1) = 0;
velocity(1) = 0;
restoreForce(1) = 0;
tangentStiffness(1) = stiffness;
% - 1.2
acceleration(1) = (force(1) - damping * velocity(1) - restoreForce(1)) / mass;
% - 1.3 // time step has been given
% - 1.4
newmarkConst1 = mass / (beta * timeStep.^2) + gamma * damping / (beta * timeStep);
newmarkConst2 = mass / (beta * timeStep) + damping * (gamma / beta - 1);
newmarkConst3 = mass * (1 / (2 * beta) - 1) + timeStep * damping * (gamma / (2 * beta) - 1);
i = 1;
while i <= numSteps
% Step 2.0 Calculations for each time instant
% - 2.1
displacement(i + 1) = displacement(i);
restoreForce(i + 1) = restoreForce(i);
tangentStiffness(i + 1) = tangentStiffness(i);
% - 2.2
adjustedForce(i + 1) = force(i + 1) + newmarkConst1 * displacement(i) + newmarkConst2 * velocity(i) + newmarkConst3 * acceleration(i);
% Step 3.0 Modified Netwon-Raphson iteration
j = 1;
% - 3.1
residual(i + 1, j) = adjustedForce(i + 1) - restoreForce(i + 1) - newmarkConst1 * displacement(i + 1);
% - 3.2 Check convergence
while (abs(residual(i + 1, j)) > tolerance) && (j < maxIteration)
% - 3.3 // (ex.5-5 Newton-Raphson) update at every iteration
% adjustedTangentStiffness(i + 1, j) = tangentStiffness(i + 1) + newmarkConst1;
% - 3.3 // (ex.5-6 Modifided Newtin-Raphson) update once at the first iteration
if j == 1
adjustedTangentStiffness(i + 1, j) = tangentStiffness(i + 1) + newmarkConst1;
else
end
% - 3.4
deltaDisplacement(i + 1, j) = residual(i + 1, j) / adjustedTangentStiffness(i + 1, j);
% - 3.5
displacement(i + 1) = displacement(i + 1) + deltaDisplacement(i + 1, j);
% - 3.6
restoreForce(i + 1) = restoreForce(i) + stiffness * (displacement(i + 1) - displacement(i));
tangentStiffness(i + 1) = stiffness;
if restoreForce(i + 1) > stiffness * yieldDeform
% yielding in positive direction
restoreForce(i + 1) = stiffness * yieldDeform;
tangentStiffness(i + 1) = 0;
elseif restoreForce(i + 1) < - stiffness * yieldDeform
% yielding in negative direction
restoreForce(i + 1) = - stiffness * yieldDeform;
tangentStiffness(i + 1) = 0;
end
% - 3.7 -> 3.1
j = j + 1;
residual(i + 1, j) = adjustedForce(i + 1) - restoreForce(i + 1) - newmarkConst1 * displacement(i + 1);
end
% 4.0 Calculations for velocity and acceleration
% - 4.1
velocity(i + 1) = gamma * (displacement(i + 1) - displacement(i)) / (beta * timeStep) + velocity(i) * (1 - gamma / beta) + timeStep * acceleration(i) * (1 - gamma / (2 * beta));
% - 4.2
acceleration(i + 1) = (displacement(i + 1) - displacement(i)) / (beta * timeStep.^2) - velocity(i) / (beta * timeStep) - acceleration(i) * (1 / (2 * beta) - 1);
% 5.0 Replace i by i + 1
i = i + 1;
end
%% post-processing
% show solution
solution = [time, force, residual(:, 1), tangentStiffness, adjustedTangentStiffness(:, 1), deltaDisplacement(:, 1), displacement, restoreForce, velocity, acceleration];
columnNames = {'Time', 'Force', 'Residual', 'Tangent_Stiffness', 'Adjusted_Tangent_Stiffness', 'deltaDisplacement', 'Displacement', 'Restoring_Force', 'Velocity', 'Acceleration'};
solution = array2table(solution, 'VariableNames', columnNames);
disp(solution);
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Newton-Raphson Method finden Sie in Help Center und File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by