In the following program replace the function of duhamel3 with a Newmark beta me

Question

In the following program replace the function of duhamel3 with a Newmark beta method function, the results of the program must be the same as using the duhamel3 function. Please provide the code of the Newmark beta method in matlab and a brief explanation.

%---------------------- Program VIBRA1.m ----------------------------%
%                                                                    %
% Program to calculate the response of a viscously underdamped single%
% DOF system to an arbitrary force F(t). using using recursive       %
% equation.                                                          %
%                                                                    %
%------------------ Last updated: November 12 -2017 -----------------%
clear all; close all

%----------- System properties --------------------------------------%

m = 0.85;                                  % mass of the system (N.m.s^2/rad)
c = 0.9;                                   % damping coefficient (N.m.s/rad)  
k = 140.0;                                 % stiffness coefficient (N.m/rad)
u0 = 0.0;                                  % initial angular displacement (rad)
v0 = 0.0;                                  % initial angular velocity (rad/s)
dt = 0.0002;                               % time interval of the external moment (sec)
F0 = 4;                                    % amplitude of the moment (N.m)
tf = 10;                                   % final time(sec)


%-------- Read and plot force time history ------------------------%

t = (0:dt:tf)';                       
nt = length(t);
for ii =1:nt
   ft(ii)=0;
end
for ii = 1:10000
   ft(ii) = F0*sin(3.14/2*t(ii));
end
figure; plot( t,    ft ); grid on;
xlabel('Time [sec]'); ylabel('Moment [N.m]')
   
  
%---------- Calculate the natural frequency, natural period and damping ratio -----------%
   
wn = sqrt( k / m );    
Tj = 2*pi ./ wn;
zj = c / (2* m * wn);
wd = wn *sqrt(1-zj^2);

disp(['*** The natural frequency [rad/s] is :',num2str(wn)]);
disp(['*** The natural period [sec] is :',num2str(Tj)]);
disp(['*** The damping ratio is :',num2str(zj)]);
disp(['*** The damped natural frequency [rad/s] is :',num2str(wd)]);

%---------- Calculate the modal and physical response -----------

[u , v , a] = duhamel3(wn,zj,m,dt,u0,v0,ft,t);          % Duhamel variation linear force



%---------- Plot the response time histories --------------------

disp(['*** The maximum rotation [rad.] is :',num2str(max(u))]);

figure; plot( t,u, 'MarkerSize',10 ); grid on
legend('Angular Displacement')
title('Angular Displacement time history ');
ylabel('Angular displacement [rad]'); xlabel('Time [sec]')

figure; plot( t,v, 'MarkerSize',10 ); grid on
legend('Angular Velocity')
title('Angular velocity time history ');
ylabel('Angular velocity [rad/s]'); xlabel('Time [sec]')

figure; plot( t,a, 'MarkerSize',10 ); grid on
legend('Angular Acceleration')
title('Angular acceleration time history ');
ylabel('Angular acceleration [rad/s^2]'); xlabel('Time [sec]')

-------------------------------------------------------------------------------------------------

function [u,v,a] = duhamel3( wn, z, m, dt, u0, v0, f, t )
% Program to calculate the response of a linear viscously damped
% SDOF system subjected to an excitation f(t) with arbitrary time
% variation sampled at constant intervals of duration 'dt'.
%
% The excitation is assumed to be constant in each time interval.
%
% d^2(q(t))/dt^2 + 2*z*w*d(q(t))/dt + (w^2)*q(t) = f(t)
%
% d^2(q(t))/dt^2 = f(t) - 2*z*w*d(q(t))/dt - (w^2)*q(t)
%
% Input data:
% ----------
% w = Natural frequency [rad/sec].
% z = Damping ratio for underdamped case.
% m = Mass of oscillator [consistent units].
% dt = Time interval [sec].
% u0 = Initial displacement [Length].
% v0 = Initial velocity [Length/sec].
% f = Vector with the sampled excitation [consistent units].
% t = Vector with discrete time [sec].
%

% Output data:
% -----------
% u = Displacement time history.
% v = Velocity time history.
% a = Acceleration time history.

% Initialize:
nt = length(t);
Z = zeros(3,nt);

% Initial conditions:

Z(1,1) = u0;
Z(2,1) = v0;
Z(3,1) = ( f(1) / m ) - (wn^2)*u0 - 2*z*wn*v0;

% Calculate coefficients:

wd = wn * sqrt(1-z^2);
ex = exp(-z*wn*dt);
co = cos(wd*dt);
si = sin(wd*dt);

a11 = ( co + z/sqrt(1-z^2) * si ) * ex;
a12 = 1 / wd * si * ex;
a21 = -wn^2 * a12;
a22 = ( co - z/sqrt(1-z^2) * si ) * ex;
SS = sqrt(1-z^2);
kas = 1/(m*dt*SS*wn^3);
b11 = kas*(2*z*SS+((1-2*z^2-z*wn*dt)*si-(2*z*SS+wd*dt)*co)*ex);
b12 = kas*(wd*dt-2*z*SS+(-(1-2*z^2)*si+2*z*SS*co)*ex);
b21 = kas*(-wd+(wn*(z+wn*dt)*si+wd*co)*ex);
b22 = kas*(wd-(wn*z*si+wd*co)*ex);

A1 = [a11 , a12 ; a21 , a22] ;
B1 = [b11 , b12 ; b21 , b22] ;

% Compute response:

for i = 2 : nt-1
    
     ff = [f(i-1) ; f(i)];
    
     Z(1:2,i) = A1 * Z(1:2,i-1) + B1 * ff;

     Z(3,i)   = ( f(i) / m ) - (wn^2)*Z(1,i) - 2*z*wn*Z(2,i);

end

u = Z(1,:);
v = Z(2,:);
a = Z(3,:);

return

Explanation / Answer

function[At,V,D]=newm2(xm,hvec,s,beta,z,dt,ndev)

%newmark beta for MDOF

%xm: mass matrix

%hvec: each floor damping [h(n);h(n-1);...]

%s: stiffness matrix

%beta

%z: ground acceleration

%dt: time interval

%ndev: time interval division

siz=length(z);

nsiz=length(xm);

a(:,1)=zeros(nsiz,1)-z(1,1);

v(:,1)=zeros(nsiz,1);

d(:,1)=zeros(nsiz,1);

sa=zeros(nsiz,1);

sv=zeros(nsiz,1);

sd=zeros(nsiz,1);

temp=inv(xm)*s;

[vec,ei]=eig(temp);

h=zeros(nsiz,nsiz);

h=h+diag(hvec);

c=2*h*sqrt(ei)*xm;

for j=2:length(z)

dz=z(j,1)-z(j-1,1);

ddt=dt/ndev;

ddz=dz/ndev;

for k=1:ndev

% Newmark beta integration

ap=-xm*eye(nsiz,1)*ddz+xm*(v/(beta*ddt)+a/(2*beta))+c*(v/(2*beta)+(.25/beta-1)*a*ddt);

ak=s+c/(2*beta*ddt)+1/(beta*ddt*ddt)*xm;

dd=inv(ak)*ap;

dv=dd/(2*beta*ddt)-v/(2*beta)-(0.25/beta-1)*a*ddt;

da=dd/(beta*ddt*ddt)-v/(beta*ddt)-a/(2*beta);

%=============================

zz=z(j)+ddz*k;

d=d+dd;

v=v+dv;

a=a+da;

aa=a+zz;

At((j-1)*ndev+k,:)=aa';

D((j-1)*ndev+k,:)=d';

V((j-1)*ndev+k,:)=v';

end

end

Related Questions

Navigate

• Browse (All)
• Browse I
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.