Time history response by mode superposition in relation to structural dynamics during earthquakes Abstract: In this page, the mode superposition method using the mode displacement and the mode acceleration methods is used to find the response of the structure with classical damping. In addition, numerical methods discussed in previous pages for a single-degree-of-freedom system are extended to find the dynamic response for a multiple-degrees-of-freedom system. Relevant programs in MATLAB are also given.
Time history response by mode superposition in relation to structural dynamics during earthquakes
Abstract: In this page, the mode superposition method using the mode displacement and the mode acceleration methods is used to find the response of the structure with classical damping. In addition, numerical methods discussed in previous pages for a single-degree-of-freedom system are extended to find the dynamic response for a multiple-degrees-of-freedom system. Relevant programs in MATLAB are also given.
Key words: normal coordinate, classical damping, time history, response spectrum, SRSS rule, CQC rule.
Introduction
In the multi-storey building or in any system having n degrees of freedom, if principalmodes (normal modes) are used as generalized coordinates, the n dynamic equilibrium equations will be uncoupled. Hence each uncoupled equation can be solved independently as each equation contains one degree of freedom only. One can apply any one of the numerical methods to determine the response of a single-degree-of-freedom system. The response of a multiple-degrees-of-freedom (MDOF) system is obtained by mode superposition by summing the response of the individual modes. This procedure of dynamic analysis is referred to as the mode superposition method, normal mode method or simply modal analysis.
Limitations
Modal analysis is valid for linear systems only. Damping in the system must be proportional to mass and stiffness and this damping is known as classical damping. There are two forms of mode superposition method:
• mode displacement method;
• mode acceleration method.
The mode displacement method is based on stiffness formulation and the mode acceleration method is based on the flexibility approach. Hence the second method is not widely used. We will discuss mode displacement method first followed by the mode acceleration method.
It is observed in structural or mechanical systems that in most types of dynamic loadings, the contributions of the various modes to the dynamic response are the greatest for lowest frequencies and tend to decrease for higher frequencies. Hence for practical applications, it is usually not necessary to include any of the higher modes of vibration in the superposition process.
Mode displacement method for uncoupled system
Method 1
The dynamic equation of motion for an n degree of freedom system may be written as
[m]{u} + [k ]{u} = {F(t )} = {F} f (t ) ---- 12.1
where k and m are the symmetric stiffness and mass matrices respectively, u , u˙˙ are displacement and acceleration in physical coordinates respectively, F(t) is the external force vector. Equation 12.1 is a coupled equation since it cannot be solved independently, only simultaneously. The coupling may be due to stiffness coupling or inertial coupling. The mass matrix will be uncoupled for a lumped mass system and coupled for a consistent mass system. The mode superposition method cannot be applied to the coupled form of equation such as Eq. 12.1.
It is therefore necessary to find a coordinate system which will exhibit neither stiffness coupling nor mass coupling which is the main essence of mode superposition method. The coordinates that enable the decoupling of the equations of motion are called ‘principal coordinates’ or ‘normal coordinates’.
To uncouple the equations we introduce an alternative set of coordinates y such that
y = y(u1, u2,…, un) --- 12.2
i.e. the actual generalized coordinate u can be written in terms of new coordinates y as
{u} = [ψ]{y} - - - - 12.3
where [ψ] is the modal matrix determined from the solution of eigenvalue problem. Each column vector of denotes the normalized eigenvector corresponding to that particular mode. Hence u can be written as
{u} = [ψ]{y} --- 12.4
Hence transforming the equations of motion Eq. 12.1 from physical coordinates to normal coordinates as
P(t) is called modal force vector or modal participation factor.
Since the modal mass = I and modal stiffness matrices = [ ω n2 ] are diagonal can be written in a uncoupled equations as
The complete response for the rth mode can be expressed as the sum of response due to initial conditions and the modal response due to Pr(t).
The solution of Eq. 12.9a is written as
yr = ycr + ypr --- 12.10
where
Thus complete response for the rth mode in normal coordinate can be represented by the Duhamel integral expression as
The procedure discussed in previous pages for evaluating an integral is also applicable to Eq. 12.17. However for random dynamic excitation, it is generally necessary to employ one of the numerical techniques discussed in earlier chapters to obtain the time history response.
Now the exact response in physical coordinate is obtained by summing up all the individual modal responses in normal coordinates given by
Modal analysis is effective when only a few modes of the system are required to render an accurate solution. In the case of tall buildings having hundreds of degrees of freedom, the entire building possesses 100 eigen pairs (i.e. 100 eigenvalues and 100 eigenvectors) that describe the normal vibration mode. If it is known that frequency content of excitation force is in the vicinity of the lowest few frequencies of the building, the higher modes will not be excited and the force response can be determined by superposition of only these few low frequency modes. Hence the displacement u(t) is given as
where p ≤ N. Suppose the excited force vector is given as
Modal participation factor
The forced response of an MDOF system is expressed in Eq. 12.21b as
�'r is called modal participation factor and is particular useful when used in conjunction with the response spectrum analysis of MDOF systems.
Method 2 due to Chopra
Dynamic equation of motion of an n-degrees-of-freedom system is given by
[ m ]{u} + [ k ]{u} = { F} f ( t ) --- 12.24
Assume the force vector is written in terms of normalized eigenvector as
{F} = �' �'ι[m]{φ}I --- 12.25
Premultiplying both sides with {φ}nT we get
�'n = {φ} nT { F} --- 12.26
where �'n is known as the modal participation factor corresponding to mode n. Similarly all the modal participation factors can be found.
Now the force contribution for each mode can be found out as
[F] = [m][φ][�'] -- - 12.27
where [�'] is a diagonal matrix consisting of all modal participation factors. [F] contains n �' columns = number of degrees of freedom of the system. The first column of [F] may be viewed as an expression of the distribution of {F} of applied force in terms of force distribution {F1} associated with natural period ωn1. The force vector {Fnf (t)} produces a response only in the nth mode and no response in other modes. The dynamic response in thenth mode is entirely due to partial force vector {Fnf (t)}
Mode superposition solution for systems with classical damping
Consider a multi-storey frame shown in Fig. 12.8. Usually damping which is inherent property of the structure is present. In addition to stiffness, and masses the damping of each floor is given by c1, c2, c3 respectively.
The equation of motion for a general as shown in viscously damped MDOF Fig. 12.8 is expressed as
[m ]{u˙˙} + [c ]{u˙} + [k ]{u} = { f (t )} --- 12.28
where M is = mass matrix
c is = damping matrix k is = stiffness matrix
Writing generalized coordinates in terms of ‘normal’ coordinates as
{u} = [φ]{y} --- 12.29
Substituting in terms of normalized coordinates Eq. 12.28 becomes
[m ][φ ]{y} + [c ][φ ]{y} + [k ][φ ]{y} = {F (t )} ---- 12.30
Premultiplying with [φ] we get
where M, C and K are the generalized mass, generalized damping and generalized stiffness matrices. If [φ] is the normalized eigenvector
[φ]T[m][φ] = [I]
[φ]T[K][φ] = [ω2] ---- 12.33
C is an n �- n symmetric matrix that is diagonal only for a special case of [c]. This special case of [c] is referred to a classical damping or proportional damping for which [c] is proportional to m and k as
where [C] is a diagonal matrix and is called the modal damping matrix, in which,
cr = α + β ω r2 ---- 12.35
Cr is called the modal damping coefficient for rth mode. Assume modal damping factor is given by
The solution of Eq. 12.38 for the rth mode can be expressed by the same form as a single-degree-of-freedom (SDOF) system and the solution is given by
Frequently, in the absence of more definite information about damping, a reasonable value for modal damping as Eq. 12.38 is simply assumed to be valid.
If the solution of ‘n’ modal equations given by Eq. 12.38 are substituted in Eq. 12.29
Then the exact system response u(t) in physical coordinates for all n mode is obtained. However, if only p modes are retained in the solution then
rather than Eq. 12.41 is used to define the truncated response. The resulting mode superposition solution omits the contribution of the modes (p + 1) to n.
The damping discussed above is called Rayleigh damping, which is just one example of proportional damping. This damping is also referred to as classical damping, orthogonal damping and modal damping.
In Rayleigh damping the modal damping matrix [c] is defined as
Rayleigh damping is therefore defined for an MDOF system by specifying for two different unequal frequencies of vibration and solving by simultaneous solution of the equation
Solving Eqs 12.45a and 12.45b we get α and β and hence [c] matrix is defined. It is also seen in Eq. 12.44 that contribution of mass in damping is
inversely proportional to ωr and contribution of stiffness is proportional to
ωr.
Therefore for large ωr the stiffness proportional term determines the system damping. This trend generally leads to unrealistically high damping ratios in higher modes. For a large MDOF system fortunately the modal superposition methods, we consider p < < n modes only and hence damping in the higher modes is generally not a critical issue. In modal analysis it is not necessary to formulate or explicitly give the equation for [C] displacement since the values of ζr are required. But [C] is required for dynamic response of MDOF system by direct numerical integration procedure.
Numerical evaluation of modal response
The numerical integration procedures discussed in earlier chapters for an SDOF system can readily be applied to MDOF systems. Using modal analysis we get n independent (uncoupled) differential equations which are then evaluated at discrete time intervals by any numerical procedure. These normal coordinates are transformed back to physical coordinatesu and updated after each time step.
One has to resort to numerical integration procedure since the forcing function is arbitrary and cannot be expressed by means of a simple analytical expression. The algorithm is given in Table 12.1. As has been already seen for accuracy t ≤ Tn/10 should be selected where Tn is the natural period corresponding to highest ratio mode associated with the system. When fewer modes are considered ∆t ≤ Tp/10 where Tp is the period of pth mode.
Program 12.1: MATLAB program for dynamic response using modal superposition
% program to get dynamic response of MDOF using modal superposition %force.dat contains data time, force1, force2, force3
clc; close all;
m=[18348 0 0;0 13761 0;0 0 9174]; [nd,nd]=size(m);
disp(‘ mass matrix’) m
%if forces are acting at degrees of freedom force=‘force.dat’
f=load(force);
%if base ground acceleration is given
%dis=‘disp.dat’
%di=load(dis);
%% convert to equivalent nodal loads
%for i=1:nd
%f(:,i)=-di*m(i,i);
%end
%you can give stiffness matrix disp(‘ stiffness matrix’)
k=[.74504e8 -.64228e8 0;-.64228e8 .14746e9 -.8324e8;0 -.8324e8 .8324e8]; 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
cc=a*m;
[ms,ns]=size(m);
%eigen values and eigen vectors [V,D]=eig(cc);
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 ph(i,j)=p1(i,j)/sqrt(L(j,j));
end end
disp(‘ Natural frequencies in rad/sec’) disp(omega)
disp(‘ normalized modal vector’) disp(ph)
ph’*m*ph
%give alpha and beta alpha=0.6978; beta=9.4e-4;
disp(‘ damping matrix’) cd=alpha*m+beta*k; cd
for i=1:ms zeta(i)=0.5*(alpha/omega(i)+beta*omega(i));
end
%give initial displacements and velocities of all degrees of freedom xin=[0; 0.0; 0.0];
vin=[.0; 0.0; 0];
disp(‘ initial displacements of all degrees of freedom’) xin
disp(‘ initial velocities of all degrees of freedom’) vin
xmin=ph’*m*xin; vmin=ph’*m*vin;
%define forces at various times at all degrees of freedom dt=0.02;
ttot=5.0;
p=f*ph; for ii=1:nd
x(1,ii)=xmin(ii,1);
x1(1,ii)=vmin(ii,1);
c=2.0*zeta(ii)*omega(ii); x2(1,ii)=p(1,ii)-c*x1(1,ii)-omega(ii)^2*x(1,ii);
ks=omega(ii)^2+2.0*c/dt+4.0/(dt^2);
ic=1;
for t=0:dt:ttot-dt dps=p(ic+1,ii)+(4.0*x(ic,ii)/dt^2+4.0*x1(ic,ii)/dt...
+x2(ic,ii)+c*2.0*x(ic,ii)/dt+x1(ic,ii));
x(ic+1,ii)=dps/ks; x2(ic+1,ii)=(4.0/dt^2)*(x(ic+1,ii)-x(ic,ii)-...
dt*x1(ic,ii))-x2(ic,ii); x1(ic+1,ii)=(2.0/dt)*(x(ic+1,ii)-x(ic,ii))-x1(ic,ii); ic=ic+1;
end end u=x*ph’; v=x1*ph’;
acn=x2*ph’; tt=linspace(0,5,251); figure(1); plot(tt,u(:,3),‘k’) xlabel(‘ time in secs’);
ylabel(‘ roof displacement in m’);
title(‘ displacement response of the roof’) figure(2);
plot(tt,v(:,3),‘k’) xlabel(‘ time in secs’);
ylabel(‘ roof velocity in m/sec’); title(‘ velocity response of the roof’) figure(3);
plot(tt,acn(:,3),‘k’) xlabel(‘ time in secs’);
ylabel(‘ roof acceleration in m/sq.sec’); title(‘ acceleration response of the roof’)
OUTPUT OF MATLAB
mass matrix
m =
18348 0 0
0 13761 0
0 0 9174
force = force.dat
stiffness matrix
k =
74504000 �'64228000 0
�'64228000 147460000 �'83240000
0 �'83240000 83240000
flexibility matrix
a =
1.0e-006 *
0.0974 0.0974 0.0974
0.0974 0.1130 0.1130
0.0974 0.1130 0.1250
damping matrix
cd =
1.0e+005 *
0.8284 �'0.6037 0
�'0.6037 1.4821 �'0.7825
0 �'0.7825 0.8465
Natural frequencies in rad/sec
15.3317 74.7946 134.2407
normalized modal vector
�'0.0046 �'0.0055 0.0016
�'0.0051 0.0024 �'0.0064
�'0.0052 0.0063 0.0065
ans =
1.0000 0.0000 �'0.0000
0.0000 1.0000 0.0000
�'0.0000 0.0000 1.0000
initial displacements of all degrees of freedom
xin =
0
0
0
initial velocities of all degrees of freedom
vin =
0
0
0
Dynamic analysis using direct integration methods
In Chapter 7, we saw the dynamic analysis of a SDOF system using the direct integration method. Here Newmark’s constant and linear acceleration methods can be extended to get dynamic response for MDOF systems. Table 12.2 gives the algorithm for the dynamic analysis of MDOF systems using Newmark’s method. This algorithm can be modified so as to apply Wilson-θ or any other numerical method discussed in previous pges.
The program developed in MATLAB is applied to Example 12.4 and the dynamic displacement, velocity and acceleration response of the roof are obtained as shown in Fig. 12.12.
Program 12.2: MATLAB program for finding dynamic response of MDOF using direct integration method (Newmark’s method)
%dynamic analysis using direct integration method
%input mass matrix
*force.dat contains data for forces as time, force1, force2,force3 clc;
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 force=‘force.dat’
f=load(force);
%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
%input damping matrix
c=100000*[0.82836 -0.6037 0;-0.6037 1.4821 -0.7824;0 -0.7824 0.84647]; disp(‘ damping matrix’)
%input stiffness matrix
k=[74504937 -64228394 0;-64228394 147468394 -83240000;0 -83240000 83240000];
disp(‘ stiffness matrix’) k
format long; kim=inv(k)*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 0.0 0.0];
disp(‘ initial displacements’) u0
%give initial condition for velocities v0=[0. 0. 0.];
disp(‘ initial velocities’) v0
for i=1:ns a0=inv(m)*(f(1,:)’-c*v0'-k*u0');
end dt=0.02;
kba=k+(gama/(beta*dt))*c+(1/(beta*dt*dt))*m;
kin=inv(kba);
aa=(1/(beta*dt))*m+(gama/beta)*c; bb=(1/(2.0*beta))*m+dt*(gama/(2.0*beta)-1)*c; u(1,:)=u0;
v(1,:)=v0;
a(1,:)=a0; for i=2:251
df(i,:)=f(i,:)-f(i-1,:)+v(i-1,:)*aa’+a(i-1,:)*bb’; du(i,:)=df(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 tt=linspace(0,5,251); figure(1); plot(tt,u(:,3),’k’); xlabel(‘ time in secs’);
ylabel(‘ roof displacement’);
title(‘ displacement response of the roof’); figure(2);
plot(tt,v(:,3),’k’); xlabel(‘ time in secs’); ylabel(‘ roof velocity’);
title(‘velocity response of the roof’); figure(3);
plot(tt,a(:,3),’k’); xlabel(‘ time in secs’);
ylabel(‘ roof acceleration’);
title(‘ acceleration response of the roof’)
force =
force.dat
damping matrix
c =
82836 �'60370 0
�'60370 148210 �'78240
0 �'78240 84647
stiffness matrix
k =
74504937 �'64228394 0
�'64228394 147468394 �'83240000
0 �'83240000 83240000
natural frequencies
omega =
1.0e+002 *
0.15323852988388 0.7469278455770 1.34226456460257
initial displacements
u0 =
0 0 0
initial velocities
v0 =
0 0 0
Normal mode response to support motions
In many structures there is a great interest in studying the response of MDOF systems due to support motions such as the earthquake response of structures (see Fig. 12.13).
Let ug be the ground displacement due to earthquake motion. The system equilibrium equations is physical components is written as
The procedure is similar to the force acting at the coordinates.
Example 12.5
The base of the three storeyed shear frame building shown in Fig.12.9 is excited by horizontal ground acceleration u˙˙g shown in Fig. 12.14. Determine the dynamic response of the structure by mode superposition method. Plot the histories for relative displacement, relative velocity and total absolute acceleration for the top level of the building for the time interval 0 ≤ t ≤ 5.
Solution
The displacement, velocity and total acceleration history at the roof are shown in Fig. 12.15.
Response spectrum analysis
In Previous pages we discussed the response spectra or shock spectra for a specific dynamic disturbance particular for an SDOF system. The technique is also applicable to MDOF systems. Usually the design engineer is interested in finding the maximum response of a system to a specified input. In such cases response spectrum analysis is useful. The maximum modal response of rth mode of an MDOF system to a specified input is expressed in physical coordinates as
An upper bound for the total response may be obtained by summing numerically the maximum response for each mode. Hence maximum response can be expressed as
The numerical sum as given in Eq. 12.49 is a very conservative estimate. This shows that maximum response for all modes occur simultaneously.
This is neither reasonable nor correct. Many statistical methods for combining the maximum modal response have been developed and several of these modal combinations will be discussed later. A popular method is square root of sum of squares (SRSS) given by
The SRSS method for combining modal maximum has been shown to render acceleration approximation to two-dimensional structural systems exhibiting well-separated vibration frequencies. For three-dimensional systems and other systems with closely spaced modes, the complete quadratic combination (CQC) method renders significant improvement in estimating the response. (Two consecutive modes are assumed as closely spaced of their corresponding frequencies differ from each other by 10% or less of the lower frequency.) The CQC combination is expressed as
where Rr and Rs are the peak values or the particular response for the rth and sth modes respectively. For constant damping ζ = ζ1 = ζ2 = … = ζp and Prs is given is
Mode acceleration method
An alternative form of mode superposition method is mode acceleration method. Since it is a flexibility approach, it is not usually employed in dynamics. Compared with mode superposition analysis it exhibits improved convergence characteristics. In this method we require fewer vibration modes.
The equation of motion can be written as
where [K]�'1 = [a]. If the acceleration vector is approximately uˆ then truncated mode acceleration solution is given by
But u˙˙(t ) (physical coordinates) can be written in terms of normalized coordinates as
The first term of Eq. 12.60 represents the pseudo-static response and the second term is the mode acceleration, which gives the method its name. Since ω is in the denominator, this method provides quadratic convergence. Even though it is enough to consider ‘p’ modes, it is also important to consider all the roots whose frequencies are close proximity to any excitation frequency.
Summary
In this page, time history analysis has been carried out by modal superposition and numerical integration methods. In addition, the modal acceleration method is also applied to find the dynamic response of the structures.