Page 284 - Numerical Analysis and Modelling in Geomechanics

ENRICO PRIOLO 265
Issues specific to the implementation
The order of the polynomials can be set arbitrarily. Usually, orders equal to six
or eight are chosen, which, for this kind of formulation, have been shown to be
the best compromise between accuracy, computational efficiency, and memory
requirement (Padovani et al., 1994). The method needs a low number of grid
points per wavelength (G), and the accuracy does not degrade even for very long
propagation times. For example, values of G=4.5 and G=5.2 are typically used
with polynomial orders N=8 and N=6, respectively. Material inhomogeneity is
modelled simply by defining different material parameters for adjacent elements.
Material characteristics do not vary inside each element. The free-surface
boundary condition is obtained simply by imposing no constraints at boundary
nodes (Hughes, 1987). The inner reflection of the outgoing wavefield at external
boundaries is eliminated through absorbing strips where the wavefield is
smoothly attenuated (Cerjan et al., 1985). Attenuation Q is introduced in a
simplified form (Graves, 1996) with a single value for both P and S waves. In
this approximation, Q is set for a reference frequency, and it depends linearly on
the frequency. It can be shown that this form approximates the constant Q fairly
well for frequencies near the reference frequency. The value of Q is set on the
shear wave velocities, which are the dominant waves in the wavefield generated
by an earthquake.
The method presented here is 2-D. This means that both model structure and
source extend infinitely in the direction perpendicular to the vertical model plane.
The earthquake source is a point shear dislocation with no torsion, and it is
equivalent to a point double-couple. It is introduced through a field of external
forces. The 2-D hypothesis implies that the point source is actually a 3-D line
source orthogonal to the model plane. The time history used to simulate
earthquakes, i.e. the slip velocity function, is a unit Ohnaka impulse (Herrmann,
1996) of the form:

This function is defined by only one parameter α, which relates to the corner
2
frequency f by α=2πf . The Fourier amplitude spectrum has a typical ω decay.
c
c
The source is scaled by setting the slip magnitude and duration through some
empirical relations (Wells and Coppersmith, 1994; Somerville et al., 1999).
The use of irregular grids and isoparametric high-order elements makes it
possible to adjust the element size to the minimum wave velocity locally, and
follow all structure interfaces exactly (Figure 9.1). Technically, the mesh of the
computational model is generated by two steps. First, the model is decomposed
into several sub-regions, which are meshed one at a time. Then the complete
mesh is obtained by sticking together all the meshed sub-regions. Irregular
regions are first meshed (Priolo, 2001) into triangles using an element size that is
double the final one; then the triangles are joined two by two; uncoupled