Home | | **Structural Dynamics and Earthquake Engineering** | Free vibration of multiple degrees of freedom

Free vibration of multiple degrees of freedom in relation to structural dynamics during earthquakes:
Abstract: In this chapter, free and forced vibration of the multiple-degrees-of-freedom (MDOF) system are discussed. The generalized equations are obtained by applying either Newtonâ€™s second law or solving Lagrange equations. The orthogonality principle between modes is illustrated. It is also shown how modes are normalized. The stiffness method is used to derive the stiffness matrix of a structure. Relevant programs in MATHEMATICA and MATLAB are also given.

**Free vibration of multiple degrees of freedom in relation to structural dynamics during earthquakes**

**Abstract: **In this chapter, free and forced vibration of the multiple-degrees-of-freedom (MDOF) system are discussed. The generalized equations are obtained by applying either Newtonâ€™s second law or solving Lagrange equations. The orthogonality principle between modes is illustrated. It is also shown how modes are normalized. The stiffness method is used to derive the stiffness matrix of a structure. Relevant programs in MATHEMATICA and MATLAB are also given.

**Key words: **Cordonâ€™s solution, orthogonality, normalization, stiffness** **method, static condensation.

**Introduction**

After having seen the behaviour of two degrees of freedom, we now proceed to a more general treatment of the problem. Analysis of two degrees of freedom is more involved than a single degree of freedom. When more masses are considered, the mathematical formulation becomes much more complicated. Most engineering systems are continuous and have an infinite number of degrees of freedom. The vibration analysis of continuous systems requires the solution of partial differential equations which is quite difficult. In fact a solution is not available for many partial differential equations. The analysis of multiple degrees of freedom (MDOF), on the other hand, requires the solution of a set of ordinary differential equation which is relatively simple. Hence, for simplicity of analysis, continuous systems are often approximated as MDOF. All the concepts introduced in the previous chapters can be directly extended to the case of MDOF. One can write one differential equation of equilibrium for each generalized coordinate. The equations can be obtained either by applying Newtonâ€™s second law or solving Lagrangeâ€™s equations.

There are *N* natural frequencies, each associated with its own mode shape for a system having *n* degrees of freedom. Similar to two degrees of freedom, natural frequencies can be obtained by finding the roots of the nonlinear equations. Two degrees of freedom involves quadratic equations in terms of *Ï‰* n2 . Since all the frequencies are real numbers, a direct solution such as *Cordonâ€™s solution *is available for solving a cubic equation. However, as the number of degrees of freedom increases, the solution of the characteristic equation becomes complex. The entire mode shapes exhibit a property known as *orthogonal* property which often enables us to simplify the analysis of such degrees of freedom systems.

**Modelling of a continuous system as an MDOF system**

Consider a three storeyed reinforced concrete (RC) building as shown in Fig. 10.1. One can replace the distributed mass or inertia of the system by a finite number of lumped masses or rigid bodies. The masses are assumed to be connected by mass-less elastic damping members. Linear or angular coordinates are used to describe the motion of the lumped masses. Such models are called lumped parameter, lumped mass or discrete mass systems. The minimum number of coordinates necessary to describe the deformed shape of the structure is called the *generalized degrees of freedom*. Considering the three storeyed RC building as shown in Fig. 10.1a suggests a three degree lumped-mass model as indicated in Fig. 10.1b.

Another popular method of approximating a continuous system as an MDOF system involves replacing the geometry of the system by a large number of elements. By assuming a simple solution within each element, the principle of compatibility and equilibrium are used to find the approximate solution to the original system. This method is the finite element method, which will be discussed later.

**Equations of motion of an MDOF system**

Consider a three storeyed building modelled as shown in Fig. 10.1b. If one wants to study the dynamic behaviour of a multi-storey frame in the *x* direction one can consider the generalized coordinates *U*1, *U*2, *U*3 as shown.

Consider the free body diagram of three masses as shown in Fig. 10.2 and applying Dâ€™Alembertâ€™s principle, the equations of motion are written as

where [*m*] = mass matrix [*c*] = damping matrix [*k*] = stiffness matrix

If the system has *n* degrees of freedom, the size of [*m*], [*c*] and [*k*] is *n *Ã-* n*.

