Dynamic Shear Building Model

Simulating the dynamics of a shear building using MATLAB

This assignment requires MATLAB programming, so you need to type your source codes and detailed workings and include final results (including plots if required) in each question.

Introduction

This assignment is a small project which gives you an opportunity to work on a small-scale engineering design problem in MATLAB. The engineering system that you will be working on is a shear building model, which is commonly simplified as a spring-damper-mass model with multiple degrees of freedom. The engineering design problem is to find the natural frequencies and choose the right parameters using MATLAB programming. The first step is to write ODEs describing the behaviors of the MDOF model and then write a MATLAB program to model this system. You will then use the simulation results to explore the response of the building under an earthquake and find the best possible parameters.

The model and the parameters

Figure 1 is a schematic of the four story shear building.  Each floor is represented by its mass and the superstructure supporting each floor is idealized by a spring constant representing resistance to lateral motion and a damping coefficient providing frictional energy losses. External forces applied to each floor are ignored in this model, but the base (or the ground) may move during an earthquake. Hence, only four degrees of freedom are needed to describe total displacements of the structure, described as Figure 2.

Figure 1: A four-story shear building excited by the base motion.

Figure 2: A simplified four DOF model for a shear building.

This assignment consists of five main tasks below.

Task a: If the analysis is limited to horizontal motion of the structure, derive the general ODEs for this system. Assume there are no external forces applied to each mass. Your equations should include displacement of the base y(t) (You may need to look at section 9.5 in week 9 or the 2016 assignment for mechanical students) and damping factors.  Then re-write your second-order ODEs into a group of first-order ODEs.

Task b:  Develop a general code named four_dof.m  which includes the base motion and damping. Use the function ode45. Follow the example of  two_dof.m.

Use the code developed in Task b to simulate the dynamics of this structure from t = 0s to t = 6s. The initial conditions are,,.,,and . Other parameters are,,,,,,and. Ignore the base motion and damping in this calculation. Present your results as time-series plots of displacement of each floor in the same plot.

Task c: Write your differential equations in a matrix form. Develop a code named four_dof_matraix.m to find the natural frequencies of the building if the damping between different floors are ignored. Use the same data for the masses and stiffness in Task b. Follow the example of  two_dof_matraix.m.

Print the corresponding eigenvectors and natural circular frequencies of vibrations. Show the four mode shapes of vibration and associated natural frequencies of vibration.

Task d: In a region, when an earthquake occurs, the ground motion can be simply modelledby , where is the circular frequency of the ground motion. The typical  in that region is 30 and amplitude is . Determine the response of the system if the initial displacements and initial velocity of each floor are zero. Assume ,,and. You need to modify the code in Task b.  Change  values and find out what will happen if  values are close to the natural frequencies.

Task e: Base isolation is a technique developed to prevent or minimise damage to buildings during an earthquake. It might be thought that more rigid attachment of a building to its foundation will result in less damage in an earthquake. However, if the foundation is rigidly attached to the building or any other structure, all of the earthquake forces will be transmitted directly to the rest of the building without a change in frequency. If the earthquake has natural frequencies with high energy that match the natural frequencies of the building, it will cause the building to oscillate violently in harmony with the earthquake frequency. However, if the natural frequency of the building can be changed to a frequency that does not coincide with that of earthquakes, the building is less likely to fail. The base isolator reduces the stiffness of the structure and thereby lowers its natural frequency. In this condition, the building’s superstructure will respond to the vibrations as a rigid unit instead of resonating with the vibrations.

Suppose you are designing this building using the base isolation technique, and you have options for stiffness value  from the range of and damping valuefrom the range of . Use  = 30 and. Other parameters are not changed. What parameters and will you choose?

There may be different criteria for the best design. Here in this problem, just simply take the maximum displacement difference between the 4th and 3rd floor during certain time period () as a rough indicator. So for each group of parameters, develop the algorithm to determine . Then compare these values among the different groups of parameters. The best design (the best combination of parameters) gives the minimum during an earthquake. You may need multiple loops to circulate the parameters. Discuss your results.

Task a: Deriving differential equations for the system

Second-order equations;

First-order equations.

Task b:Writing a MATLAB codefour_dof.m  using ode4

four_dof.m;

rhsfunctionand base motion;

Simulating the motions of each floor.

The code of four_dof_matraix.m;

Calculation of natural frequencies;

Plotting four modal shapes of vibration.

Task d:  Studying the response of the building under an earthquakes

Calculating and plotting the response under an earthquakes;

Resonance studies.

Task e: Finding the best parameters

Parameter scan;

Discussions.

Solution

four_dof.m

% four_dof

functionfour_dof

t = 0:0.01:6 ; % time scale

c1 = 0 ;

c2 = 0 ;

c3 = 0 ;

c4 = 0 ;

k1 = 3400000 ;

k2 = 2800000 ;

k3 = 2300000 ;

k4 = 1800000 ;

m1 = 3900 ;

m2 = 3600 ;

m3 = 3500 ;

m4 = 2500 ;

