Clear Filters
Clear Filters

Hello everyone i am solving multi degree freedom of system by using Newmarks beta method in which loading function F(x,v). When i am ruing the code, the result i

8 views (last 30 days)
% dynamic analysis using direct integration method
% % the porbleme is: Mx''+Cx'+Kx=(C1x'+K1x) *eps*cos(Omeg*t)+K2x* eps^2cos(Omeg*t)^2
clc;
clear all
format short e
%close all
m=[18438.6 0 0;0 13761 0;0 0 9174];
disp(' mass matrix')
m;
[ns,ms]=size(m);
% if forces are acting at degrees of freedom
m=[40000 0 0;0 20000 0;0 0 20000];
c0=[0 0 0;0 -8000 8000;0 -8000 8000];
c1=[0 0 0;0 -4000 8000;0 -8000 4000];
k0=[30000 -10000 0;-10000 20000 -10000;0 -10000 10000];
k1=[30000 -10000 0;-10000 50000 -10000;0 -10000 50000];
k2=[30000 -10000 0;-10000 20000 -40000;0 -40000 10000];
% disp(' force at various degrees of freedom');
% f;
%if base ground acceleration is given
% dis='disp.dat'
% di=load(dis)
% % convert to equivalent nodal loads
% for i=1:ns
% f(:,i)=-di*m(i,i)
% end
disp(' damping matrix')
c0;
disp(' stiffness matrix')
k0;
format long;
kim=inv(k0)*m;
[evec,ev]=eig(kim);
for i=1:ns
omega(i)=1/sqrt(ev(i,i));
end
disp(' natural frequencies')
omega;
% give gamma=0.5 and beta=0.25 for Newmark average accln method
%gama=0.5;
%beta=0.25;
%give gamma=0.5 and beta=0.1667 for Newmark linear accln method
gama=0.5;
beta=0.167;
%give initial conditions for displacements
u0=[0 0.01 0.05];
disp(' initial displacements')
u0;
%give initial condition for velocities
v0=[0. 0. 0.];
%y0=[0;0.01;0.05;0;0;0];
disp(' initial velocities')
v0;
om=5; eps=0.01;
IM=inv(m);
X=u0'; Z=v0';
tt=0;
% S1=-IM*k0+IM*k1*eps.*cos(om*tt)+IM*k2*eps^2.*cos(om*tt)^2;
% S2=-IM*c0+IM*c1*eps.*cos(om*tt);
dt=0.02;
S=k2*dt^2*beta*eps^2*cos(om*tt)^2+k1*dt^2*beta*eps*cos(om*tt)+...
c1*gama*dt*eps*cos(om*tt)+k0*dt^2*beta+c0*gama*dt+m;
S1=k1*eps.*cos(om*tt)+k2*eps^2.*cos(om*tt)^2;
S2=c1*eps.*cos(om*tt);
f(1,:)=(S1*X+S2*Z);
%for i=1:ns
a0=-inv(m)*(f(1,:)'+c0*v0'+k0*u0');
%end
kba=k0+(gama/(beta*dt))*c0+(1/(beta*dt*dt))*m;
kin=inv(kba);
aa=(1/(beta*dt))*m+(gama/beta)*c0;
bb=(1/(2.0*beta))*m+dt*(gama/(2.0*beta)-1)*c0;
u(1,:)=u0;
v(1,:)=v0;
a(1,:)=a0;
t=linspace(0,5,251);
for i=2:10%251
X=u(i-1,:)'; Z=v(i-1,:)';
S1=k1*eps.*cos(om*tt)+k2*eps^2.*cos(om*tt)^2;
S2=c1*eps.*cos(om*tt);
f(i,:)=(S1*X+S2*Z); %%%%% ??????
df=f(i,:)-f(i-1,:);
dfb(i,:)=df+v(i-1,:)*aa'+a(i-1,:)*bb';
du(i,:)=dfb(i,:)*kin;
dv(i,:)=(gama/(beta*dt))*du(i,:)-(gama/beta)*v(i-1,:)+dt*...
(1-gama/(2.0*beta))*a(i-1,:);
da(i,:)=(1/(beta*dt^2))*du(i,:)-(1/(beta*dt))*v(i-1,:)-(1/(2.0*beta))*a(i-1,:);
u(i,:)=u(i-1,:)+du(i,:);
v(i,:)=v(i-1,:)+dv(i,:);
a(i,:)=a(i-1,:)+da(i,:);
end
%figure(1);
hold on
plot(tt,u(:,1),'k');
xlabel(' time in secs');
ylabel(' roof displacement');
title(' displacement response of the roof');
%figure(2);
hold on
plot(tt,u(:,2),'k');
xlabel(' time in secs');
ylabel(' roof velocity');
title('velocity response of the roof');
%figure(3);
hold on
plot(tt,u(:,3),'k');
xlabel(' time in secs');
ylabel(' roof acceleration');
title(' acceleration response of the roof')
  1 Comment
Ashok
Ashok on 26 Aug 2024
Edited: Ashok on 26 Aug 2024
Hi abdellhadi,
Your question seems to be incomplete. I observed that the figure is blank. Is that the error you are referring to?

Sign in to comment.

Answers (1)

Zinea
Zinea on 2 Sep 2024
The code needs the following changes to give accurate results for the Newmark-beta method for solving the dynamic response of a multi-degree-of-freedom system with a specific loading function:
  1. The iteration loop goes from i’ = 2 to 10, but the time vector ‘t’ is defined from 0 to 5 seconds with 251 points. The loop needs to be adjusted to start from ‘i’ = 2 to 251 to cover the entire time range.
  2. tt’ must be updated within the loop to reflect the current time step: tt = t(i)
  3. The plotting commands are currently inside the loop and may not properly display the results. They should be moved outside the loop along with the usage of separate figures for displacement, velocity and acceleration.
  4. It would be good to use a vectorized approach to compute results for multiple degrees of freedom simultaneously.
  5. Pre-allocating arrays (u,v,a) for storing results and initializing them with zeros or initial conditions in a single step.
Here is the MATLAB code for your reference:
% Define mass, damping, and stiffness matrices
m = [40000 0 0; 0 20000 0; 0 0 20000];
c0 = [0 0 0; 0 -8000 8000; 0 -8000 8000];
c1 = [0 0 0; 0 -4000 8000; 0 -8000 4000];
k0 = [30000 -10000 0; -10000 20000 -10000; 0 -10000 10000];
k1 = [30000 -10000 0; -10000 50000 -10000; 0 -10000 50000];
k2 = [30000 -10000 0; -10000 20000 -40000; 0 -40000 10000];
% Calculate natural frequencies
kim = inv(k0) * m;
[evec, ev] = eig(kim);
omega = 1 ./ sqrt(diag(ev))'; % Vectorized calculation of natural frequencies
disp('Natural frequencies:');
disp(omega);
% Newmark-beta method parameters
gama = 0.5;
beta = 0.167;
% Initial conditions
u0 = [0 0.01 0.05];
v0 = [0 0 0];
% System parameters
om = 5;
eps = 0.01;
IM = inv(m);
% Time settings
dt = 0.02;
t = linspace(0, 5, 251);
% Pre-allocate response arrays
u = zeros(length(t), length(u0));
v = zeros(length(t), length(v0));
a = zeros(length(t), length(u0));
% Initial conditions
u(1, :) = u0;
v(1, :) = v0;
% Initial acceleration
a(1, :) = (-IM * (c0 * v0' + k0 * u0'))';
% Calculate effective stiffness matrix
kba = k0 + (gama / (beta * dt)) * c0 + (1 / (beta * dt^2)) * m;
kin = inv(kba);
aa = (1 / (beta * dt)) * m + (gama / beta) * c0;
bb = (1 / (2.0 * beta)) * m + dt * (gama / (2.0 * beta) - 1) * c0;
% Precompute constant matrices
S1_base = eps * k1;
S2_base = eps * c1;
S1_quad = eps^2 * k2;
% Time integration loop
for i = 2:length(t)
tt = t(i);
% Update force vector using vectorized operations
cos_om_tt = cos(om * tt);
S1 = S1_base * cos_om_tt + S1_quad * cos_om_tt^2;
S2 = S2_base * cos_om_tt;
% Calculate effective force increment
f = S1 * u(i-1, :)' + S2 * v(i-1, :)';
dfb = f + aa * v(i-1, :)' + bb * a(i-1, :)';
% Newmark-beta updates, vectorized
du = kin * dfb;
dv = (gama / (beta * dt)) * du - (gama / beta) * v(i-1, :)' + dt * (1 - gama / (2.0 * beta)) * a(i-1, :)';
da = (1 / (beta * dt^2)) * du - (1 / (beta * dt)) * v(i-1, :)' - (1 / (2.0 * beta)) * a(i-1, :)';
% Update responses
u(i, :) = u(i-1, :) + du'; % Transpose du to match 1-by-3
v(i, :) = v(i-1, :) + dv';
a(i, :) = a(i-1, :) + da';
end
% Plotting results for displacement (moved outside loop)
figure;
plot(t, u(:, 1), 'k');
xlabel('Time (s)');
ylabel('Displacement');
title('Displacement Response of the First Degree of Freedom');
grid on
% Simlarly plot for velocity and acceleration
Here is the output from the above code and the plotted figure:
You may refer to the following documentation links for more information on vectorization and pre-allocation in MATLAB :
  1. https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html
  2. https://www.mathworks.com/help/matlab/matlab_prog/preallocating-arrays.html
Hope this helps!

Tags

Community Treasure Hunt

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

Start Hunting!