â€¢ **Mass matrix: **There are two ways of forming the mass matrix of the** **structure: lumped mass and consistent mass. The mass matrix shown in Eq. 10.2 is a lumped mass matrix. If we consider rotational degrees of freedom, *rotary inertia* must be considered. The consistent mass matrix is *fully populated* whereas lumped mass matrix is *diagonal.*

â€¢ **Stiffness matrix: **The stiffness matrix [*k*] is symmetric. The total stiffness** **matrix will be formed by assembling all the elements together to form a structure. The stiffness matrix for dynamic analysis is to be calculated by using standard structural application procedures.

â€¢ **Damping matrix: **The damping values are to be obtained experimentally.** **Usually a concrete structure has 5% of critical damping and steel structure and 2% of critical damping. [*c*] matrix will be reduced to simpler forms to facilitate the analysis.

â€¢ **Load vector: **The dynamic loads are assumed to act at nodal points.

**Free undamped vibration of an MDOF system**

The equation of motion for an undamped system can be written as

The determinant is known as frequency determinant or characteristic equation. Equation 10.7 can also be written as

Hence for a non-trivial solution to exist the determinant of the following matrix must be equal to zero

Equation 10.8b leads to the *n*th order nonlinear equation in terms of *Ï‰* *n*2 whereas Eq. 10.10 involves the *n*th order nonlinear equation in terms of 1/*Ï‰* *n*2 . Equation 10.9 is also written as

where [*D*] = [*k*]-1[*m*] = [*a*][*m*] is known as the dynamic matrix. Equation 10.12 is in the form of a typical eigenvalue problem. The vector {*Ïˆ*} is called an eigenvector and *Î»* = 1/ *Ï‰* *n*2 is called the eigenvalue.

One of the various ways of obtain the solution is to make the following determinant zero.

The determinant with respect to eigenvalue is plotted in Fig. 10.3.

The solution of the polynomial equation Eq. 10.14 is known as characteristic equation or frequency equation that will yield *n* values of *Î»*. Once *Î»* are

determined {*Ïˆ*} can be determined. It may be noted that a unique solution for {*Ïˆ*} does not exist. We obtain only the ratios among {*Ïˆ*}s.

**Orthogonality relationship**

For an MDOF system if *Ï‰* *n*( *r* ) , *Ï‰* *n*( *s* ) are two natural frequencies, {*Ïˆ*}(*r*), {*Ïˆ*}(*s*) are the corresponding modal vectors.

For the *r*th mode Eq. 10.7 can be written as

Equation 10.27 can then be interpreted as the work done by inertia forces occurring in the first mode in going through displacements of the second mode are equal to zero. This is known as the *orthogonality* relationship.

**Normalization of modes**

It has been seen that the normal modes indicate the ratio between displacements. As such, the different elements may be varied in such a way that a constant ratio is maintained. There are an infinite numbers of such possibilities. Scaling of the normal modes is sometimes done to standardize their elements associated with amplitudes in various degrees of freedom which is known as normalization. For example for the two storey frame discussed

Let us normalize the vector {*Ïˆ*}(0) by dividing each of the elements by *M*1 and call it as {*Ï†*}(1) denoted as the normalized eigenvector corresponding to first mode. It is given by

In future chapters we will use a normalized vector {*Ï†*} instead of an eigenvector {*Ïˆ*}.

**Influence coefficient method**

The equations of motion of an MDOF system can also be written in terms of influence coefficients which are extensively used in structural engineering. Basically one set of influence coefficients can be associated with each of the matrices involved in the equation of motion. For some problems it is easier to find the flexibility matrix (inverse of stiffness matrix) rather than the stiffness matrix.

Assume three masses are attached to a string as shown in Fig. 10.4. The tension in the string can be assumed to be equal to *P*. We can develop dynamic equation of equilibrium as shown below.

First develop the flexibility matrix. The flexibility influence coefficient *aij *is defined as the deflection at* *â€˜*i*â€™* *due to unit force at* *â€˜*j*â€™. Apply unit force* *at 1, the vertical displacement at 1 is found out by resolving the tension in the direction of 1 as

For a non-trivial solution to exist, the determinant of the above matrix should be equal to zero.

The characteristic equation is

For civil engineering problems all roots are real. If roots are real one can apply Cordonâ€™s solution to find the roots of a characteristic equation if it is cubic.

1 Cordonâ€™s solution

Consider the following cubic equation

*Î»*3 + *b**Î»*2 + *c**Î»* + *d* = 0 ,,,, ,,,, 10.67

*Î±* = 0.407rad = 23.346Â o

*y*1* *= 2* *Ã-* *1.7635 cos (0.407/3) = 3.527

*Y*2* *=* *-2* *Ã-* *1.7635 cos [(0.407 +* **Ï€*)/3] = 2

*Y*2* *= 2* *Ã-* *1.7635 cos [(0.407 +* **Ï€*)/3] = 1.16

*Î»*1 = *y*1 - *b*/3 = 3.527 + 3.33 = 6.85

