EARTHQUAKE ENGINEERING AND STRUCTURAL DYNAMICS, VOL.

12, 4 8 1 4 9 7 (1984)

THE GENERATION OF SPECTRUM COMPATIBLE

ACCELEROGRAMS FOR THE DESIGN OF NUCLEAR POWER
PLANTS

ANDRE PREUMONT*
c/o Belgonucleaire, Rue du Champ de Mars, 25 1050, Brussels, Belgium

SUMMARY
Based on the main features of the computer program THGE, the paper reviews the techniques which are available for the
construction of an accelerogram whose response spectrum matches a design spectrum. First, a sample accelerogram is
generated as the product of the stationary random sequence by a deterministic shape function. It is assumed that the
spectral properties of the stationary process have been determined in a previous step, in order to lead to the design
spectrum as expected maximum responses of a set of single degree of freedom oscillators. The principal properties of the
Fast Fourier Transform which is used for this first step are reviewed. The paper then describes the procedures which are
available, both in the frequency domain and in the time domain, to improve the agreement between the response
spectrum and the target. It also discusses some related issues, such as the response spectrum calculations, the statistical
dependence between the three earthquake components, the duration of the time history, the variability of the secondary
response to various samples and the generation of an accelerogram whose response spectra envelop a set of design
spectra. The point of view adopted is the one of the structural engineer.

1. INTRODUCTION
Design response spectra constitute the traditional means of specifying ground motion for seismic design,
’
especially for nuclear power plants. In many instances they can be used directly to estimate the maximum
response of a linear structure. There are a number of cases, however, where this approach is no longer
possible and a time history analysis is necessary. In these cases, the generation of accelerograms of given
response spectra is necessary prior to the dynamic analysis. Examples of applications are
(i) floor response spectra calculations,
(ii) qualification of sensitive equipments, and
(iii) time integration of models involving non-linearities.
The problem of constructing an accelerogram from a given response spectrum has been treated by many
authors in various (additional references can be found in Reference 43). Basically, all the available
methods start from a sample accelerogram, the spectrum of which is as close as possible to the target
spectrum. This sample is then altered iteratively, in order to improve the agreement between its spectrum and
the target spectrum. The first sample can be obtained in one of the following ways.
(1) To select an accelerogram from a library”, whose spectrum is close to the target spectrum. This
is difficult for smoothed (broadened) spectra and floor response spectra.
(2) To use Hudson’s analogy” and consider the zero damping velocity spectrum as if it were the Fourier
spectrum of the accelerogram (that is, the modulus of its Fourier transform).?
(3) Using the random vibration theory to connect the response spectrum to the power spectral density
(PSD) of the accelerogram. This point of view has been pioneered in Reference 7; it has also been
adopted in this study.

* Consultant.
t It is demonstrated in Reference 12 that the zero damping velocity spectrum is in fact an upper bound to the Fourier spectrum.
0098-8847/84/04048 1-1 7$01:70 Received 2 June 1983
@ 1984 by John Wiley & Sons, Ltd. Revised 21 October 1983 and I 1 January 1984
482 A. PREUMONT

