The mathematical problem under
consideration is one of earthquake-induced waves traveling from a
homogeneous halfspace of rock into a basin with heterogeneous soils
and irregular geometries. To completely specify the problem,
consider first a halfspace made up of the basement rock material.
Suppose there is a transient incident plane SH-wave that produces
an incident plus reflected displacement field 
, which is a
function of position and time, and is polarized perpendicularly to
the cross-sectional plane of the valley. This is the free-field
displacement. Suppose now that the valley, idealized as an
infinitely long cylinder with an irregular cross section, is
inserted in the free surface, as shown in Fig. 2. The
problem, then, is to determine the total displacement (antiplane)
field, u, within and outside the valley due to the incident field.
Since the soil material is assumed to be piecewise uniform, u satisfies the wave equation
in each subdomain and in the exterior region occupied by the rock. The displacement u and the traction are required to be continuous across each subdomain, and the traction must vanish at the free surface. The system is initially at rest.
To solve this problem, we will use the finite element method with a mesh tailored to the local wavelength of the propagating waves. Two important issues must be considered for solving wave propagation problems in infinite domains by this method. One, mentioned already in the Introduction, is the need to render the domain of computation finite and to limit the occurrence of spurious reflections. This is accomplished here by introducing the outermost circular segment I-I shown in Fig. 3a as an artificial boundary. On this arc we impose an absorbing boundary condition that is local in both space and time (Bayliss and Turkel, 1980). We use the implementation procedure proposed by Kallivokas et al (1991), which produces symmetric element damping and stiffness matrices. These can be readily incorporated into standard finite element software. Physically, this boundary condition can be interpreted as a set of springs and dashpots distributed along the arc.
The second point that requires attention is the need
to incorporate the excitation into the model when the earthquake source
is located outside the computational domain. This is carried out by
means of a method developed by Bielak and Christiano (1984)])
and Cremonini et al (1988) that expresses
the earthquake excitation
in the form of applied nodal forces acting on a strip of finite
elements located outside the region of interest. The interface
between the two regions is chosen as the circular arc II-II.
The basic idea is straightforward. One introduces a new variable w
which in the region interior to II-II coincides with the total
displacement field u and on the exterior is equal to the scattered
field 
. Since 
is a solution of
the wave equation over the entire halfspace, then w is also governed
by Eq. 1 in each subdomain. To ensure that the
solution for w is
equivalent to that of the original problem, continuity of the total
displacement u and the corresponding traction must be imposed across
II-II. When expressed in terms of w, these conditions introduce
the free-field displacement 
and the corresponding traction
on II-II explicitly into the formulation. It is these terms that upon
discretization give rise to the equivalent nodal forces on exterior
elements containing nodes on II-II.
Besides yielding the equivalent seismic forces, which are exact
within discretization error, another advantage of formulating the
problem in terms of the variable w is that
the motion exterior to II-II will be purely outgoing. The
absorbing boundary I-I need then transmit only outgoing waves.
The reduced problem defined over the bounded domain shown on Fig. 3a is solved with a 2D version of the parallel elastic finite element wave propagation code mentioned in the preceding section. The components of the system include a 2D and a 3D mesh generator, a mesh partitioner, a parceler, and a parallel code generator, as well as parallel numerical methods for applying the seismic forces, incorporating the absorbing boundaries, and solving the discretized wave propagation problem (Bao et al, 1996). The finite element mesh, shown in Fig. 3b, is generated by Triangle, a 2D quality mesh generator (Shewchuk, 1996). Triangle operates on geometric boundaries and interfaces defined by piecewise straight segments, as illustrated in Fig. 3a. The length of each segment is chosen according to the local material properties and maximum frequency of interest, so that Triangle produces a well-tailored unstructured mesh, as shown in Fig. 3b. Even though this problem is small and can be solved on a single processor, the simulations were performed on 64 processing elements of an Intel Paragon and a Cray T3D as a test of the parallel implementation. The corresponding partitioned mesh is presented in Fig. 3c. Each subdomain is mapped onto a processor of the parallel machine.
Quadratic six-node subparametric elements and standard Galerkin ideas are used for the spatial discretization of the governing wave equation over the triangular mesh in Fig. 3b. This leads to a system of ordinary differential equations of the form:
in which 
is the vector of time-dependent nodal displacements,
an overdot denotes differentiation with respect to time,
M, C, and K are the
global mass, damping, and stiffness matrices, and 
is the vector
of equivalent nodal forces. The stiffness and damping matrices
K and C contain terms representing the stiffness and
material damping in the soil as well as those arising from the
absorbing boundary condition. Damping in the soil is assumed to
be of the Rayleigh type, that is, within each finite element the
damping matrix is of the form:
in which 
and 
are arbitrary scalars,
is a reference frequency, and M 
and
K 
are the element mass and stiffness matrices. With this
choice of damping, the mass proportional and stiffness proportional
terms in Eq. 2 yield a damping ratio that is inversely
and directly proportional to frequency, respectively. Since the damping
ratio in soils is usually assumed to remain constant independently of the
rate of application of the load, Rayleigh damping can represent this
behavior only approximately. This approximation is achieved here by
evaluating 
and 
for each soil material such that
the squared difference between the target damping ratio within each
soil and the actual damping ratio specified over a prescribed
frequency range of interest [ 
, 
] is
minimized; that is
The range [ 
, 
] is selected so as to include
the dominant frequencies of the earthquake excitation and the dominant
resonant frequencies of the valley. The same approach can be used if
is frequency-dependent.
After assembly of the individual mass, damping, and stiffness matrices, and of the equivalent seismic forces, Eq. 2 is solved numerically by the unconditionally stable Newmark's trapezoidal step-by-step method. An iterative method, Jacobi preconditioned conjugate gradients, is used to solve the resulting algebraic equations at each time step, as this is more efficient than a direct solver on a parallel computer.
It is important to emphasize that while Eq. 2 provides an approximate solution for the total displacement field within the valley and the surrounding region due to an arbitrary incident SH-wave generated outside the computational domain, this equation strictly characterizes the response of a bounded domain. Since M and K are symmetric and positive definite, it then follows that the reduced system that governs the undamped free vibrations of the elastically supported valley, which is obtained from Eq. 2 by setting to zero the damping and forcing terms, possesses undamped natural frequencies and corresponding orthogonal modes. For the damped case these modes are not classical, in general, due both to the type of soil damping considered and to the contribution to C from the absorbing boundary condition. However, the damped system has complex conjugate eigenvalues and associated complex conjugate orthogonal eigenvectors (Foss, 1958). The eigenvalues represent damped natural frequencies and exponential attenuation. The response to a prescribed excitation, including the incident seismic wave, can be expressed as a linear combination of the individual modes. The extent to which each individual complex conjugate eigenvector pair of the valley model (mode shapes) contributes to the total response depends, of course, on the particular spatial and temporal distribution of the seismic excitation.