*Î»*2 = *y*2 - *b*/3 = -1.33 + 3.33 = 2 *Î»*3 = *y*3 - *b*/3 = -2.17 + 3.33 = 1.16

Check *Î»*1 + *Î»*2 + *Î»*3 = 10 = 3 + 3 + 4 = trace of the matrix.

We can get the characteristic equation and the solution of the characteristic equation using MATHEMATICA as shown below.

**Program 10.2: MATLAB program to find the frequencies and normalized mode shapes**

% program to get normalized vectors and eigen values

clc;

close all;

m=[1 0 0;0 1 0;0 0 1];

disp(â€˜ mass matrixâ€™)

m

%you can give stiffness matrix

% disp(â€˜ stiffness matrixâ€™)

% k=[2 -1 0;-1 2 -1;0 -1 1];

% k

% a=inv(k);

% or you can given flexibility matrix directly

a=[.75 .5 .25;.5 1 .5;.25 .5 .75]; disp(â€˜ flexibility matrixâ€™)

a c=a*m;

[ms,ns]=size(m);

%eigen values and eigen vectors [V,D]=eig(c);

for i=1:ms e(i)=1/D(i,i);

end Qh=max(e)+0.001; Ql=0;

for i=1:ms for j=1:ms

if e(j) > Ql & e(j) < Qh kk=j;

Qh=e(j); else

end end Ql=Qh;

Qh=max(e)+0.001;

om1(i)=e(kk);

omega(i)=sqrt(e(kk)); for l=1:ms

p1(l,i)=V(l,kk); end

end

%Normalizing the mode shape L=p1'*m*p1;

%develop modal matrix for i=1:ms

for j=1:ms p(i,j)=p1(i,j)/sqrt(L(j,j));

end end

disp(â€˜ Natural frequencies in rad/secâ€™) disp(omega)

disp(â€˜ normalized modal vectorâ€™) disp(p)

pâ€™*m*p

OUTPUT

mass matrix

m =

1 0 0

0 1 0

0 0 1

a =

0.7500 0.5000 0.2500

0.5000 1.0000 0.5000

0.2500 0.5000 0.7500

Natural frequencies in rad/sec

0.7654 1.4142 1.8478

normalized modal vector

0.5000 -0.7071 -0.5000

0.7071 0.0000 0.7071

0.5000 0.7071 -0.5000

ans = proof of normality principle

1.0000 0.0000 0.0000

0.0000 1.0000 0.0000

0.0000 0.0000 1.0000

Example 10.1

Determine the stiffness matrix and write the dynamic equation of equilibrium for the frame shown in Fig. 10.6. Neglect the effect of axial stiffness of the members AB, BC and CD.

Solution

In total there are six degrees of freedom, three at B and three at C. But we neglect axial deformations in all the three members.

Generalized coordinate = total coordinates - constraint equations = 6 - 3 = 3

The three generalized coordinates are *q*1, *q*2, the rotations at B and C and horizontal sway *q*3 at C. For every member let us consider the moment stress resultants as shown in Fig. 10.7.

The deformations are denoted by *Î´*. We can establish the compatibility matrix which gives relationship between deformations and generalized coordinates (see Fig. 10.8)

{*Î´*} = [*Î²*]{*q*}

