Newmark's Method for Nonlinear Systems

148 ビュー (過去 30 日間)
Markus Landwehr
Markus Landwehr 2017 年 5 月 23 日
編集済み: Shivang Bhatt 2022 年 8 月 28 日
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

回答 (4 件)

Christopher Wong
Christopher Wong 2018 年 7 月 6 日
編集済み: Christopher Wong 2018 年 9 月 24 日
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 件のコメント
Shivang Bhatt
Shivang Bhatt 2022 年 8 月 28 日
編集済み: Shivang Bhatt 2022 年 8 月 28 日
Hey Christopher.
Very impressed by your work.
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!

サインインしてコメントする。


jiawei jiang
jiawei jiang 2017 年 8 月 20 日
hello,I'm researching on the newmark-βmethod,have you resolve this problem ,could you send me a code of your。so appreciate。my email is : jjwsmile@163.com

victorroda
victorroda 2018 年 11 月 29 日
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,

civil s
civil s 2020 年 4 月 1 日
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 件のコメント
civil s
civil s 2020 年 4 月 1 日
thanks for your kind response.
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
Ya Su 2021 年 4 月 12 日
Have you found the solution for negative stiffness?

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeNewton-Raphson Method についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by