init_x1 = 0 ;

init_v1 = 0.012 ;

init_x2 = 0;

init_v2 = 0.032 ;

init_x3 = 0 ;

init_v3 = 0.055 ;

init_x4 = 0 ;

init_v4 = 0.071 ;

init_y1 = init_x1 ;

init_y2 = init_v1 ;

init_y3 = init_x2 ;

init_y4 = init_v2 ;

init_y5 = init_x3 ;

init_y6 = init_v3 ;

init_y7 = init_x4 ;

init_y8 = init_v4 ;

initial = [init_y1 init_y2 init_y3 init_y4 …

init_y5 init_y6 init_y7 init_y8];

[t,y]=ode45(@rhs,t,initial);

plot(t,y(:,1), t,y(:,3),t,y(:,5), t,y(:,7));

xlabel(‘t’); ylabel(‘x’) ;

legend(‘x1 num’,’x2 num’,’x3 num’, ‘x4 num’);

functiondydt=rhs(t,y)

dydt_1 = y(2);

dydt_2 = -(c1+c2)*y(2)/m1 +c2*y(4)/m1-(k1+k2)*y(1)/m1+k2*y(3)/m1 …

+k1*y_base(t)/m1+c1*y_base_dif(t)/m1;

dydt_3 = y(4);

dydt_4 = c2*y(2)/m2-(c2+c3)*y(4)/m2+k2*y(1)/m2-(k2+k3)*y(3)/m2 …

+k3*y(5)/m2+c3*y(6)/m2 ;

dydt_5 = y(6);

dydt_6 = c3*y(4)/m3-(c3+c4)*y(6)/m3+k3*y(3)/m3-(k3+k4)*y(5)/m3 …

+k4*y(7)/m3+c4*y(8)/m3 ;

dydt_7 = y(8);

dydt_8 = k4*y(5)/m4- k4*y(7)/m4 +c4*y(6)/m4- c4*y(8)/m4 ;

dydt=[dydt_1;dydt_2;dydt_3;dydt_4;dydt_5;dydt_6;dydt_7;dydt_8];

end

functiony_t = y_base(t)

y_t = 0;

end

functiony_t_dif = y_base_dif(t)

y_t_dif = 0;

end

end

four_dof_matrix.m

% four_dof_matrix

functionfour_dof_matrix

k1 = 3400000 ;

k2 = 2800000 ;

k3 = 2300000 ;

k4 = 1800000 ;

m1 = 3900 ;

m2 = 3600 ;

m3 = 3500 ;

m4 = 2000 ;

M = [m1 0 0 0; 0 m2 0 0; 0 0 m3 0; 0 0 0 m4];

K = [k1+k2 -k2 0 0; -k2 k2+k3 -k3 0; 0 -k3 k3+k4 -k4; 0 0 -k4 k4];

B = inv(M)*K ;

[v d]

= eig (B);

omega1 = sqrt(d(1,1)) ;

omega2 = sqrt(d(2,2)) ;

omega3 = sqrt(d(3,3)) ;

omega4 = sqrt(d(4,4)) ;

fprintf(‘ the 1st natural frequency is %12.7f\n’,omega1) ;

fprintf(‘ the 2nd natural frequency is %12.7f\n’,omega2) ;

fprintf(‘ the 3rd natural frequency is %12.7f\n’,omega3) ;

fprintf(‘ the 4th natural frequency is %12.7f\n’,omega4) ;

fprintf(‘ the eigenvactor matrix is: \n’) ;

disp(v) ;

fprintf(‘ the eigenvalue matrix is: \n’) ;

disp(d) ;

end

plot_mode_shapes.m

x1 = [0 -0.5457 0.6688 -0.4380 0.2511] ;

x2 = [0 -0.4916 0.0655 0.5709 -0.6543] ;

x3 = [0 -0.4697 -0.5148 0.0771 0.7131] ;

x4 = [0 0.1988 0.4069 0.5839 0.6737] ;

y = [0 1 2 3 4] ;

plot (x1,y,’o–‘,x2,y,’o-.’,x3,y, ‘o-‘,x4,y, ‘o-‘) ;

legend(‘mode 1’, ‘mode 2’, ‘mode 3’, ‘mode 4’);

xlabel(‘amplitude’); ylabel (‘mass position’) ;

% four_dof

t = 0:0.01:6 ; % time scale

c1 = 0.13 ;

c2 = 0.15 ;

c3 = 0.14 ;

c4 = 0.18 ;

k1 = 3400000 ;

k2 = 2800000 ;

k3 = 2300000 ;

k4 = 1800000 ;

m1 = 3900 ;

m2 = 3600 ;

m3 = 3500 ;

m4 = 2500 ;

omega = 30;

y0 = 0.15;

init_x1 = 0 ;

init_v1 = 0 ;

init_x2 = 0 ;

init_v2 = 0 ;

init_x3 = 0 ;

init_v3 = 0 ;

init_x4 = 0 ;

init_v4 = 0 ;

init_y1 = init_x1 ;