The first, second and third columns of [*Î²*] matrix are established by applying unit values for the generalized coordinates as

From the principle of virtual work it can be proved that the generalized forces are written in terms of element stress resultants as

{*Q*} = [*Î²*]*t*{*p*} = [*Î²*]*t*[*k*][*Î²*] = [*K*]{*q*}

so the structural stiffness matrix can be written in terms of element stiffness matrix as

The stiffness matrix could also be obtained straight away by applying the influence coefficient method.

The mass matrix is given by

Hence *dynamic equilibrium* equations can be written as

[ *M* ]{*q***Ë™Ë™**} + [ *K* ]{*q*} = {*Q*}

Since the structure has three degrees of freedom, Cordonâ€™s solution can be applied to evaluate the frequencies and mode shapes.

Example 10.2

Determine the stiffness matrix and dynamic equation of equilibrium for the frame shown in Fig. 10.9. Neglect axial deformations.

Solution

Let us consider total six displacements, three at each node B and C as shown in Fig. 10.9d.

Generalized coordinates = total coordinates - constraint equations

= 6 - 3 = 3

viz. *q*1, *q*2, *q*3.

Consider any member *IJ* as shown in Fig. 10.10. Axial deformation can be written as

âˆ†*IJ* = (*UJ* - *UI*)*C* + (*VJ* - *VI*)*S*

Where

We can write the three constraint equations corresponding to axial deformations of the three members as

âˆ†*AB* = *l*(*U*1 + *U*3) = 0

âˆ†*BC* = *L*(*U*2 - *U*1) + 0(*U*4 - *U*3) = 0; i.e (*U*2 - *U*1) = 0

âˆ†*CD* = *l*(-*U*2 + *U*4) = 0 i.e. (-*U*2 + *U*4) = 0

Writing this in matrix form

Selecting constraint coordinate as *U*1 , the other three displacements can be written as in terms of constrained coordinate as

solving we get *U*2 = *U*1; *U*3 = -*U*1; *U*4 = *U*1.

We can write relation between constrained coordinates and generalized coordinates as

Now we can get the relationship between generalized displacements and deformations (see Fig. 10.11). The relationship between element deformations and generalized displacements is given as (*l* = 0.707). (Assume lengths of all the members are equal to *L* = 1).

Since there are only three degrees of freedom Cordonâ€™s solution can be applied to evaluate the frequencies.

**10.10 ****Program 10.3: MATLAB program for solving structural problem by the stiffness method**

% Analysis of structures by stiffness method - semi direct approach clc;

clear all; r=input(â€˜ENTER Ã'0Ã' FOR DATA FROM FILE,Ã'1Ã' FOR DATA FROM TERMINAL,Ã'2Ã' FOR INTERACTIVEâ€™); ff1=fopen(â€˜stout.datâ€™,â€™wâ€™);

if r<1 ff=fopen(â€˜st.datâ€™,â€™râ€™); x1=fscanf(ff,â€™%sâ€™,1); disp(x1); rb=fscanf(ff,â€™%fâ€™,1); cb=fscanf(ff,â€™%fâ€™,1); for i=1:rb

for j=1:cb c1=fscanf(ff,â€™%fâ€™,1); b(i,j)=c1;

end end disp(b);

x2=fscanf(ff,â€™%sâ€™,1); disp(x2); rk=fscanf(ff,â€™%fâ€™,1); ck=fscanf(ff,â€™%fâ€™,1); for i=1:rk

for j=1:ck c2=fscanf(ff,â€™%fâ€™,1); k(i,j)=c2;

end end disp(k);

x3=fscanf(ff,â€™%sâ€™,1); disp(x3); rdo=fscanf(ff,â€™%fâ€™,1); cdo=fscanf(ff,â€™%fâ€™,1); for i=1:rdo

for j=1:cdo c3=fscanf(ff,â€™%fâ€™,1); do(i,j)=c3;

end end

disp(do); x4=fscanf(ff,â€™%sâ€™,1); disp(x4); rfo=fscanf(ff,â€™%fâ€™,1); cfo=fscanf(ff,â€™%fâ€™,1);

for i=1:rfo for j=1:cfo

c4=fscanf(ff,â€™%fâ€™,1); fo(i,j)=c4;