It is well known that the response spectrum offers a very intricate and incomplete account of the statistical
properties of the underlying stochastic process. However, from the narrow-band properties of the oscillator
involved in its definition, the response spectrum has strong links with the spectral content of the ground
acceleration. It has been pointed out7*l 3 that the statistical properties of the ground acceleration can be
entirely deduced from the response spectrum under the assumptions that the ground acceleration is
(i) Gaussian, and
(ii) quasi-stationary [i.e. the product of a stationary process xo(t)by a given deterministic shape function
f’(t)l.
In this case, the stochastic process is completely determined by the PSD Oxx(w) of the stationary process
xo(t). An algorithm for the computation of (Dxx(w)from the response spectrum S,(o) has been presented in
Reference 13; it will not be reproduced here (see also Reference 9).
Obviously, when facing the task of generating a stationary random sequence of given spectral content, it is
convenient to work in the frequency domain, i.e. with a Fourier approach. As a result, the time domain
(ARMA) models for earthquake simulation will not be reviewed in this paper. The reader is referred to
References 14, 15 and 16 and the work of Box and Jenkins. l 7 The aim of this paper is to review the techniques
which are available for the construction of a sample accelerogram of given PSD and for the improvement of
the agreement between its response spectrum and the target spectrum. Aspects which are believed to be
unfamiliar to many structural engineers are briefly reviewed. The paper also discusses some related issues,
such as the response spectrum calculations, the statistical dependence between the three earthquake
components, the duration of the time history, the variability of the secondary response to various samples
and the generation of an accelerogram whose response spectra envelop a set of design spectra.

2. GENERATION OF GROUND ACCELERATION

2.1. Shape function
In what follows, it will be assumed that
(i) the ground acceleration random process is the product of a Gaussian stationary process xo(t) and a
deterministic shape function f(t),
(ii) the PSD OXx(w) of the stationary process xo(t) has been determined according to the procedure
described in Reference 13 to correspond to the target response spectrum.
Possible shape functions are
o<u<p
f ( t ) = u,(e-“‘-e-fl’),
O,<t,<T,
or

In Sections 2.2 to 2.4 the background for the application of the Fast Fourier Transform to the generation
of a sample of a stationary Gaussian process xo of given PSD is briefly reviewed. Readers familiar with this
field are advised to skip to Section 2.5.

2.2. Sampling theorem

The sampling (or Shannon’s) theorem states that if a signal does not contain any frequency component
above a cut-off frequency fc, all the information is contained in the value of the signal sampled at a rate
f,22fc. This means that the original signal could be recovered from the sampled signal by a band-limited
interpolation (see Reference 18, p. 50). This theorem fixes the minimum length of the discrete sequence which
is necessary to represent a ground acceleration of duration To and cut-off frequency f,:
N 2 2 h To (3)
 SPECTRUM COMPATIBLE ACCELEROGRAMS 483

2.3. Digital Fourier Transform ( D F T )

As observed earlier, in order to generate a stationary random sequence of given spectral content, it is
convenient to work in the frequency domain, that is with Fourier series, or using the Digital Fourier
Transform (DFT) which is the discrete counterpart of the Fourier transform. This is not the place to discuss
the DFT, nor the Fast Fourier Transform (FFT) algorithms which reduce drastically the number of
arithmetic operations to be carried out to compute a DFT sequence; the reader is referred to the standard
literature on Digital Signal Processing." 23 It may be useful, however, to remind the reader of the following.
~

(a) Dejinition. The D F T of the finite length sequence x(m), m = 0, ..., N - 1, is defined as
1 N-1
C,(k) = - C ~ ( mWkm,) k = 0, ..., N - 1 (4)
N m=O
where W = e-jZrriN.
The inverse relationship (IDFT) is
N-1
x(rn) = C,(k) W P k m ,rn = 0, ..., N - 1
k=O

The sequences C,(k) and x(m) defined by equations (4) and (5) are both N-periodic [C,(k) = C,(k+N),
+
x(m) = x(m N)].
(b) Connection with the Fourier series and the Fourier transform. Unlike the Fourier transform, the DFT is
performed on finite length (truncated) sequences. Under the two following conditions:
(1) the signal is not modified by the truncation;
(2) the signal is band-limited and has been sampled according to Shannon's theorem;
there is a complete identity between the Fourier series coefficients ak, the DFT coefficients C,(k) and the
sampled values (in the frequency domain) of the continuous Fourier transform of the continuous signal

TO
where ak is the kth coefficient of the Fourier series under standard exponential form and x (.) is the
continuous Fourier transform of the continuous signal. Strictly speaking, since a finite duration (< To)signal
cannot be band limited, the above equality is verified only if the signal is periodic and if the duration of the
truncating function is an integer multiple of the periodicity. Any departure from the foregoing conditions
introduces leakage and aliasing errors in equation (6).
(c) Some properties. The DFT has the following properties.
If the sequence x(m) is real and if N is even

C -+i
:x
( ) =C*--i,
x(; ) i=O, ...,N/2 (7)

where * stands for the complex conjugate.

C,(O) is the mean of the sequence x(m).
The D F T spectrum is given by I C,(k) 1;' it provides a measure of the spectral content of the sampled signal
x(ml Parseval's theorem reads
(8)

(d) The Fast Fourier Transform. The FFT algorithm consists of a particular organization of the
computation of the series (4) and (5). For N = 2", the number of arithmetic operations is reduced from N2
(complex multiplication +addition) to a number proportional to only (N/2)log2N. For N = 1,024, this
amounts to a reduction by a factor of about 200! The algorithm has been generalized for any value of N . For
real sequences, additional simplifications come from equation (7), which allows computation of the DFT of a
2N points real sequence with an N points FFT algorithm (see Reference 19, p. 167). FFT subroutines are
available in the 1iterat~x-e.~~
484 A. PREUMONT

From the identity between the DFT and the Fourier series stated by equation ( 6 )and from the tremendous
reduction of the amount of computation that can be gained from the FFT algorithm, the DFT appears as a
very effective tool for the generation of a stationary random sequence of given power spectral density. From
Parseval’s theorem, [equation (S)] and the corresponding form for Qxx(o), in order to fit the given PSD
Qxx(o), the module of the D F T coefficients of the sequence must be taken according to
C,(O) = 0 (zero mean)
(9)
I C,(k)I = [ @ x x ( k ~ ~ o ) ~ ~ >k o=] *1,, ..., N / 2
where coo = 2n/To is the fundamental circular frequency. The phases of the C,(k) will now be chosen in order
to assure that the stationary sequence be Gaussian.

2.4. Central limit theorem

A simple version, known as the equal components case of the central limit theorem,24 states that if we
consider a set of identically distributed, mutually independent random variables, the probability distribution
function of the sum of those random variables, properly normalized, tends to be that of a Gaussian random
variable with zero mean and unit variance. The theorem can be extended to the sum of independent random
variables having different probability distribution functions, providing none of them contributes predomi-
nantly to the overall behaviour of the p h e n ~ m e n o n . ~ ~
If we choose the phases of the DFT coefficients as independent random variables of uniform distribution
between 0 and 2n:
independent random variables
19, = arg[C,(k)] = uniformly distributed in [0,2n], k = 1, ..., N / 2

for any given t, all the frequency components of the signal will be mutually independent random variables.
According to the central limit theorem, the sum of a large number of such components will be Gaussian. It is
to be noted that the Gaussian distribution will be achieved whatever the distribution of the phases. The
uniform distribution, however, agrees with observations on actual records.26 Computer independent
random number generators are available in the literat~re.~’. 48

2.5. Non-stationary ground motion

Equations (3), (7), (9) and (10) above give a complete characterization of the D F T of a real stationary
Gaussian random sequence of specified power spectral density. Equation (3) assures an adequate sampling
rate to avoid aliasing, equation (7) assures that the sequence is real, equation (9) gives the proper frequency
content, while relation (10) assures a Gaussian distribution. Any set of independent random phases (10) will,
after IDFT, lead to a sample of a stationary sequence of appropriate statistical properties, xo(t).The non-
stationary accelerogram is obtained simply by multiplying the stationary sequence by the given shape
function f(t) [equation (1) or (2)]
i(t) = X0(t).f’(t) (11)
However, observations of actual records tend to show that not only the amplitude but also the frequency
content of the ground motion changes with time. More specifically, it has been suggested3’ that the rate of
zero up-crossings vo and the rate of maxima v1 as given by Rice’s formulae32 vary according to
v i = qiexp(-yit), i = 0, 1
where q iand y i are positive constants which can vary considerably from one record to another. These results
indicate that high frequency components are concentrated at the beginning of the accelerogram. It is difficult
to assess how this feature affects the response of structures and, from the present knowledge of the issue, it
seems that accelerograms of the quasi-stationary form given by equation (11) are quite sufficient for
engineering purposes. Non-stationary frequency content has been modelled in various ways.16.3 1 , 51- 5 2 In
connection with the present discussion, it can be simulated very simply in the following way.
 SPECTRUM COMPATIBLE ACCELEROGRAMS 485

Several shape functions fi(t) are selected, each of them being appropriate for a subset [mi-I, oil of the
overall frequency range [0, o,].For each of those intervals, a DFT sequence is generated by selecting the
DFT coefficients C,(k) according to equations (9), (10) and (7) for wave numbers corresponding to circular
frequencies ko,pertaining to the interval [mi- mi]and C,(k) = 0 outside the interval. The IDFT provides a
stationary sequence which is multiplied by the shape function f;(t). The final accelerogram is obtained by
adding the contributions of the various intervals [ m i - oil,as a result of the linearity of the DFT.

2.6. Additional features

Two additional modifications have to be carried out on the accelerogram as obtained by equation (1 1).
(1) Maximum acceleration. The maximum acceleration of the non-stationary sequence constructed by
equation (11) varies from one sample to another and does not match exactly the maximum acceleration amax
which, in the current design practices for nuclear power plants, is specified for the site (high frequency
asymptote of the response spectrum, for all damping values). In order to set the maximum value of I x(t) I to
amaxwith the slightest possible alteration to the frequency content of x(t), the following approach is used,
depending on whether 1 x(t)Imax is larger than amaxor not.
If I x(t) l m ax < amax,all the values of the sequence in the lobe where I i ( t ) I is maximum are scaled up
according to the ratio amax/lx ( t ) lmax (Figure 1).
If 1 ?(t) Imax > amax,the values of the sequence, in each lobe where I x(t) I exceeds amax,are scaled down in
order to reduce the maximum value in this lobe to amax.

-amax

Figure 1. Procedure for setting the maximum acceleration to amax

( 2 ) Baseline correction. Although leading to accelerograms which are very similar to actual records, the
above procedure leads to velocity and displacement diagrams which may depart considerably from those of
real earthquakes. In particular, it may lead to unrealistic maximum and end values for the displacement and
the velocity (the end velocity should obviously be zero). In order to remove this anomaly without much
alteration to the accelerogram, a polynomial baseline correction is added to the acceleration
+
q t ) = X(t) a, + a , t -t a2 t 2 (12)
where the coefficients a,, a , and a2 are computed in order to minimize the integral of the square of the
velocity k,(t).2 7 , 2 8 This correction induces important changes in the ground velocity and displacement
486 A. PREUMONT

diagrams, without any noticeable alteration to the accelerogram, or to the response spectrum (except for very
low frequenciesz9).The end velocity obtained by this method is generally very small. Besides, the values of
the non-dimensional parameter r defined as
am ax L a x
r=
omax

obtained with artificial accelerograms after baseline correction, are well within the range of the values
obtained for real earthquakes, 5 < r < 15 (Reference 30, p. 233).

2.7. Example
As an example of the above procedure, Figure 2 shows an artificially generated accelerogram, together
with its 2 per cent damping response spectrum. The shape function has been taken according to equation (2),
with t , = 2, t z = 15, n = 1 and CL = 3/5;the overall duration is To = 20 s. The power spectral density (Dxx(w)
is shown in Figure 3; it has been determined by the procedure developed in Reference 13 to correspond to the
2 per cent damping NRC horizontal spectrum,' with a maximum acceleration of lg. The NRC spectrum is
also shown in Figure 2; one sees that the artificially generated accelerogram is fairly close to the target.

10- Period (sec) 1u

Figure 2. Artificially generated accelerogram and its 2 per cent damping pseudovelocity response spectrum, together with the 2 per cent
damping NRC horizontal spectrum

About the PSD of Figure 3, the following is worth noting.

(1) The behaviour of the PSD for periods above 1.8s has no physical meaning. This is a direct consequence
of the stationary approach used in connecting the response spectrum to the PSD.13 This approach
holds providing the time constant controlling the non-stationary response due to zero initial condition
is much smaller than the duration of the accelerogram, that is if (2c0,)- << To. In the present case
(To = 20s, ( = 0.02), this implies natural periods satisfying T,<<4n<T0= 5.03s. Note that this spurious
 SPECTRUM COMPATIBLE ACCELEROGRAMS 487

Figure 3. Power spectral density used to generate the accelerogram of Figure 2

behaviour of QXx(o)
has no practical influence on the generated accelerogram, for only a small number
( - 10) of Fourier components pertain to that frequency range.
(2) In relating the response spectrum to the PSD, the duration of the accelerogram is implicitly contained
in the peak factor (see Reference 44, for example), whose value depends on the number of equivalent
half-cycles (that is, twice the product of the duration by the central frequency) and the bandwidth of the
system (that is, the damping 5).

3. ITERATIVE IMPROVEMENT OF THE SOLUTION

In the previous section, we have seen how the random vibration theory can be applied to generate a sample
accelerogram of specified statistical properties. We now turn to the ways this accelerogram can be modified
in order to improve the agreement between its response spectrum Sv and the target spectrum S:arge'.

3.1. Step 1: in the frequency domain

Considering the narrow-band character of the oscillator involved in the response spectrum definition and
the connection between the probabilistic sense of this latter to the so-called first crossing problem (see
Reference 13), some improvement can be obtained by working in the frequency domain: in fact, from the
white noise approximation of the variance mo of the response of a lightly damped single degree of freedom
(s.d.0.f.) oscillator to a wide band excitation aXx(w), we have
SV(0,) -%(df
- @xx(wJf -C,(W (13)
Consequently, each frequency component of the accelerogram will be modified according to the ratio of the
target spectrum to the response spectrum. This means that a new DFT sequence will be constructed by
488 A. PREUMONT

The new time sequence x(m) will result by IDFT. In computing (14), the response spectrum ordinates S,(w,)
need not be calculated for every frequency w;, they can be interpolated, providing available spectrum values
are close enough to make this interpolation meaningful. This point is connected to the bandwidth of the
oscillator, 2t0,; it will be discussed in Section 4.
Figure 4 shows the result of three successive applications of equation (14) to the accelerogram of Figure 2.
The baseline correction has been applied each time. The procedure improves considerably the agreement
between the response spectrum and its target, particularly in the intermediate frequency range (say from 1 to
lOHz), where relation (13) holds best. It leaves, however, some peaks and troughs that cannot be removed,
mainly in the low frequency range, where the DFT coefficients are very scarce (their frequency spacing is
equal to the inverse of the total duration of the signal). The agreement can be further improved by the
following procedure, first introduced by Kaul. 3 3

S"

lo" lo" Period (XC)

I
10'
Figure 4. Accelerogram and response spectrum after three iterations of step 1 of the improvement procedure

3.2. Step 2: in the time domuin

So far, we have considered only the response spectrum ordinates; two additional pieces of information will
now be considered: the maximum deformation time, that is, the time at which the absolute relative response
of the oscillator is equal to the response spectrum ordinate, and the sign of the maximum relative response.
Figure 5 shows a plot of the maximum deformation time corresponding to Figure 4. It is worth noting that
for large periods, the maximum deformation time can be larger than the strong motion duration. Similar
observations have been made on actual records.42
From now on, for each frequency, the response spectrum will be characterized by the following
information:
( 1 ) S,(UI,)= displacement response spectrum [S,(o,) = orSd(ok)];
( 2 ) t,(w,) = maximum deformation time;
 SPECTRUM COMPATIBLE ACCELEROGRAMS 489

in0
lo-' IU 10'
Period(sec)
Figure 5. Maximum deformation time corresponding to Figure 4

(3) s(o,)=integer whose value is 0 if the maximum relative displacement is positive and 1 if it is
negative.
The procedure is based on the assumption that, for each natural frequency, the maximum deformation time
is not affected by small perturbations AX of the accelerogram. After having selected the frequencies
wi (i = 1, ..., n) at which the response spectrum has to be corrected, the amount by which the displacement
must be altered at time t,(wi) is computed by
Ad(wi) = [S~rge'(oi)-Sd(Wi)](-l)S(W') (i = 1, ...)n) (15)
Then, the problem can be formulated as follows. Find a correcting wave A i ( t ) leading to

Ai(t).h[t,(w,)-t] dt = Ad(wi) (i = 1, ..., n)

where
1
h(t) = --exp ( -mi Ct)sin odi
t
wdi

is the impulse response of the s.d.0.f. oscillator of frequency wi and damping ([adi= w i J ( l - (*)I.The left
side of equation (16) is the well-known convolution integral, giving the oscillator response at time t,(wi). The
right side represents the desired displacement correction, given by equation (15). The integral equation (16)
can be transformed into a system of linear algebraic equations by choosing the perturbation AX as a linear
combination of appropriate independent time functions g j t ) :
n
AX(t) = C1 a j g j ( t )
j=

Introducing this form into equation (16) leads to a system of algebraic equations for the coefficients a j . In the
choice of the functions gj(t), the following aspects have to be considered.
490 A. PREUMONT

(a) AX will be as small as possible.

(b) When applied to the oscillator of frequency wP the function y j t ) will lead to a response which is
maximum at time t = tm(wj).
(c) When applied to an oscillator of frequency wi # w j , the function g,(t) will lead to a response which is as
small as possible, so that the system of equations is as close as possible to the diagonal form. Ideally,
each spectral ordinate Ad(oi)could be tuned independently through the corresponding coefficient ai.
(d) The choice of gi(t)should allow the matrix coefficients of the linear system to be computed analytically.
Taking these aspects into account, a good choice for g i ( t ) is the impulse response of the s.d.0.f. oscillator of
circular frequency mi and damping 5, back in time, starting from the maximum deformation time tm(wi)
(Figure 6):
gi(t) = h[t,(oJ - tl, o< t 6 t,(wJ
(19)
=o outside

Figure 6. g,(~),ith individual component of the correcting wave, together with the function involved in the convolution leading to the
coefficients A,,

The reader is referred to D r e n i ~ kfor

~ an
~ .extensive
~~ discussion of the critical feature of this excitation. The
corresponding linear system is

j= 1
f k i j a j= A~(WJ, i = I , ..., n (20)
where
m i n [ h , tmj]
A , .= ~

exp [ - <mi ( t m i - 211 exp K - 5wj(tmj - 711

I
Qdi.Wdj S0
x sin wdi(tmi- t)sin wdj(tmj
- t)dz (21)
and tmi has been used in short for tm(wi).From equation (21), we see that the matrix A , is symmetric and has
its dominant terms on the diagonal. In fact, a function gi(t) taken according to equation (19) will lead to a
large response for the oscillator of frequency wi;its contribution to the response of oscillators having either
their natural frequency or their maximum deformation time far apart is expected to be small (Figure 6).
In implementing the above procedure, we have found it useful to multiply the correcting wave (18) by a
window function w(t) (Figure 7):

lcr 0.1 To6 t 6 To

 SPECTRUM COMPATIBLE ACCELEROGRAMS 49 1

Figure 7. Window function multiplying the computed correcting wave

whose aim is
(i) to provide a smooth start, and
(ii) to offer the possibility to attenuate the computed correcting wave by means of a relaxation coefficient
cI< 1.
This latter has been prompted by the need to prevent any strong departure from the basic assumption of
unchanged maximum deformation time. In fact, such departures inevitably occur when AX(t) becomes too
large; they may lead to the divergence of the whole process. The corrected accelerogram reads
X(t)new = Z i t ) + ~ ( t AX(t)
). (23)
The procedure can be repeated until satisfactory agreement between the response spectrum and the target is
reached. A typical correcting wave is seen in Figure 8; note that its maximum acceleration is one order of
magnitude smaller than the accelerogram to be corrected. The result of several applications of the above
procedure on the accelerogram of Figure 4 is shown in Figure 9. The procedure has been found very effective
in removing the peaks and troughs from the response spectrum. It is the author's experience, however, that it
is important to have some control between the successive applications of the procedure. In that way, all the
information contained in Sd(wk),tm(uk) and s ( q ) can be used in a sensible way in the choice of the frequencies
to be corrected and the relaxation parameter 8 . Control of the successive applications of the procedure is
particularly useful when one has to generate an accelerogram matching smoothed response spectra for two or
more damping values. Table I gives the evolution of the error between the response spectrum and the target
for the various examples discussed previously. It is to be emphasized that these examples should by no means
be considered as 'good examples'. No doubt the initial agreement (Figure 2) could be better by chance. The
agreement between the spectrum and the target in Figure 9 could also be improved by additional application
of the correction procedures.

0.6,

4. RESPONSE SPECTRUM CALCULATIONS

The calculation of the response spectrum involves the time integration of the differential equation
i + 2 ( 0 , k +0,'x = -Xo(t). (24)
492 A. PREUMONT

Figure 9. Accelerogram and response spectrum after several applications of step 2 of the improvement procedure

Table I. Relative error associated with the artificially generated

accelerograms
% points Maximum RMS
error error error
>lo% (%I (%)
Sample generated according t n 56 84.5 26.9
random vibration theory
(Figure 2)
Iteration 1 of step 1 45 55-5 15.6
Iteration 2 of step 1 39 57.2 14.5
Iteration 3 of step 1 (Figure 4) 33 56 14
After several iterations of 4.5* 12.7 4.7
step 2 (Figure 9)
*All corresponding to frequencies below 0.5 Hz.

Most of the operators available for this purpose can be written, in terms of the state variables x and i,
as

For Duhamel's formula, which is the most popular algorithm for response spectrum calculation,2R the
homogeneous part of equation (25), A(w,, 5, dt), results from a rigorous (analytic) integration of equation (24).
As a result, it is devoid of stability problems, as well as of periodicity and amplitude errors in the usual sense
(see Reference 34, Chapter 9, for example). The error associated with Duhamel's formula is connected with
the non-homogeneous part of equation (25), B(w,, <,dt), that is, with the assumption made on the variation
of x o ( t )between the sampling points. It has been studied in Reference 35, where the details of the matrices
 SPECTRUM COMPATIBLE ACCELEROGRAMS 493

A ( . )and B ( . )can be found for the popular piecewise linear assumption. This latter leads to what appears as
the best one-step scheme for response spectrum calculation. It must be pointed out, however, that the
sampling rate prescribed by Shannon's theorem [equation (3)] gives an upper bound to the time step in order
to avoid aliasing; it does not ensure an accurate result of the integration operator (25). Results of Reference 35
suggest that a sampling rate 3 to 5 times as large may be necessary from the point of view of accuracy. In the
present study, the sampling rate of the accelerograms was either f, = 102 Hz (j& = 3.1) or f, = 204 Hz
( f , l f , = 6.2); the effect of the sampling rate on the calculated response of high frequency correcting
components g i ( t ) (say f i > 25 Hz) was clearly noticeable.
The second important aspect when calculating response spectra is the spacing between control frequencies.
This aspect has been discussed in Reference 5, where a minimum of 100 frequency points was recommended.
Control frequencies are also recommended in Appendix N of the ASME code;j6 they are conspicuously
independent of the damping. The problem can be stated simply as follows: what is the maximum spacing
between control frequencies ensuring that the response spectrum does not contain any hidden spike
involving a major departure from the target? As the s.d.0.f. oscillator response to a wide band random
excitation is controlled mainly by the energy of the excitation contained within the bandwidth 2 5 0 , of the
oscillator, large differences in the response from one control frequency to the next will be prevented by
overlapping oscillators, that is, if the control frequencies are chosen according to
wn < W n + 1 < 0Al+ 25) (26)
From this result, one sees that
(1) the control frequencies should be evenly spaced on a logarithmic scale,
(2) the total number of control frequencies, M , should be taken according to

where fo and f c are the extreme frequencies of the spectrum and 5 is the damping. The number of
control frequencies should be proportional to the inverse of the damping. For f c = 50 Hz, fo = 0.1 Hz
and 4 = 0.02, equation (27) leads to M > 150 (in the examples of Figures 2, 4 and 9, 200 control
frequencies have been used); for 5 = 0.05, the minimum number is reduced to M > 65.

5. RELATED ISSUES
Numerous guidelines have been published for the generation of synthetic accelerograms for use in the design
of nuclear power 3 6 - 3 8 Alth ough some of them are still controversial, the following aspects are
worth pointing out.

5.1. Three-dimensional ground motion

In general, two sets of response spectra are defined for any given site; the first one is appropriate for any of
the horizontal directions and the second corresponds to the vertical direction.' From them, three orthogonal
components of the ground acceleration must be derived. Based on the low correlation observed for actual
record^,^' 40 common design practices require treating the three components of the input time history as
statistically independent. This ensures the applicability of the SRSS (square root of the sum of the squares)
rule in the combination of the responses due to the three input components4' and is also essential for the
analysis of complex systems where coupling between the three directions may be significant. In the present
method, the statistical independence of time histories is automatically achieved as a result of the random
phases of the DFT coefficients [equation (lo)]. Correlation coefficients at zero lag between samples obtained
with different sets of random phases are well below 0.16, which is the maximum value recommended in
Reference 40 for time histories to be considered as statistically independent.
The foregoing current design practice deserves some comment.
494 A. PREUMONT

The zero lag cross-correlation matrix between the three components of actual records allow a set of
principal axes to be defined, that is, a change of orthogonal coordinates leading to a diagonal cross-
correlation matrix. Recent studies49*5 0 seem to indicate that this set of principal axes (along which the
earthquake components are truly uncorrelated) has a tendency to be orientated in the following way: the
major principal axis is directed towards the expected epicentre and the minor principal axis is directed
vertically. The assumption that the entire cross-correlation structure of the three earthquake components is
contained in the cross-correlation matrix at zero time lag may be questioned. However, the physical meaning
of the principal axes become more apparent if we notice that, for stationary ergodic processes, this matrix is
proportional to the intensity tensor at a point as defined by Arias.s3 Therefore, the principal axes coincide
with the principal directions of the intensity.
When designing a plant against future earthquakes, the direction of the slip fault zone, that is, the major
principal axis, is not known beforehand. However, if the one-dimensional reference ground motion specified
by the horizontal response spectrum is assumed to refer to the major principal axis, the current design
procedure described above (that is to say, assigning this reference excitation to two arbitrary orthogonal
directions x and y in the horizontal plane) is believed to be conservative, in the sense that it leads to a larger
value of the intensity on the horizontal plane53

This is a consequence of the fact that 1, is invariant with respect to the rotation of the reference frame around
the z axis.
Note that the concept of principal directions can be considered in the frequency domain as welLS4The
value of the coherence function between the various components of the earthquake excitation has also been
proposed as a criterion for the acceptability of time histories;38 actual figures, however, do not display
remarkable values (that is neither close to 1, nor close to 0).
5.2. Duration of the time history
The amount of computation involved in the time history analysis of a system is proportional to the
duration of the time history. Consequently, if it is not specified from the seismo-tectonic conditions
governing site seismicity and if fatigue or liquefaction are not involved in the analysis, the strong motion
duration should be made as short as possible, bearing in mind that it should at least allow the steady-state
structural response to be achieved. This involves a time constant z = ( 2 5 0 , ) - ’ , where w, is the lowest
significant natural frequency of the structure and 5 is the corresponding damping. A minimum strong motion
duration of 6 s (coupled with a build-up duration of about 4 s) is recommended for nuclear power plants.36.3 7
It is appropriate to recall that the maximum response of long period structures may occur at times larger
than the interval where high intensity acceleration pulses occur42 (see also Figures 4 and 5 which show that
the response spectrum diagram would be affected if the last 3 s of the accelerogram were truncated). This
implies that the overall duration may also be significant.’
As mentioned before, in relating the response spectrum to the PSD, the duration of the accelerogram is
implicitly contained in the peak factor. One knows, however, that the model fails for long periods.13

5.3. Vuriubility of the structural response to vurious samples

Since an artifically generated time history is no more than a sample from the infinite population of all the
time histories matching the response spectrum, it is instructive to investigate the variability of the structural
response to various samples. This problem has been considered in Reference 55, where the responses of an
uncoupled two degrees of freedom model to 35 artificial accelerograms have been computed. The study
indicates that, although the response spectrum of each sample is very close to the target spectrum, the
variability of the response of the secondary system (acceleration and relative displacement) is very large,
especially when the systems are tuned. Ratios between the maximum and the minimum responses exceeding
5 have been observed. From the foregoing, the conservatism associated with floor response spectra (even
broadeneds6) obtained with one single time history may be questioned.
 SPECTRUM COMPATIBLE ACCELEROGRAMS 495

5.4. Accelerogram compatible with several response spectra

It is in general not possible to generate an accelerogram compatible with a set of response spectra. The
problem must therefore be changed into the generation of an accelerogram whose response spectra envelop a
set of design spectra. As far as the author is aware, there is no general method to perform such a task. A
procedure can be outlined as follows:
(1) Compute the power spectral density corresponding to each of the response spectra and use the
envelope to generate the first sample.
(2) The first step of the improvement procedure [formula (14)] can be applied by using the maximum value
of the ratio S 7 r g e t ( ~ k ) / S Vfor
( ~the
k ) various damping ratios.
(3) The second step of the improvement procedure can be applied, in much the same way as discussed
before, to remove specific peaks and troughs. Here, the fact that one peak can be removed for one
damping value without affecting the response spectra for other damping values depends on the
corresponding maximum deformation times. This step needs engineering judgement and may be time
consuming.
The foregoing procedure is currently being implemented in the computer program THGE. Enveloping
several response spectra will involve a penalty which can be measured by the average r.m.s. error (always by
excess, except for very few control frequencies). A sample problem involving five response spectra
corresponding to damping values ranging from 0.02 to 0.4 has led to an average r.m.s. error of about 30 per
cent.

6. CONCLUSIONS
(1) Under the assumption that the ground acceleration is a Gaussian quasi-stationary random process, it
is possible to determine the power spectral density which corresponds to the response spectrum. l 3
From this power spectral density, a sample accelerogram can be generated, using the Fast Fourier
Transform. Its response spectrum is fairly close to the target spectrum.
(2) The method can be extended in a straightforward manner to take into account the non-stationary
frequency content.31 From the present state of knowledge of the effect on the structural response,
however, it is believed that the assumption of quasi-stationarity is sufficient for design purposes.
(3) The sample can be subsequently altered in order to improve the agreement between its response
spectrum and the target. Both frequency domain and time domain33procedures are available. Their
combination leads to a very good agreement, throughout the whole frequency range.
(4) The sampling rate required in order to reach a satisfactory accuracy in the response spectrum
calculations may be 3 to 5 times as large as the one required by Shannon’s theorem.35
( 5 ) Control frequencies for response spectrum calculations should be evenly spaced on a logarithmic scale.
Their number ought to be connected to the damping ratio. This aspect is not taken into account in the
current recommendation^.^^
(6) Although principal directions for the intensity have been suggested,49s5 3 the current practice of
assigning statistically independent excitations with identical response spectra to two arbitrary
orthogonal axes in the horizontal plane is believed to be appropriate for design purposes. In any event,
it is conservative if the design spectrum is assumed to refer to the major principal axis.
(7) When it is not specified, one must be careful in choosing the strong motion duration of the
accelerogram, in order to allow the structural response to reach its steady state. The effect of the total
duration may also be significant on the response of long period structure^.^'
(8) The structural response to various samples having the same response spectrum may vary consider-
ably.55 Using one single time history for the generation of floor response spectra (as is currently
practised) should therefore be considered with care.
(9) In general, it is not possible to generate an accelerogram matching several design spectra. The problem
turns into the generation of an accelerogram whose response spectra envelop a set of design spectra.
The way the procedure must be changed for this purpose has been briefly outlined.
496 A. PREUMONT

REFERENCES
1. U.S. NRC, Regulatory Guide 1.60.
2. N. C. Tsai, ‘Spectrum compatible motions for design purposes’, J . eng. mech. diu. ASCE, 98, 345-356 (1972).
3. P. C. Rizzo et a/., ‘Development of real synthetic time histories to match smooth design spectra’, S M l R T - 2 , London (1973), Paper
K1/5*.
4. R. H. Scanlan and J. Sachs, ‘Earthquake time histories and response spectra’, J . eng. mech. diu. ASCE 100, 635-655 (1974).
5. D. E. Shaw et a/., ‘Proposed guidelines for synthetic accelerogram generation methods’, S M l R T - 3 , London (197% Paper K1/4.
6. S. Levy and J. P. D. Wilkinson, ‘Generation of artificial time histories, rich in all frequencies, from given response spectra’, Nucl. eng.
des. 38, 241-251 (1976).
7. E. H. Vanmarcke and D. A. Gasparini, ‘Simulated earthquake ground motions’, S M l R T - 4 , San Francisco (1977), Paper K1/9.
8. G. Kost, T. Tellkamp, H. Kamil, A. Gantayat and F. Weber, ‘Automated generation of spectrum compatible artificial time histories’,
Nucl. eng. des. 45, 243-249 (1978).
9. A. Preumont, ‘A method for the generation of artificial earthquake aecelerograms’, Nucl. eng. des. 59, 357-368 (1980).
10. D. E. Hudson et a!., ‘Strong motion earthquake accelerograms. Digitized and plotted data’, E E R L Report, Pasadena, CA August
1972.
11. CNEN-ENEL Commission on Seismic Problems associated with the Installation of Nuclear Plants, ‘Contribution to the study of
Friuli earthquake, digitized and plotted data’, Rome, Italy, July (1976).
12. D. E. Hudson, ‘Some problems in the application of spectrum techniques to strong motion earthquake analysis’, Bull. seism. soc. Am.
s 52, 417430 (1962).
13. A. Preumont, ‘On the connection between the response spectrum and the spectral properties of earthquake accelerograms’, ICASP-
4, Florence, June (1983).
14. N. W. Polhemus and A. S . Cakmak, ‘Simulation of earthquake ground motions using autoregressive moving average (ARMA)
models’, Earthquake eng. struct. dyn. 9, 343-354 (1981).
15. M. K. Chang et al., ‘ARMA models for earthquake ground motions’, Earthquake eng. struct. dyn. 10, 651-662 (1982).
16. M. Hoshiya and T. Chiba, ‘Simulation methods of multi-dimensional non-stationary processes by time domain models’, Proc.
JSCE, No 296, 121-130, April (1980).
17. G. P. Box and G. M. Jenkins, Time Series Analysis, Forecasting and Control, Holden-Day, San Francisco, 1970.
18. A. Papoulis, The Fourier Integral and its Applications, McGraw-Hill, New York, 1962.
19. E. 0. Brigham, T h e Fast Fourier Transform, Prentice-Hall, Englewood Cliffs, N.J., 1974.
20. A. Oppenheim and R. Schafer, Digital Signal Processing, Prentice-Hall, Englewood Cliffs, N.J., 1975.
21. B. Gold and C. Rader, Digital Processing ofSignals, McGraw-Hill, New York, 1969.
22. G. D. Bergland, ‘A guided tour of the fast Fourier transform’, IEEE Spectrum, July (1969).
23. N. Ahmed and K. R. Rao, Orthogonal Transforms f o r Digital Signal Processing, Springer-Verlag, 1975.
24. W. B. Davenport, Probability and Random Processes, McGraw-Hill, New York, 1970.
25. Y. K. Lin, Probabilistic Theory of Structural Dynamics, McGraw-Hill, New York, 1967.
26. Y. Ohsaki, ‘On the significance of phase content in earthquake ground motions’, Earthquake eng. struct. dyn. 7 , 4 2 7 4 3 9 (1979).
27. P. Ruiz and J. Penzien, ‘Probabilistic study of the behaviour of structures during earthquakes’, EERC Report, College of
Engineering, University of California, Berkeley, CA, March 1969.
28. N. C. Nigam and P. C. Jennings, ‘Digital calculation of response spectra from strong motion earthquake records’, E E R L Report,
Pasadena, CA, June 1968.
29. A. K. Chopra and 0. A. Lopez, ‘Evaluation of simulated ground motions for predicting elastic response of long period structures
and inelastic response of structures’, Earthquake eng. struct. dyn. 7 , 383-402 (1979).
30. N. M. Newmark and E. Rosenblueth, Fundamentals o/ Earthquake Engineering, Prentice-Hall, Englewood Cliffs, N.J., 1971.
31. G . R. Saragoni and G . C. Hart, ‘Simulation of artificial earthquakes’, Earthquake eng. struct. dyn. 2, 249-267 (1974).
32. S. 0. Rice, ‘Mathematical analysis of random noise’, reprinted in Selected Papers on Noise and Stochastic Processes (Ed. N. Wax),
Dover, New York, 1954.
33. M. K. Kaul, ‘Spectrum consistent time history generation’, J . eng. mech. diu. ASCE, 104, 781-788, (1978).
34. K.J. Bathe and E. L. Wilson, Numerical Methods in Finite Element Analysis, Prentice-Hall, Englewood Cliffs, N.J., 1976.
35. A. Preumont, ‘Frequency domain analysis of time integration operators’, Earthquake eng struct. dyn. 10, 691-697 (1982).
36. ASME Code, Section 111, Appendix N, Dynamic Analysis Methods.
37. C . W. Lin, ‘Criteria for generation of spectra consistant time histories’, S M l R T - 4 , San Francisco (1977), Paper Klj11.
38. C. W. Lin, ‘Time history input development for the seismic analysis of piping systems’, J . press. uess. techn., trans. A S M E , 102,212-
218 (1980).
39. N. Newmark, J. Blume and K. Kapur, ‘Seismic design spectra for nuclear power plants’, J . power dio. A S C E 99, 287-303 (1973).
40. C. Chen, ‘Definition of statistically independent time histories’, J . struct. diu. ASCE 101, 449-451 (1975).
41. U.S. NRC, Regulatory Guide 1.92, December 1974.
42. M. Amin and A. H . 4 . Ang, “on-stationary stochastic model of earthquake motions’, J:eng. mech. dio. A S C E 94, 559-583 (1968).
43. P.-T. D. Spanos, ‘Digital synthesis of response-design spectrum compatible earthquake records for dynamic analyses’, Shock uib.
dig. 15(3), 21-30 (1983).
44. E. H. Vanmarcke, ‘Structural response, to earthquakes’, in Seismic Risk and Engineering Decisions (Eds. C. Lomnitz and E.
Rosenblueth), Elsevier, Amsterdam, 1976, Chapter 8.
45. R . F. Drenick, ‘Model-free design of aseismic structures’, J . eng. mech. diu. ASCE 96, 4 8 3 4 9 3 (1970).
46. R. F. Drenick, ‘Aseismic design by way of critical excitation’, J . eng. mech. dio. ASCE 99,649-667 (1973).
47. E. J. McGrath and D. C. Irving, ‘Techniques for efficient Monte Carlo simulation’, Vol. 11: ‘Random number generation for selected
probability distributions’, Oak Ridge National Laboratory, ORNL-RSlC-38, December 1974.
48. WEE Acoustics, Speech, and Signal Processing Society, Programs for Digital Signal Processing, IEEE Press, 1979.
 SPECTRUM COMPATIBLE ACCELEROGRAMS 497

49. J. Penzien and M. Watabe, ‘Characteristics of three-dimensional earthquake grounc motions’, Earthquake eng struct. dyn. 3, 365-
373 (1975).
50. T. Kubo and J. Penzien, ‘Analysis of three-dimensional strong ground motion along principal axes, San Fernando earthquake’,
Earthquake eng. struct. dyn. 7, 265-278 (1979).
51. T. Kubo and J. Penzien, ‘Simulation of three-dimensional strong ground motion along principal axes, San Fernando earthquake’,
Earthquake eng. struct. dyn. 7, 279-294 (1979).
52. H. L. Wong and M. D. Trifunac, ‘Generation of artificial strong motion accelerograms’, Earthquake eng. struct. dyn. 7, 509-527
(1979).
53. A. Arias, ‘A measure of earthquake intensity’, in Seismic Designfor Nuclear Power Plants (Ed. R. J. Hansen) MIT Press, Cambridge,
Mass., 1970, pp. 438483.
54. J. Bendat and A. G. Piersol, Engineering Applications ofCorrelation and Spectral Analysis, Wiley, New York, 1980, Chapter 7.
55. A. Kurosaki and M. Kozeki, ‘Statistical uncertainty of response characteristic of building-appendage system for spectrum
compatible artificial earthquake motion’, SMIRT-6, Paris (1981), Paper K7/7.
56. U S . NRC, Regulatory Guide 1.122, February 1978.