Abstract: There are two basic approaches to numerically evaluate the dynamic response. The first approach is numerical interpolation of the excitation and the second is numerical integration of the equation of motion. Both approaches are applicable to linear systems but the second approach is related to non-linear systems. In this page and previous page, various numerical techniques based on interpolation, finite difference equation and assumed acceleration are employed to arrive at the dynamic response due to force and base excitation.
Response to base
excitation
Evaluating the dynamic response of structures due to arbitrary
base or support motion can generally be facilitated by numerical
integration methods. The response of structures to base excitation is
the most important consideration in earthquake engineering. Consider the one
storey shear frame shown in Fig. 7.13. The structure experiences any arbitrary
ground displacements or acceleration u˙˙_{g} ( t
). It is assumed that the shear frame is attached to a rigid base that moves
with the ground. In analysing the structural response there are several
components of motion that must be considered; specifically u_{i}.
The relative displacement of the structure and u_{T}
the total or absolute displacement of the structure measured from reference
axis.
The zero on the right hand side
of Eq. 7.61 would suggest that the structure is not subjected to any external
load F(t). This is not entirely true since the ground motion
creates the inertia of forces in the structure.
Thus, noting total displacement of the mass u_{t}
is given by
The term mu˙˙_{g} can be thought
of as an effective load F_{eff} (t) applied externally in
the mass of the structure shown in Fig. 7.14. Therefore the equation of motion
Note that in Eq. 7.69. u _{i}
, u˙ _{i} , u˙˙_{i}
represent relative displacement, velocity and acceleration of the mass
respectively. We will see in later chapters that we can convert an ‘n’-degrees-of-freedom
system to ‘n-single-degree-of-freedom system’ for linear vibration with
proportional damping.
Equation 7.69 can be integrated
by any of the methods discussed in the chapter. We will apply the Wilson-θ method. This is also useful in
establishing response spectra.
Example 7.5
A single storey shear frame shown in Fig 7.15a is subjected to
El Centro ground excitation as shown in Fig. 7.15c. The simplified model is
shown in
Fig. 7.15b. The rigid girders
support a load of 25.57kN/m. Assuming a damping factor ρ = 0.02 for steel frame, E
= 200GPa. Write a computer program for the Wilson-θ method to evaluate dynamic
response of the frame and plot u(t), v(t) and a_{t}(t)
in the interval.
Solution
ω_{n} =
9.53rad/s
T_{n}
=
0.659s
ρ = 0.02
The absolute acceleration can be obtained. The program using
Wilson-θ method is given below for ground movement.
The displacement, velocity and acceleration (total) response are shown in Fig
7.16.
Program 7.8: MATLAB program for dynamic response to base
excitation using Wilson-θ
method
% Response due
to ground excitation using Wilson-Theta method
%**********************************************************
ma=1;
k=90.829;
wn=sqrt(k/ma);
tn=6.283/wn;
theta=1.42;
r=0.02;
c=2.0*r*sqrt(k*ma);
u(1)=0;
v(1)=0;
tt=50.0;
n=2500;
n1=n+1;
dt=tt/n;
d=xlsread(‘eqdata’);
for i=1:n1;
ug(i)=d(i,2);
p(i)=-ug(i)*9.81;
end;
an(1)=(p(1)-c*v(1)-k*u(1))/ma;
kh=k+3.0*c/(theta*dt)+6.0*ma/(theta*dt)^2;
a=6.0*ma/(theta*dt)+3.0*c;
b=3.0*ma+theta*dt*c/2.0;
for i=1:n1;
s(i)=(i-1)*dt;
end;
for i=2:n1;
ww=(p(i)-p(i-1))*theta+a*v(i-1)+b*an(i-1);
xx=ww/kh;
zz=(6.0*xx/((theta*dt)^2)-6.0*v(i-1)/(theta*dt)-3.0*an(i-1))/theta;
yy=dt*an(i-1)+dt*zz/2.0;
v(i)=v(i-1)+yy;
an(i)=an(i-1)+zz;
vv=dt*v(i-1)+dt*dt*(3.0*an(i-1)+zz)/6.0;
u(i)=u(i-1)+vv;
end;
% Find total
acceleration
for i=1:n1;
an(i)=an(i)+ug(i)*9.81;
end;
figure(1);
plot(s,u);
xlabel(‘
time (t) in seconds’);
ylabel(‘ Response displacement (relative) u in m’); title(‘
dynamic response’);
figure(2);
plot(s,v);
xlabel(‘
time (t) in seconds’);
ylabel(‘ Response velocity (relative) v in m/sec’); title(‘
dynamic response’);
figure(3);
plot(s,an);
xlabel(‘
time (t) in seconds’);
ylabel(‘ Response acceleration (total) a in m/sec’); title(‘
dynamic response’);
figure(4);
plot(s,ug);
xlabel(‘ time (t) in seconds’);
ylabel(‘ ground acceleration / g’); title(‘ Elcentro NS’);
Wilson’s procedure
(recommended)
Damped free vibration due to
initial conditions
The equation of motion is written as
in which initial nodal
displacements u_{0} and velocity u˙_{0}
are specified due to previous loading acting on the structure. Note that the
functions s(t) and c(t) given in Table 7.13 are
solutions to Eq. 7.71a.
There are many different methods
available to solve the typical modal equations. However, according to Wilson,
the use of the exact solutions for a linear load over a small time increment
has been found to be the most economical and accurate method to numerically
solve the equations using a computer program. This method does not have
problems such as stability and does not introduce numerical damping. Since the
most seismic ground motion is defined a linear within 0.005s intervals, the
method is exact for the type of loading for all frequencies. All modal
equations are converted to ‘n’ uncoupled equations.
Using the equation as
where time ‘t’ refers to the start of time step. Now
the exact solution within the time step can be written as
where all functions are as defined in Table 7.14.
Equations 7.75, 7.76a and 7.76b are evaluated at the end of
time increment ∆t and the
following displacement, velocity and acceleration at the end of ith time
step are given by the recurrence relation.
The constants A_{1} to A_{12}
are summarized in Table 7.14 and they need to be calculated only once.
Therefore for each time increment only 12 multiplications and 9 conditions are
required. Hence the computer time required solving for 200 steps per second for
50s duration earthquake approximately 0.01s or 100 modal equations can be
solved in 1s of computer time. Therefore, there is no need to consider other
numerical methods (as per Wilson) such as approximate fast Fourier transform,
or numerically evolution of the Duhamel integral to solve the equation.
Because of the speed of this
exact piecewise linear technique, it can also be used to develop accurate
earthquake response spectra using a very small amount of computer time.
Example 7.6
Solve Example 7.4 by Wilson’s
proposed procedure with different data as shown below.
m = 0.065 k = 7.738 ρ = 0.07 F_{0} =
10
Solution
The equation of motion is written as
A program for Wilson’s method is written in MATLAB as shown
below. The force time curves are shown in Fig. 7.18a. The displacement,
velocity and acceleration response are shown in Fig.7.18b, c and d
respectively. Wilson’s method is the best to solve problems involving base excitation.
Program 7.9: MATLAB program
for dynamic response by Wilson’s general method
%Matlab program using Wilson’s
general method %********************************************* ma=0.065;
k=7.738;
wn=sqrt(k/ma)
r=0.07; wd=wn*sqrt(1-r^2); c=2.0*r*sqrt(k*ma); wnb=wn*r;
rb=r/sqrt(1-r^2); tt=3.0;
n=300
n1=n+1
dt=tt/n
td=.75;
jk=td/dt;
a0=2.0*r/(wn*dt);
a1=1+a0 a2=-1/dt;
a3=-rb*a1-a2/wd; a4=-a1;
a5=-a0; a6=-a2;
a7=-rb*a5-a6/wd; a8=-a5;
a9=wd^2-wn^2; a10=2.0*wnb*wd; u(1)=0;
v(1)=0;
for m=1:n1; pa(n)=0.0 p(m)=0.0;
end;
jk1=jk+1; for n=1:jk1;
t=(n-1)*dt;
p(n)=10.0*(1-t/td)*exp(-2.0*t/td);
p(n)=pa(n)/ma;
end; s=exp(-r*wn*dt)*sin(wd*dt);
c=exp(-r*wn*dt)*cos(wd*dt); sd=-wnb*s+wd*c; cd=-wnb*c-wd*s; sdd=-a9*s-a10*c cdd=-a9*c+a10*s;
ca1=c+rb*s;
ca2=s/wd;
ca3=(a1+a2*dt+a3*s+a4*c)/(wn^2);
ca4=(a5+a6*dt+a7*s+a8*c)/(wn^2)
ca5=cd+rb*sd;
ca6=sd/wd;
ca7=(a2+a3*sd+a4*cd)/(wn^2);
ca8=(a6+a7*sd+a8*cd)/(wn^2);
ca9=cdd+rb*sdd;
ca10=sdd/wd;
ca11=(a3*sdd+a4*cdd)/(wn^2);
ca12=(a7*sdd+a8*cdd)/(wn^2);
an(1)=(p(1)-2.0*r*wn*v(1)-(wn^2)*u(1)); for i=1:n1
ti(i)=(i-1)*dt end
for i=2:n1
u(i)=ca1*u(i-1)+ca2*v(i-1)+ca3*p(i-1)+ca4*p(i);
v(i)=ca5*u(i-1)+ca6*v(i-1)+ca7*p(i-1)+ca8*p(i);
an(i)=an(i-1)+ca9*u(i-1)+ca10*v(i-1)+ca11*p(i-1)+ca12*p(i);
end;
figure(1);
plot(ti,u);
xlabel(‘
time (t) in seconds’)
ylabel(‘ Response displacement u in m’) title(‘ dynamic
response’)
figure(2);
plot(ti,v);
xlabel(‘
time (t) in seconds’)
ylabel(‘ Response velocity v in m/sec’) title(‘ dynamic
response’)
figure(3);
plot(ti,an);
xlabel(‘
time (t) in seconds’)
ylabel(‘ Response acceleration a in m/sec’) title(‘ dynamic
response’)
figure(4);
plot(ti,p)
xlabel(‘ time (t) in seconds’)
ylabel(‘ force in Newtons’) title(‘ Excitation Force’)
Response of
elasto-plastic SDOF system
When a steel or reinforced
concrete building is subjected to extreme loading it undergoes elasto-plastic
behaviour. Usually excursions beyond the elastic range are not permitted under normal
conditions of operation; the extent of the permanent damage the structure may
sustain when subjected to extreme conditions such as blast or earthquake
loading is frequently of interest to the design engineer.
Consider a one storey structural steel shear frame subjected
to a horizontal static force F as shown in Fig. 7.19. Assume the girder
to be infinitely rigid
compared with the column, when
the load is applied numerically with the frame. Plastic hinges will eventually
form at the end of the columns. The plot of resistance versus displacement
relationship is given in Fig. 7.20a. The behaviour is linear up to the point ‘a’
corresponding to resistance R_{y} where first yielding in the
cross section occurs. As the load is increased, the resistance curve becomes
nonlinear as the column cross-sections plasticize under the system softens. The
full plastification of the cross-section occurs at point ‘b’
corresponding to maximum resistance R_{m}. Upon unloading, the
system rebounds elastically along the line ‘bc’ parallel to initial
linear portion ‘Oa’. The system remains elastic until first yielding
again attained at point ‘d’ corresponding to resistance R_{y}.
As the load is increased, further plastic hinges reform at - R_{n}
corresponding to point ‘e’. Unloading will be linear elastic parallel to
line ‘cd’.
If the maximum positive and
negative resisting forces R_{m} and -R_{m} are
numerically equal, the hysteresis loop formula by the cyclic loading will be
symmetric with respect to origin. For each cycle, energy is dissipated by an
amount that is proportional to the area within the hysteresis loop. The
behaviour illustrated in Fig. 7.20a is often simplified by assuming linear
behaviour up to the point of plastification. This type of behaviour is referred
to as ‘elasto-plastic’. This slope of the elastic loading and unloading
curve is proportional to stiffness. The elasto-plastic behaviour can be
idealized shown in Fig 7.20b. One hysteresis loop is discussed below.
Stage 1 Elastic loading
This is defined by segment ‘oa’
on the resistance displacement curve 0 ≤
u ≤ uel and u˙
> 0 where uel = Rm/k.
The resisting force as the stages is given by
Rm = Kx -- - - - - -7.78
Unloading in the stage occurs when u˙ < 0.
Stage 2 Plastic loading
This stage is represented by the segment ab on the
resistance-displacement
curve and corresponds to the
condition uel < u <u_{max} and u˙
> 0 where u _{max} is the displacement in hysteresis loop.
The resisting force in this stage is
given by
F_{s} = R_{m - -- - - - - -}7.79
Stage 3 Elastic rebound
This stage is defined by the segment bc on the
resistance-displacement
curve and corresponds to a condition (u_{max}-2
uel) < u < u_{max} and u˙ <
0. The resistance is given by
It is to be noted that load
reversal in this stage occurs than u < ( u
_{max} �' 2u˙)
and u˙ < 0 .
Stage 4 Plastic loading
The system response in this stage
is represented by segment ‘cd’ on the resistance displacement curve and
corresponds to condition u_{min} < u < (u_{max}
-uel) and u˙ < 0 where u_{min}
is the minimum displacement as the stage. The system resistance is given by
Stage 5 Elastic rebound
Once the cycle of hysteresis is completed, the system unloads
elastically along segment ‘de’. This stage corresponds to the condition u_{min}
< u < (u_{min} + 2uel) and u˙
> 0. This resisting force is given by
F_{s} = k(u
- u_{min}) - R_{m -- }7.82
_{ }
Program 7.10: MATLAB
program for dynamic response for elasto-plastic SDOF system
%program
for elasto-plastic analysis
%simulates nonlinear response of
SDOF using %elasto-plastic hysteresis loop to model
%spring resistance. The program
uses Newmark-B integration scheme clc;
k=1897251;
m=43848;
c=7767.7;
%for
elasto-plastic rm=66825.6 and
%for elastic response rm is
increased such that rm=6682500.6 rm=66825.6;
tend=30.0;
h=0.02;
nfor=1500;
%earthquake
data
%d=xlsread(‘eldat’);
%for
i=1:nfor
% tt(i)=d(i,1);
% ft(i)=d(i,2);
% ft(i)=m*ft(i)*9.81;
%end
%force data d=xlsread(‘forcedat’)
for i=1:nfor
tt(i)=d(i,1)
ft(i)=d(i,2) end
ic=1;
for t=0:h:tend for i=1:nfor-1
if (t >= tt(i)) & (t < tt(i+1))
p(ic)=ft(i)+(ft(i+1)-ft(i))*(t-tt(i))/(tt(i+1)-tt(i));
ic=ic+1; continue
end continue continue
end end x(1)=0; x1(1)=0;
x2(1)=p(1)/m;
xel=rm/k;
xlim=xel; xmin=-xel; a1=3/h; a2=6/h; a3=h/2; a4=6/(h^2);
kelas=k+a4*m+a1*c
kplas=a4*m+a1*c
ic=2;
loop=1;
for t=h:h:tend-30*h if loop==1
[x,x1]=Respond(kelas,p,x,x1,x2,m,c,ic,a2,a3,a1);
r=-rm-(xmin-x(ic))*k; x2(ic)=(p(ic)-c*x1(ic)-r)/m;
if x(ic) >= xlim loop=2;
end ic=ic+1; continue
elseif(loop==2)
loop;
[x,x1]=Respond(kplas,p,x,x1,x2,m,c,ic,a2,a3,a1);
r=rm; x2(ic)=(p(ic)-c*x1(ic)-r)/m; if x1(ic)<=0
loop=3;
xmax=x(ic); xlim=x(ic)-2*xel; end
ic=ic+1; continue elseif(loop==3)
loop;
[x,x1]=Respond(kelas,p,x,x1,x2,m,c,ic,a2,a3,a1);
r=rm-(xmax-x(ic))*k; x2(ic)=(p(ic)-c*x1(ic)-r)/m;
if x(ic)<=xlim loop=4;
end ic=ic+1;
continue elseif(loop==4)
loop;
[x,x1]=Respond(kplas,p,x,x1,x2,m,c,ic,a2,a3,a1); r=-rm;
x2(ic)=(p(ic)-c*x1(ic)-r)/m; if x1(ic)>=0
loop=1
xlim=x(ic)+2.0*xel;
xmin=x(ic); end ic=ic+1;
continue end
end ic=ic-1; for i=1:ic
tx(i)=(i-1)*h; xx(i)=x(i);
end plot(tx,xx) hold on xlabel(‘
time’)
ylabel(‘ displacement in m’)
title(‘ Displacement response history ‘)
function
[x,x1]=Respond(k,p,x,x1,x2,m,c,ic,a2,a3,a1)
dps=p(ic)-p(ic-1)+x1(ic-1)*(a2*m+3*c)+x2(ic-1)*(3*m+a3*c); dx=dps/k;
dx1=a1*dx-3*x1(ic-1)-a3*x2(ic-1);
x(ic)=x(ic-1)+dx; x1(ic)=x1(ic-1)+dx1;
Example 7.7
A shear frame structure shown in
Fig. 7.4a is subjected to time varying force shown in Fig. 7.21. Evaluate the
elastic and elasto-plastic response of the structure by the Newmark method
without equilibrium iterations.
Solution
k =
1897251N/m
m = 43848kg
c = 34605.4Ns/m
R_{m} =
66825.6N
∆t = 0.001
Figure 7.22 shows the
displacement due to loading for elastic and elasto-plastic cases.
Response spectra by
numerical integration
Construction of response spectrum by analytical evaluation of
the Duhamel integral is quite tedious. To develop response spectrum by
numerical integration,
consider a family of SDOF oscillations shown in Fig. 7.23.
Each oscillator has different natural period and frequency:
T1< T2 < T3ÉÉTn 7.83
Specify a function F(t).
The dynamic magnification factor can be calculated as
Finally DMF is plotted against desired spectrum curve.
Example 7.8
Construct a response spectrum for
the symmetric triangle shown in Fig. 7.24a. Plot DMF vs td/T in
the integral 0 ≤ td/T
≤ 10. Assume ρ = 0, td = 2s. A MATLAB program for drawing response
spectrum is given in Chapter 6. The response curve is shown in Fig. 7.24b.
Solution
A family of response spectrum curves or response spectra can
be produced for a specific load case by evaluating response maxima for several
values of damping ρ. Hence
maximum response of a linear SDOF to a specified time depends on ω_{n} and ρ.
Numerical method for
evaluation of the Duhamel integral
1.For an undamped system
The response of an undamped SDOF system subjected to a general
type of forcing function as given by the Duhamel integral as
The variation of displacement
with time is of interest. Time is divided into a number of equal intervals.
Each duration can be taken as ∆t and the
response at these sequences of time can be evaluated.
Applying Simpson’s rule
A(t -
2∆t) is the
value of the integral at time instant (t - 2∆t) by
summation of previous values. B(t) can be obtained in a
similar way.
Example 7.9
A water tank shown in Fig. 7.25a
is subjected to a dynamic load shown in Fig. 7.25b. Evaluate numerically using
the Duhamel integral for the response.
Solution
mass = 400 �-
1000/9.81 = 40774.7kg k = 35000 �- 1000N/m
The values of A and B are evaluated as shown in
Table 7.15.
2.For an under-damped system
Again numerical integration by
Simpson’s rule can be carried out to find the values of A(t) and B(t)
and hence u(t).
Selection of direct
integration method
It is apparent that a large
number of different numerical integration methods are possible by specifying
different integration parameters. A few of the most commonly used methods are
summarized in Table 7.16.
For SDOF systems the central
difference method is most accurate, and the linear acceleration method is more
accurate than the average acceleration method. It appears that the modified
average acceleration method, with a minimum addition of proportional damping is
a general procedure that can be used for the dynamic analysis of all structural
systems.
The basic Newmark constant
acceleration method can be extend to nonlinear dynamic analysis. This requires
that iterations must be performed at each time step in order to satisfy
equilibrium. For multiple degrees of freedom, which will be seen in later
chapters incremental matrix must be formed and triangularized at each iteration
or at selective point of time. Many different numerical tricks including
element by element methods have been developed in order to minimize
computational requirements.
Table 7.16 Summary of Newmark
method modified by δ factor
_{ }
Summary
For earthquake analysis of linear
structures it should be noted that direct integration of the dynamic
equilibrium is normally not numerically efficient as compared to mode
superposition method (for MDOF). The Newmark constant acceleration method with
the addition of very small amount of stiffness proportional damping is
recommended for dynamic analysis of nonlinear structure systems.
In the area of nonlinear dynamic
analysis one cannot prove that any one method will always converge. One should
always check the error in the conservation of energy for every solution
obtained.