end end

disp(fo); x5=fscanf(ff,â€™%sâ€™,1); disp(x5); rpo=fscanf(ff,â€™%fâ€™,1); cpo=fscanf(ff,â€™%fâ€™,1); for i=1:rpo

for j=1:cpo c5=fscanf(ff,â€™%fâ€™,1); po(i,j)=c5;

end end

disp(po); x6=fscanf(ff,â€™%sâ€™,1); disp(x6); rsk=fscanf(ff,â€™%fâ€™,1); csk=fscanf(ff,â€™%fâ€™,1); for i=1:rsk

for j=1:csk

c5=fscanf(ff,â€™%fâ€™,1); ks(i,j)=c5;

end end

disp(ks); else

if r<2

%INPUT FOR BETA MATRIX

b=[-1.414 0 0;-1.414 1 0;2 1 0;2 0 1;-1.414 0 1;-1.414 0 0]; %disp(b);

%INPUT FOR ELEMENT STIFFNESS MATRIX k

k=[4 2 0 0 0 0; 2 4 0 0 0 0;0 0 4 2 0 0;0 0 2 4 0 0;0 0 0 0 4 2;0 0 0 0 2 4]; disp(k);

%INPUT FOR INITIAL STRAIN DO do=[0;0;0;0;0;0];

disp(do);

%INPUT FOR APPLIED FORCE fo=[0;0;0];

disp(fo);

%INPUT FOR INITIAL STRESS pO po=[0;0;0;0;0;0];

disp(po);

%input spring stiffness ks=[0,0,0;0,0,0;0,0,0]; disp(ks);

else

rb=input(â€˜NUMBER OF ROWS IN B MATRIXâ€™); cb=input(â€˜NUMBER OF COLUMNS IN B MATRIXâ€™); for i = 1:rb

for j = 1:cb

bm = input(â€˜ â€˜); b(i,j)=bm;

end end disp(b);

rk=input(â€˜NUMBER OF ROWS IN ELEMENT STIFFNESS MATRIXâ€™); ck=input(â€˜NUMBER OF COLUMNS IN ELEMENT STIFFNESS MATRIXâ€™); for i=1:rk

for j=1:ck c2=input(â€˜â€™); k(i,j)=c2; end

end

disp(k);

rdo=input(â€˜NUMBER OF ROWS IN INITIAL STRAIN MATRIXâ€™); cdo=input(â€˜NUMBER OF COLUMNS IN INITIAL STRAIN MATRIXâ€™); for i=1:rdo

for j=1:cdo c3=input(â€˜â€™); do(i,j)=c3; end

end disp(do);

rfo=input(â€˜NUMBER OF ROWS IN APPLIED FORCESâ€™); cfo=input(â€˜NUMBER OF COLUMNS IN APPLIED FORCESâ€™); for i=1:rfo

for j=1:cfo c4=input(â€˜â€™); fo(i,j)=c4; end

end disp(fo);

rpo=input(â€˜NUMBER OF ROWS IN INITIAL STRESS MATRIXâ€™); cpo=input(â€˜NUMBER OF COLUMNS IN INITIAL STRESS MATRIXâ€™); for i=1:rpo

for j=1:cpo c5=input(â€˜â€™); po(i,j)=c5; end

end disp(po); end

end K1=(bâ€™*k*b); K=K1+ks;

ffi=fo-(bâ€™*po)+(bâ€™*k*do); u=inv(K)*ffi;

p=(k*b*u); pf=po+p-(k*do);