init_y2 = init_v1 ;

init_y3 = init_x2 ;

init_y4 = init_v2 ;

init_y5 = init_x3 ;

init_y6 = init_v3 ;

init_y7 = init_x4 ;

init_y8 = init_v4 ;

initial = [init_y1 init_y2 init_y3 init_y4 …

init_y5 init_y6 init_y7 init_y8];

[t,y]=ode45(@rhs,t,initial);

plot(t,y(:,1), t,y(:,3),t,y(:,5), t,y(:,7));

xlabel(‘t’); ylabel(‘x’) ;

legend(‘x1 num’,’x2 num’,’x3 num’, ‘x4 num’);

functiondydt=rhs(t,y)

dydt_1 = y(2);

dydt_2 = -(c1+c2)*y(2)/m1 +c2*y(4)/m1-(k1+k2)*y(1)/m1+k2*y(3)/m1 …

+k1*y_base(t)/m1+c1*y_base_dif(t)/m1;

dydt_3 = y(4);

dydt_4 = c2*y(2)/m2-(c2+c3)*y(4)/m2+k2*y(1)/m2-(k2+k3)*y(3)/m2 …

+k3*y(5)/m2+c3*y(6)/m2 ;

dydt_5 = y(6);

dydt_6 = c3*y(4)/m3-(c3+c4)*y(6)/m3+k3*y(3)/m3-(k3+k4)*y(5)/m3 …

+k4*y(7)/m3+c4*y(8)/m3 ;

dydt_7 = y(8);

dydt_8 = k4*y(5)/m4- k4*y(7)/m4 +c4*y(6)/m4- c4*y(8)/m4 ;

dydt=[dydt_1;dydt_2;dydt_3;dydt_4;dydt_5;dydt_6;dydt_7;dydt_8];

end

functiony_t = y_base(t)

y_t = y0*sin(omega*t);

end

functiony_t_dif = y_base_dif(t)

y_t_dif = y0*omega*cos(omega*t);

end

end

% four_dof

t = 0:0.01:6 ; % time scale

c2 = 0.15 ;

c3 = 0.14 ;

c4 = 0.18 ;

k2 = 2800000 ;

k3 = 2300000 ;

k4 = 1800000 ;

m1 = 3900 ;

m2 = 3600 ;

m3 = 3500 ;

m4 = 2500 ;

omega = 30;

y0 = 0.15;

init_x1 = 0 ;

init_v1 = 0 ;

init_x2 = 0 ;

init_v2 = 0 ;

init_x3 = 0 ;

init_v3 = 0 ;

init_x4 = 0 ;

init_v4 = 0 ;

init_y1 = init_x1 ;

init_y2 = init_v1 ;

init_y3 = init_x2 ;

init_y4 = init_v2 ;

init_y5 = init_x3 ;

init_y6 = init_v3 ;

init_y7 = init_x4 ;

init_y8 = init_v4 ;

initial = [init_y1 init_y2 init_y3 init_y4 …

init_y5 init_y6 init_y7 init_y8];

min_x43m = 1e10 ;

step_size_k = 200000 ;

step_size_c = 0.1 ;

for k1 = 20000:step_size_k:4000000

for c1 = 0.1:step_size_c:10

[t,y]=ode45(@rhs,t,initial);

x43m = max(abs(y(:,7) – y(:,5))) ;

if x43m < min_x43m

min_x43m = x43m ;

k1_best = k1 ;

c1_best = c1 ;

end

end

end

fprintf(‘ the best k1 is %12.4f\n’,k1_best) ;

fprintf(‘ the best c1 is %12.4f\n’,c1_best) ;

fprintf(‘ the least disp of the 4th floor is %12.4f\n’,min_x43m) ;

functiondydt=rhs(t,y)

dydt_1 = y(2);

dydt_2 = -(c1+c2)*y(2)/m1 +c2*y(4)/m1-(k1+k2)*y(1)/m1+k2*y(3)/m1 …

+k1*y_base(t)/m1+c1*y_base_dif(t)/m1;

dydt_3 = y(4);

dydt_4 = c2*y(2)/m2-(c2+c3)*y(4)/m2+k2*y(1)/m2-(k2+k3)*y(3)/m2 …

+k3*y(5)/m2+c3*y(6)/m2 ;

dydt_5 = y(6);

dydt_6 = c3*y(4)/m3-(c3+c4)*y(6)/m3+k3*y(3)/m3-(k3+k4)*y(5)/m3 …

+k4*y(7)/m3+c4*y(8)/m3 ;

dydt_7 = y(8);

dydt_8 = k4*y(5)/m4- k4*y(7)/m4 +c4*y(6)/m4- c4*y(8)/m4 ;

dydt=[dydt_1;dydt_2;dydt_3;dydt_4;dydt_5;dydt_6;dydt_7;dydt_8];

end

functiony_t = y_base(t)

y_t = y0*sin(omega*t);

end

functiony_t_dif = y_base_dif(t)

y_t_dif = y0*omega*cos(omega*t);

end

end