# Newmark's Method for Nonlinear Systems

323 views (last 30 days)

Show older comments

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

Werkstoffmodell.fnUpdateStateVariables;

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 Comments

### Answers (4)

Christopher Wong
on 6 Jul 2018

Edited: Christopher Wong
on 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')

##### 0 Comments

jiawei jiang
on 20 Aug 2017

##### 0 Comments

victorroda
on 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 Comments

civil s
on 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 Comments

### See Also

### Community Treasure Hunt

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

Start Hunting!