q=input(â€˜ENTER Ã'0Ã' FOR O/P IN FILE,Ã'1" for O/P IN TERMINALâ€™) ; if q>0

â€˜DISPLACEMENTS AT GENERALIZED COORDINATESâ€™ disp(u);

â€˜ELEMENT FORCES DUE TO APPLIED FORCESâ€™ disp(p);

â€˜FINAL ELEMENT FORCESâ€™

disp(pf); else

fprintf(ff1,â€™\r\nDISPLACEMENTS AT GENERALIZED COORDINATESâ€™); fprintf(ff1,â€™\r\n%10.4f â€˜, u);

fprintf(ff1,â€™\r\nELEMENT FORCES DUE TO APPLIED FORCESâ€™); fprintf(ff1,â€™\r\n%10.4fâ€™, p);

fprintf(ff1,â€™\r\nFINAL ELEMENT FORCESâ€™); fprintf(ff1,â€™\r\n%10.4fâ€™,pf);

end

Output

>> K

K =

95.9855 3.5160 3.5160

3.5160 8.0000 2.0000

3.5160 2.0000 8.0000

**Static condensation of stiffness matrix**

In most of the structural dynamics problems, mass moment of inertia is only included corresponding to translational degrees of freedom. However, a more generalized formulation of stiffness matrix includes terms corresponding to rotational degrees of freedom. It is necessary to eliminate the extraneous degree of freedom associated with rotation from the stiffness matrix. This method is called *static condensation.* This enables the number of degrees of freedom to be reduced to a desired level.

To establish the equations used in static condensation, assume that [*K*] and the corresponding force vector are partitioned as

The stiffness is now suitable for use with the lumped mass matrix to write equations of motion for translational vibration. Once translational degrees of freedom are obtained, rotational degrees of freedom can be obtained from Eq. 10.64.

**General viscous damping**

The differential equation governing the free vibration of an MDOF system with general viscous damping is given by

If the damping matrix is arbitrary, it is somewhat complex to solve the equations. The above equation is reformulated as 2*n* first order differential equations by writing

The values of *Î³* occur in complex conjugate pairs. The system is stable only if all eigenvalues have non-negative real parts. Eigenvectors correspond to complex conjugates of one another. Eigenvectors correspond to eigenvalues which are not complex conjugates satisfy the orthogonality relationship

Example 10.3

Plot the free vibration response of the shear frame shown in Fig. 10.12 under initial conditions.

All the three displacements are plotted as shown in Fig. 10.13.

**Program 10.4: MATLAB program for free vibration of MDOF with general damping**

% FREE VIBRATION OF GENERAL DAMPED THREE DEGREE OF FREEDOM SYSTEM

clc; clear all; digits(5)

% mass matrix disp(â€˜ mass matrixâ€™)

m=[40000 0 0;0 20000 0;0 0 20000]; m

[ms,ns]=size(m); %damping matrix disp(â€˜ damping matrixâ€™)

c=[0 0 0;0 -8000 8000;0 -8000 8000]; c

%stiffness matrix disp(â€˜ stiffness matrixâ€™)

k=[30000 -10000 0;-10000 20000 -10000;0 -10000 10000]; k

z=zeros(ms); MT=[z m;m c]; KT=[-m z;z k]; Z=inv(MT)*(KT); [V,D]=eig(Z); disp(â€˜ eigenvaluesâ€™) for i=1:2*ms

DS(i)=D(i,i); end

DS

disp(â€˜ eigenvectorsâ€™) V

disp(â€˜ initial conditionsâ€™) x0=[0;0;0;0;0.01;0.05];

disp(â€˜ constants of integrationâ€™) S=inv(V)*x0 tk=linspace(0,40,2001);

%Evaluation of time dependent response %recall that x1=y4; and x2=y5; x3=y6 for k=1:2001

t=tk(k);

for i=ms+1:2*ms x(k,i-ms)=0; for j=1:2*ms

x(k,i-ms)=x(k,i-ms)+(real(S(j))*real(V(i,j))-...

imag(S(j))*imag(V(i,j)))*cos(imag(D(j,j))*t); x(k,i-ms)=x(k,i-ms)+(imag(S(j))*real(V(i,j))-...

real(S(j))*imag(V(i,j)))*sin(imag(V(j,j))*t); x(k,i-ms)=x(k,i-ms)*exp(-real(D(j,j))*t);

end end

end figure(1)

plot(tk,x(:,1),â€™kÃ'â€™,tk,x(:,2),â€™kâ€™,tk,x(:,3),â€™k-.â€™) xlabel(â€˜ t(sec)â€™)

ylabel(â€˜ displacement in mâ€™) gtext(â€˜Ã'Ã'Ã'Ã' x1â€™) gtext(â€˜______ x2â€™) gtext(â€˜-.-.-.-.-.x3â€™)

Output mass matrix

ans =

-0.0882 - 1.2170i -0.0882 + 1.2170i 0.0367 - 0.8097i 0.0367 + 0.8097i 0.0516 - 0.3537i 0.0516 + 0.3537i

eigenvectors

V =

0.2012 - 0.0597i 0.2012 + 0.0597i -0.0206 - 0.4559i -0.0206 + 0.4559i -0.0156 - 0.0762i -0.0156 + 0.0762i

-0.6333 -0.6333 0.1004 - 0.1794i 0.1004 + 0.1794i -0.0280 - 0.1935i - 0.0280 + 0.1935i

0.0741 - 0.3841i 0.0741 + 0.3841i 0.3371 + 0.1799i 0.3371 - 0.1799i - 0.0379 - 0.2600i -0.0379 + 0.2600i

0.0608 + 0.1609i 0.0608 - 0.1609i 0.5630 0.5630 0.2172 - 0.0124i 0.2172 + 0.0124i

-0.0375 - 0.5177i -0.0375 + 0.5177i 0.2155 + 0.1337i 0.2155 - 0.1337i 0.5470 + 0.0005i 0.5470 - 0.0005i

0.3184 + 0.0378i 0.3184 - 0.0378i -0.2405 + 0.4054i -0.2405 - 0.4054i 0.7349 0.7349

initial conditions x0 =

0

0

0

0

0.0100

0.0500

constants of integration

S =

0.0049 - 0.0109i

0.0049 + 0.0109i -0.0135 - 0.0058i -0.0135 + 0.0058i 0.0237 + 0.0317i 0.0237 - 0.0317i

**Newmarkâ€™s numerical integration**

In 1959, N M Newmark developed a family of time stepping methods. The method discussed for the SDOF system in Chapter 7 may be extended to the MDOF system.

*Algorithm*

â€¢ Average acceleration method *Î³* = 1/2; *Î²* = 1/4

Linear acceleration method *Î³* = 1/2; *Î²* = 1/6

3. Repeat step 2 for next time step by replacing *i* by *i* + 1

**Program 10.5: MATLAB program for Newmarkâ€™s method of MDOF with generalized damping**

%free vibration of mdof system damped using newmarkâ€™s method

%mass matrix

clc; close all;

m=[40000 0 0;0 20000 0;0 0 20000]; disp(â€˜ mass matrixâ€™)

m

%damping matrix

c=[0 0 0;0 8000 -8000;0 -8000 8000]; disp(â€˜ damping matrixâ€™)

c

%stiffness matrix

k=[30000 -10000 0;-10000 20000 -10000;0 -10000 10000]; k

%specify integration parameters for linear acceleration method

%beta=1/6;

%gamma=0.5;

%specify integration parameters for constant acceleration method beta=1/4;

gamma=0.5;

%specify increment in time dt=0.002;

%calculate the following constants b1=1/(beta*dt^2);

%b2=1/(beta*dt);

%b3=beta-0.5;

b4=gamma*dt*b1;

%b5=1+gamma*dt*b2;

%b6=dt*(1+gamma*b3-gamma); kt=k+b1*m+b4*c; ad=(1/(beta*dt))*m+(gamma/beta)*c;

bd=(1/(2*beta))*m+dt*(gamma/(2.0*beta)-1)*c; %specify initial displacements

u0=[0 0.01 0.05]â€™; v0=[0 0 0]â€™;

a0=inv(m)*(-c*v0-k*u0); t(1)=0;

for i=1:3 u(i,1)=u0(i); v(i,1)=v0(i); a(i,1)=a0(i); end

for i=2:20000 t(i)=(i-1)*dt; for j=1:3

uu(j)=u(j,i-1); vv(j)=v(j,i-1);

aa(j)=a(j,i-1); end

fba=(ad*vvâ€™+bd*aaâ€™); kti=inv(kt); up=kti*fba;

vp=(gamma/(beta*dt))*up-(gamma/beta)*vvâ€™+(1-gamma/(2.0*beta))*aaâ€™; ap=b1*up-(1/(beta*dt))*vvâ€™-(1/(2.0*beta))*aaâ€™;

for j=1:3 u(j,i)=u(j,i-1)+up(j); v(j,i)=v(j,i-1)+vp(j); a(j,i)=a(j,i-1)+ap(j);

end end figure(1)

plot(tâ€™,u(1,:)â€™,â€™-â€™,tâ€™,u(2,:)â€™,â€™kâ€™,tâ€™,u(3,:)â€™,â€™-.â€™) gtext(â€˜- - x1â€™)

gtext(â€˜___x2â€™) gtext(â€˜-.-. x3â€™)

When Example 10.3 is solved by Newmarkâ€™s method we get the response as shown in Fig. 10.14.

**Forced response of a three-degrees-of-freedom under-damped system**

Consider forced vibration of a three-degrees-of-freedom under-damped system as shown in Fig. 10.1b. The equation of motion is written in matrix form as

Combining Eq. 10.76 and 10.77 we get six first order linear differential

equations and they are solved using the Runge-Kutta method as explained in Previous Pages.

As an example consider

A MATLAB listing of the program is given below and the displacement response is shown in Fig. 10.15.

Program 10.6 MATLAB program for forced vibration of three degrees of freedom by Rungeâ€'Kutta method

% three degrees of freedom forced vibration by rk method clc;

tspan=[0:.01:10];

y0=[0;0;0;0;0;0]; [t,y]=ode45(â€˜dfunc3fâ€™,tspan,y0); subplot(311)

plot(t,y(:,1)); xlabel(â€˜tâ€™); ylabel(â€˜x1(t)â€™); title(â€˜ x1(t) vs tâ€™);

subplot(312)

plot(t,y(:,3)); xlabel(â€˜tâ€™); ylabel(â€˜x2(t)â€™); title(â€˜x2(t) vs tâ€™); subplot(313) plot(t,y(:,5)); xlabel(â€˜tâ€™); ylabel(â€˜x3(t)â€™); title(â€˜x3(t) vs tâ€™)

function f=dfunc3f(t,y)

% four first order equations are given by f f=zeros(6,1);

%mass matrix m=[100,0,0;0,10,0;0,0,10]; %damping matrix

c=[400,-200,0;-200,400,-200;0,-200,200]; %stiffness matrix k=1000*[8,-4,0;-4,8,-4;0,-4,4];

%give forces force=50*[1;1;1]; om=5.0;

%four first order equations f(1)=y(2);

f(2)=force(1)*cos(om*t)-c(1,1)*y(2)/m(1,1)-c(1,2)*y(4)/m(1,1)...

-c(1,3)*y(6)/m(1,1)-k(1,1)*y(1)/m(1,1)-k(1,2)*y(3)/m(1,1)-k(1,3)*y(5)/ m(1,1);

f(3)=y(4); f(4)=force(2)*cos(om*t)-c(2,1)*y(2)/m(2,2)-c(2,2)*y(4)/m(2,2)...

-c(2,3)*y(6)/m(2,2)-k(2,1)*y(1)/m(2,2)-k(2,2)*y(3)/m(2,2)-k(2,3)*y(5)/m(2,2); f(5)=y(6); f(6)=force(3)*cos(om*t)-c(3,1)*y(2)/m(3,3)-c(3,2)*y(4)/m(3,3)...

-c(3,3)*y(6)/m(3,3)-k(3,1)*y(1)/m(3,3)-k(3,2)*y(3)/m(3,3)-k(3,3)*y(5)/m(3,3)

**Summary**

The vibration of an *n*-degree-of-freedom system is governed by a system of *n *differential equations. The general solution of these differential equations* *is the sum of a homogeneous solution and a particular solution. The homogeneous solution is the free vibration response, and its determination is often necessary before the forced response can be determined. The free vibration analysis of an MDOF system is significantly more complicated than the free vibration analysis of a one or two-degrees-of-freedom system.

Study Material, Lecturing Notes, Assignment, Reference, Wiki description explanation, brief detail

Civil : Structural dynamics of earthquake engineering : Free vibration of multiple degrees of freedom |

**Related Topics **

Privacy Policy, Terms and Conditions, DMCA Policy and Compliant

Copyright Â© 2018-2023 BrainKart.com; All Rights Reserved. Developed by Therithal info, Chennai.