How the Function "NewmarkNon" is working and meaning of Keyi and Keyp in this coding?
3 views (last 30 days)
Show older comments
%----------------------- Mario Paz Example 7.1 ----------------------------
clear;
%--------------------------------------------------------------------------
inpf=input('File name of input file data ','s');
M=importdata([inpf,'.xlsx']);
f=M.data.Force1; % Force
t = f(:,1); %Time
F = f(:,2); %Force
Dt = t(2)-t(1); %Time interval %Time interval
m=0.2; % Mass
k =12.35; %Stiffness
c =0.274; %Damping coefficient
xi =c/(2*sqrt(m*k)); %Damping ratio
omega = sqrt(k/m); %Natural frequency (rad/sec)
Rt = 15; %Forces yielding in tension
Rc = -15; %Forces yielding in compression
%%%%Linear acceleration method (Newmark beta method[Ch.6])
gamma =1/2;
beta = 1/6;
tt= length(t);
%%%Initial calculation
u(1)=0; %Initial condition; Displ.
v(1)=0; %Initial condition; Velocity
a(1)=(F(1)-c*v(1)-k*u(1))/m;
A = m/(beta*Dt)+gamma*c/beta; %A in DFbar = DF + A*v0+B*a0 (Eq. 6.46)
B = m/(2*beta)+Dt*c*((0.5*gamma/beta)-1); %B in DFbar = DF + A*v0+B*a0 (Eq. 6.46)
kbar = k +gamma*c/(beta*Dt)+m/(beta*Dt*Dt); %Eq.6.45
key = 0; %Initial key=0 before Loop over
k_p = k; %Initial stiffness before iteration
%%%Setting up for initial value of Loop over
u0=u;
v0=v;
a0=a;
t =t(1);
%%%% Calculate initial yield points %%%%
ut = Rt/k;
uc = Rc/k;
R(1)=0;
%%%%Iteration for each time step using Newmark beta method
ua =[];
va =[];
aa=[];
ta=[];
for i = 1:(tt-1)
DF=F(i+1)-F(i);
F1 = F(i+1);
[t,u,v,a, kbar, R, keyp, key, Du, k_p] = NewmarkNon( t, DF, Dt, u0, v0, ut, uc, a0, F1, k,c, m, Rt, Rc, R, gamma, beta, key, k_p);
keyi(:,i)=key(:,1)';
keypi(:,i)=keyp(:,1)';
F1i(:,i)=F1(:,1)';
Ri(:,i)= R(:,1)';
k_pi(:,i)=k_p(1,:)';
% %%%Creating column of time,displ,velocity and acceleration
ta = [ta; t];
ua = [ua; u];
va = [va; v];
aa = [aa; a];
%%%New yielding point-Increasing displacement (Eq.6.49)
if v < 0 && v0 >0
ut = max(ua);
uc = ut-(Rt-Rc)/k;
else
ut = ut;
uc = uc;
end
%%%New yielding point-Decreasing displacement (Eq.6.50)
if v > 0 && v0 <0
uc = min(ua);
ut = uc+(Rt-Rc)/k;
else
ut = ut;
uc = uc;
end
%%%Setting up parameters for next iteration
u0 = u;
v0 = v;
a0 = a;
t = t;
R = R;
end
result = [ta,F1i',ua,keypi',keyi',va,Ri',aa,k_pi'];
%%%Plot response
figure (1)
plot (ta, ua);
grid on
xlabel ('t(sec)');
ylabel ('Displacement(in.)');
figure (2)
plot (ua,Ri);
grid on
xlabel ('Displacement(in.)');
ylabel ('Restoring Force (kip)');
function [t,u,v,a, kbar, R, keyp, key, Du, k_p] = NewmarkNon( t, DF, Dt, u0, v0, ut, uc,a0, F1, k, c, m, Rt, Rc, R, gamma, beta, key, k_p)
%%%Algorithm-For each time step:(3),(4),(5),(6),(7)
kbar = k_p +gamma*c/(beta*Dt)+m/(beta*Dt*Dt); %Eq.6.57
A = m/(beta*Dt)+gamma*c/beta; %A in DFbar (Eq. 6.46)
B = m/(2*beta)+Dt*c*((0.5*gamma/beta)-1); %B in DFbar (Eq. 6.46)
DFbar = DF + A*v0+B*a0; % I n cremental effective force (Eq.6.59)
Du = DFbar/kbar; % I n cremental displacement (Eq.6.60)
Dudot = gamma*Du/(beta*Dt)-gamma*v0/beta+ Dt*a0*(1-0.5*gamma/beta); %Incremental velocity(Eq.6.61)
u=u0+Du ; %Displacement at the end of time interval (Eq.6.62)
v=v0+Dudot; %Velocity at the end of time interval (Eq.6.63)
%%%Algorithm-For each time step:(2)-(a)
if u<ut && u>uc
keyp = 0;
elseif u > ut
keyp = 1;
else
keyp = -1;
end
%%%Algorithm-For each time step:(2)-(b) and (2)-(c)
if keyp == 1 && v > 0
key = 1;
elseif keyp == -1 && v < 0
key = -1;
else
key = 0;
end
%%%Algorithm-For each time step:(8) (Eq. 6.65)
if key ==0
if (R+(Du)*k)>=0
R = min(R+(Du)*k, Rt);
else
R = max(R+(Du)*k, Rc);
end
elseif key ==1
R = Rt;
else
R = Rc;
end
%%%Algorithm-For each time step:(3) (Eq.6.58)
if key == 0
k_p = k;
else
k_p = 0;
end
%%%Algorithm-For each time step:(8) (Eq.6.64)
a=1/m*(F1-c*v-R);
t=t+Dt;
end
0 Comments
Answers (0)
See Also
Categories
Find more on Numerical Integration and Differential Equations in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!