ISSN 09655425, Computational Mathematics and Mathematical Physics, 2013, Vol. 53, No. 10, pp. 1523–1533. © Pleiades Publishing, Ltd., 2013. Original Russian Text © V.I. Golubev, I.B. Petrov, N.I. Khokhlov, 2013, published in Zhurnal Vychislitel’noi Matematiki i Matematicheskoi Fiziki, 2013, Vol. 53, No. 10, pp. 1709–1720.
Numerical Simulation of Seismic Activity by the GridCharacteristic Method V. I. Golubev, I. B. Petrov, and N. I. Khokhlov Moscow Institute of Physics and Technology, Institutskii per. 9, Dolgoprudnyi, Moscow oblast, 141700 Russia email:
[email protected] Received February 21, 2013
Abstract—Seismic activity in homogeneous and layered enclosing rock masses is studied. A numerical mechanicalmathematical model of a hypocenter is proposed that describes the whole range of elastic perturbations propagating from the hypocenter. Synthetic beachball plots computed for various fault plane orientations are compared with the analytical solution in the case of homogeneous rock. A detailed analysis of wave patterns and synthetic seismograms is performed to compare seismic activ ities in homogeneous and layered enclosing rock masses. The influence exerted by individual compo nents of a seismic perturbation on the stability of quarry walls is analyzed. The gridcharacteristic method is used on threedimensional parallelepipedal and curvilinear structured grids with boundary conditions set on the boundaries of the integration domain and with welldefined contact conditions specified in explicit form. DOI: 10.1134/S0965542513100060 Keywords: hyperbolic equations, gridcharacteristic method, mathematical modeling, heterogeneous media.
INTRODUCTION In the modern world, progressively more attention is given to the assessment and provision of popula tion safety. Active research is being conducted concerning imitation modeling methods [1, 2]. They pro vide a comprehensive analysis of risks and make it possible to identify the most significant of them, which need to be reduced. Obviously, any simulation is based on some approximation (model) of actual physical processes. Previously, models were mainly simple and the governing equations could be solved by analyt ical methods, but, with the emergence of highperformance computers, they have changed qualitatively. Problem formulations are becoming more sophisticated. For example, they include a description of the motion of individual atoms and molecules in order to simulate macroscopic phenomena, such as water flow in rivers [3] and landslides [4]. A process attracting the attention of the scientific community is natural seismicity, which gives rise to earthquakes. Over the last two years, more than 10 large earthquakes have occurred around the world (China, Turkey, Japan), which have caused a great amount of damage and cost thousands of human lives. Thus, the importance of the study of seismic activity aimed at its prediction and fast assessment of possible destruction is beyond any doubt. For this purpose, the globe is covered with a network of seismic stations that continuously record motions of the ground. The densest networks are deployed in Californium and Japan. As a rule, a seismometer measures the ground velocity or acceleration (or parts of their projections) as a function of time at the point of its location. In the case of an earthquake, threedimensional elastic strain waves propagate from its hypocenter, which usually lies at a depth of 5–30 km. Reaching the Earth’s surface, they cause destruction. At present, numerous methods are available that can determine the geographical coordinates and depth of a hypocenter from seismometer records (see [5]). Unfortunately, these methods have a shortcom ing: seismic wave propagation is mainly described in the ray approximation. The cause is that the fullwave simulation of dynamic processes in a geological medium is of high computational complexity. In this paper, we develop numerical methods and parallel algorithms that, in an acceptable time, compute elastic wave propagation through a heterogeneous medium from a hypocenter toward the Earth’s day surface, the formation of Rayleigh and Love surface waves [6], and their effect on ground structures [7]. 1523
1524
GOLUBEV et al.
North ⎡θ
u
λ
δ
Fig. 1. Fault slip model of an earthquake’s hypocenter.
1. MECHANICALMATHEMATICAL MODEL OF AN EARTHQUAKE’S HYPOCENTER One of the research goals addressed in this work is the development of a numerical mechanicalmath ematical model of a hypocenter in order to describe the initiation of seismic activity in rock. It should be noted that there is a variety of hypocenter models in geophysics (see, e.g., [8]). However, no one can be considered unquestionably the best as based on available experimental data. Moreover, a more compli cated model contains more free parameters and, accordingly, the procedure for determining them from ground measurements becomes more complicated. In this paper, we use the simplest but physically cor rect model called a fault slip. This model relies on the geophysical idea that an earthquake arises from the gradual accumulation of elastic stresses in the Earth’s crust, which consists of separate rock blocks. Since the faults between the blocks are filled with breccia (rock composed of broken fragments larger than 1 cm in size cemented together) and the normal compressing regional stresses are fairly high, the effective frictional force between the fault edges is large. However, when a certain threshold is exceeded, boundary blocks begin to slide along the fault and, accordingly, the built up strain releases abruptly. Note that elastic, plastic, and brittleductile processes can be observed near the fault. Away from the site of the perturbation initiation, important are only elastic processes, which are represented by elastic waves propagating through the rock. Consider the proposed model of a hypocenter in more detail. Figure 1 shows the contact between two blocks (the contact is a fault). The orientation of a plane in threedimensional space is uniquely deter mined by two inclination angles. One of them is known as the strike angle (θ) and is defined as the angle made by the line of intersection of the fault plane with a horizontal plane and the north direction. The strike angle ranges from 0° to 360° and is measured in the clockwise direction when viewed from above. The second angle is called the dip (δ). It is defined as the angle between the fault plane and a horizontal plane and ranges from 0° to 90°. The model assumes that the blocks do not move in the normal direction to the fault plane. Thus, the velocity of the medium lies in the fault plane. The orientation of a vector lying in a plane can be specified by a single parameter. As such a parameter, we use the rake angle (λ) defined as the angle between the direction along the line of intersection of the fault plane with a horizontal plane (strike direction) and the direction of slipping blocks (slip direction). The rake angle ranges from –180° to +180°. At the rake angle equal to 0, the observer standing on the ground surface and looking in the strike direction sees that the right edge of the fault recedes in this direction. For the rake angle equal to plus or minus 180°, the observer sees that the right edge of the fault moves in the opposite direction. The velocity of the right fault edge has a positive vertical component for positive rake angles and a negative vertical component for negative rake angles. A characteristic of an earthquake is its magnitude (i.e., a quantity characterizing the energy released in an earthquake in the form of seismic waves), which determines the displacements of the ground surface. In the model under consideration, a nonzero medium velocity that is constant in magnitude and direction (though having an inversion discontinuity on the fault) is specified at the initial time in a bounded domain near the contact of blocks (see Fig. 2). In view of what was said above, we still have an undetermined parameter, namely, the slip velocity, which can be estimated by directly simulating the process and com paring the computed values of ground displacements with experimental data. Moreover, the size of the domain in which we specify a nonzero velocity is picked so that the map of contour lines of displacement velocities (which can be constructed from experimental data) is reproduced as accurately as possible. Of special interest is a common method for representing hypocenter information in the form of a beachball plot. Let us describe in detail the physics behind this image. A hypocenter is enclosed by a sphere of some radius such that the slips over its surface are initially zero. As the process develops, the displace ment vector, more precisely, its projection onto the position vector directed from center of the sphere is recorded at each point on the lower half of the spherical surface. As soon as this projection becomes non COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 53
No. 10
2013
NUMERICAL SIMULATION OF SEISMIC ACTIVITY
1525
−v
v
Fig. 2. Spatially constant initial velocity (with an inversion discontinuity).
zero for the first time, the given point is labeled by a plus if the projection is positive and by a minus if the projection is negative. After the perturbation has propagated beyond the observed sphere, there are three classes of points (labeled by a plus, minus, and unlabeled). The lower halfsphere is projected stereograph ically onto a plane and the points corresponding to different signs are marked with different colors. As a result, a beachball plot is produced for the given earthquake. Numerous seismometers are deployed around the globe. They measure motions of the ground and are joined in seismic networks. Special attention is given to the recording of early motion, which shows the direction of the first nonzero motion of the ground (from the arrival of a Pwave): away from the hypo center or toward it. A special procedure is used that reconstructs a beachball plot from the wave paths of a single event recorded at several seismic stations. Note that, in numerical simulation, the medium velocity is known at each point of the media at successive times. Therefore, a beachball plot can be obtained directly from its definition. 2. MATHEMATICAL MODEL OF ROCK AND A NUMERICAL METHOD FOR SOLVING THE PROBLEM The state of an infinitesimal volume of a linear elastic medium is described by the linear dynamic elas ticity equations. Consider the nonstationary elasticity equations in the case of three variables in some orthonormal of coordinate system (x1, x2, x3):
ρ ⋅ v i = ∇ j ⋅ σij , σ ij = qijkl ε kl + Fij .
(1)
Here, ρ is the density of the medium, v i are the components the displacement velocity, σ ij and ε kl are the respective components of the Cauchy stress tensor and the strain tensor, ∇ j is the covariant derivative with respect to the j th coordinate, and Fij is an additional righthand side. The form of the fourthorder tensor components qijkl is determined by the rheology of the medium. In the linear elastic case, they are given by
qijkl = λδ ij δ kl + μ(δ ik δ jl + δ il δ jk ) , which is an generalization of Hooke’s law. Here, λ and µ are the Lamé parameters and δ ij is the Kronecker delta. The first vector equation in system (1) represents three equations of motion, while the second equation represents six rheological relations. The vector of unknown functions consists of nine components and has the form
u = {v1,v 2,v 3, σ11, σ12, σ13, σ22, σ23, σ33}T. Note that the solid mechanics equations (1) can be written in the matrix form (see [9])
∂u = ∂t
3
∑A j =1
j
∂u , ∂x j
(2)
where A j is a 9 by 9 matrix. COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 53
No. 10
2013
1526
GOLUBEV et al.
The gridcharacteristic method is widely applied to numerical simulation in solid mechanics [10, 11]. First, a splitting method with respect to spatial coordinates is used to obtain three onedimensional sys tems of equations ∂u = A ∂u , j = 1,2,3. (3) j ∂t ∂x j Each of these systems is hyperbolic and has a complete set of eigenvectors with real eigenvalues. Each of the systems can be rewritten as ∂u = Ω −1Λ Ω ∂u , j j j ∂t ∂x j where Ω j is the eigenvector matrix and Λ j is a diagonal eigenvalue matrix. For all coordinates, the matrix Λ has the form Λ = diag{c p, −c p, cs , −cs , cs , −cs ,0,0,0}, where c p = (λ + 2μ)/ ρ is the longitudinal speed of sound in the medium and c s = μ/ ρ is the transverse speed of sound. After changing to the variable v = Ω u , each system in (3) splits into nine independent scalar advection equations (in what follows, the index j is omitted wherever possible) ∂v + Λ ∂v = 0. ∂t ∂x The onedimensional advection equations are solved by applying the method of characteristics or usual finitedifference schemes. After all the components v have been transferred, the solution is reconstructed as u
n +1
−1
=Ω v
n +1
.
In this work, we used the secondorder accurate TVD schemes [12]
umn+1 = umn − σ( f m+1/2 − f m−1/2 ), where f m+1/2 and f m−1/2 are antidiffusion fluxes given by
f m+1/2 = umn − 0.5φ(rm )(1 − σ)(umn +1 − umn ), f m−1/2 = umn −1 − 0.5φ(rm−1)(1 − σ)(umn − umn −1). The function φ(rm ) is called a limiter. The parameter rm is calculated using the formula u − um−1 . rm = m um+1 − um The program used involves 15 different limiters (see [13]). Most frequently, we used the superbee [14] and MC (monotonized central) limiters [15], which have proved the best: φ sb(r ) = max[0, min(2r,1), min(r,2)], φ MC (r ) = max[0, min(2r,0.5(1 + r ),2)]. Additionally, we used gridcharacteristic monotone difference schemes (see [16, 17]). Second to fourthorder accurate schemes were implemented in the program. Most computations were based on the fourthorder accurate scheme. As applied to the onedimensional linear advection equation, this scheme is given by
umn+1 = umn − σ(Δ1 − σ(Δ 2 − σ(Δ 3 − σΔ 4 )), Δ1 = 1 −2umn +2 + 16umn +1 − 16umn −1 + 2umn −2 , 24 1 Δ2 = −umn +2 + 16umn +1 − 30umn + 16umn −1 − umn −2 , 24 Δ 3 = 1 2umn +2 − 4umn +1 + 4umn − 2umn −2 , 24 Δ 4 = 1 umn +2 − 4umn +1 + 6umn − 4umn −1 + umn −2 . 24
(
)
(
)
(
(
)
)
COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 53
No. 10
2013
NUMERICAL SIMULATION OF SEISMIC ACTIVITY (a)
1527
(b)
Fig. 3. Comparison of (a) analytical and (b) numerical beachball plots for homogeneous enclosing rock. Strike = 70°, Dip = 90°, and Slip = 0°.
Moreover, we used the gridcharacteristic monotonicity criterion from [16]. For positive λ the crite rion has the form
{
}
{
}
{ {
} } {
{ {
} }
min umn , umn −1 ≤ umn+1 ≤ max umn , umn −1 . For negative λ , the criterion is symmetric. In the simplest implementation, if this criterion is violated, the solution is corrected as
⎧max umn , umn −1 , umn+1 > max umn , umn −1 , ⎪⎪ n+1 n n n+1 n n um = ⎨min um, um−1 , um < min um, um−1 , ⎪ n+1 n n n+1 n n ⎪⎩um , min um, um−1 ≤ um ≤ max um, um−1 . This limiter remains fourthorder accurate in domains where the solution is fairly smooth (the charac teristic criterion holds). In the case of high solution gradients, the order of the scheme is reduced to the third (see [16]). If the computational domain is not a parallelepiped or a combination of parallelepipeds, the possibility of using curvilinear structured grids is provisioned in the program. Curvilinear grids are constructed by mapping the physical integration domain to a parallelepiped with a uniform grid. In this case, splitting with respect to coordinates is performed in the transformed rather than physical space. The original system (2) becomes ∂u = ∂t
3
A j ∂u , ∂ξ j j =1
∑
}
A j =
3
∂ξ j
∑ ∂x i =1
{
}
Ai .
i
For each splitting step in the direction ξ j , the matrix Λ has the form
Λ = diag{ c1l j , −c1l j , c2l j , −c2l j , c2l j , −c2l j ,0,0,0}. Here, l j = v j = (v1j )2 + (v 2j )2 + (v 3j )2 and v j = ∇ ⋅ ξ j . The transformation matrices Ω and Ω −1 are also modified. To determine the coefficients v j , the inverse of the Jacobian matrix has to be computed in each cell. The coefficients can be computed analytically if the transformation is given by ξ = ξ(x) or numerically by using the secondorder accurate formulas (v ij )m = ((ξ j )m+1 – (ξ j )m−1)/ Δxi . The structure of the computational domain is defined by a set of parallelepipedal blocks with edges par
∪
m
allel to the Cartesian axes: Ω = Ω . Each of the blocks Ω k is characterized by its own set of medium k =1 k parameters and boundary and contact conditions. In each Ω k , rectangular grids are constructed with the COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 53
No. 10
2013
1528
GOLUBEV et al.
(a)
(b)
Fig. 4. Comparison of (a) analytical and (b) numerical beachball plots for homogeneous enclosing rock. Strike = 90°, Dip = 20°, and Slip = 0°.
(a)
(b)
Fig. 5. Comparison of (a) analytical and (b) numerical beachball plots for homogeneous enclosing rock. Strike = 90°, Dip = 50°, and Slip = 0°.
step hk = {hk1, hk 2, hk3}, which is constant inside that block. If there is a contact boundary in neighboring blocks, the corresponding grids have to be consistent. By the consistency of grids, we mean that their mesh sizes are identical in the plane of the contact and the coordinates of the nodes in the adjacent grids coin cide. The unknown functions on the grid boundaries in neighboring blocks were computed using two types of boundary conditions: noslip and slip conditions. 3. NUMERICAL SOLUTION OF THE PROBLEM 3.1. Comparison of Numerical Results with the Analytical Solution in the Case of Homogeneous Enclosing Rock In this paper, we performed the direct simulation of dynamic processes occurring in a threedimen sional rock mass with elastic waves propagating from a hypocenter modeled by a fault slip. The enclosing rock was assumed to be elastic, homogeneous, and isotropic. The characteristics were specified as follows: COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 53
No. 10
2013
NUMERICAL SIMULATION OF SEISMIC ACTIVITY
(a)
1529
(b)
Fig. 6. Comparison of (a) analytical and (b) numerical beachball plots for homogeneous enclosing rock. Strike = 90°, Dip = 90°, and Slip = 75°.
Fig. 7. Speed distribution for an earthquake in a homogeneous elastic medium. Strike = 90°, Dip = 0°, and Slip = 0°.
the rock density was 2000 kg/m3, the longitudinal wave velocity was 2000 m/s, and the transverse wave velocity was 1300 m/s. During the numerical simulation, a beachball plot was computed by tracing values in the lower halfsphere enclosing the hypocenter and lying completely in the elastic medium. Note that, in the given approximation, beachball plots can also be constructed using analytical methods (see [18]). Examples comparing an analytical solution with numerical results are given in Figs. 3–6. The nonunifor mity of points lying on the computed beachball plot is explained by the coarse grid used and by the features of the subsequent construction of a stereographic projection. The illustrations presented show that the numerical results agree with the analytical solution to high accuracy. It should be noted that an analytical solution describing the stress and velocity distributions in the medium as functions of coordinates and time cannot be obtained even in the simplest case of a homo geneous medium. The computed speed distribution in the rock at a fixed time is depicted in Fig. 7. A detailed analysis [7] shows that a set of elastic waves, namely, two longitudinal and two transverse waves propagate from the hypocenter toward the ground surface. Since the velocity of the former is considerably COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 53
No. 10
2013
1530
GOLUBEV et al.
Characteristics of the layered medium Layer
Thickness, m
Density, kg/m3
Longitudinal wave velocity, m/s
Transverse wave velocity, m/s
1 2 3 4
300 400 500 1600
2500 2500 2500 2500
4190 4650 5250 5850
2793 3100 3500 3900
higher than that of the latter, they gradually separate farther apart and the longitudinal waves are the first to reach the surface and initiate early motions on the wave paths. However, the amplitudes of the trans verse waves are much larger and it is them that give rise to Rayleigh and Love surface waves, which cause major destructions. (a)
(b)
Fig. 8. Comparison of wave patterns at successive times for (a) homogeneous and (b) layered enclosing rock. Strike = 90°, Dip = 0°, and Slip = 0°. COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 53
No. 10
2013
NUMERICAL SIMULATION OF SEISMIC ACTIVITY (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
1531
(b) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Fig. 9. Comparison of seismograms for (a) homogeneous and (b) layered enclosing rock. Strike = 90°, Dip = 0°, and Slip = 0°.
Fig. 10. Description of the geometry of the computational domain with the use of curvilinear structured grids.
3.2. Dynamics of Wave Processes in Rock under Earthquake: Comparison of Homogeneous and Layered Rock An advantage of computer simulation of seismicity is that the characteristics of wave processes can be analyzed in detail depending on the geometry and parameters of the enclosing rock mass. Since the veloc ities and stresses at each point of an elastic body are known in the computation, it is convenient to visualize results by using threedimensional wave patterns in which the speed of the medium at each point is shown in color. However, in field experiments, measurements are available only on the ground surface at sites where seismometers are situated. As a rule, these instruments measure the horizontal or vertical projection of the velocity or acceleration and only rarely their complete vectors. Wave paths measured by individual seis mometers are used to construct generalized graphs for several seismometers, which are known as seismograms. The vertical axis represents the time from the beginning of the records, which increases downward. In this paper, we numerically simulated an elastic perturbation propagating from a hypocenter through a homogeneous elastic medium and a layered elastic medium. The parameters of the homogeneous medium were specified as follows: the density was 2500 kg/m3, the longitudinal wave velocity was set to 5850 m/s, and the transverse wave velocity was 3900 m/s. The layered medium consisted of four layers with different elastic characteristics. Their values are given in the table. The results of the numerical simulation are displayed in Fig. 8 (vertical cross section) in the form of a sequence of wave patterns. The basic qualitative difference in the latter case is the occurrence of numerous multiply reflected waves because of the differences in the elastic characteristics of the layers. Moreover, some of the energy COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 53
No. 10
2013
1532
GOLUBEV et al.
Fig. 11. Wave patterns of the seismic effect on a quarry at successive times.
is carried by reflected waves back deep into the media. As a result, the waves reaching the surface become weaker. The numerical results were also used to construct synthetic seismograms for both types of media. They clearly show the difference in the arrival times of the corresponding types of waves, which is explained by the vertical variations in the elastic characteristics in the layered medium. 3.3. Effect of Individual Perturbation Components on the Stability of Quarry Walls A feature of the gridcharacteristic method used is that we can specify an arbitrary threedimensional geometry of the computational domain and boundary conditions can be formulated in explicit form. In this paper, a threedimensional curvilinear grid describing the geometry of a conical quarry was con structed (Fig. 10). It consists of 13 grids, each being a structured one. The noslip condition is used at their contacts, which ensures no artifacts in the form of smallamplitude reflected waves. The large number of grids ensures that the ratio of the lengths of the longest to shortest grid element edges is small, which improves the accuracy of the computations. The material of each grid has the following elastic character istics: the density is 2000 kg/m3, the longitudinal wave velocity is 2000 m/s, and the transverse wave veloc ity is 1300 m/s. As an initial perturbation, we used a longitudinal wave (one of the components of the seis mic signal from the hypocenter) propagating from deep in the rock to its surface. Figure 11 shows wave patterns at successive times. CONCLUSIONS A numerical mechanicalmathematical model of a hypocenter was developed, which was used in the computer simulation of seismic activity. The problem of initiating an earthquake in a homogeneous medium was solved numerically by applying the gridcharacteristic method. Synthetic and analytical beachball plots were compared for various fault plane orientations. A complete spectrum of perturbations propagating from a hypocenter to the day ground surface in the case of homogeneous and layered media was compared in fullwave simulation. An analysis of synthetic seismograms showed that the transmitted waves become weaker due to the multiple reflections from the boundaries of the layers having different acoustic rigidity. COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 53
No. 10
2013
NUMERICAL SIMULATION OF SEISMIC ACTIVITY
1533
The numerical method used on parallelepipedal and curvilinear structured grids makes it possible to simulate wave processes in domains of complex geometry with boundary and contact conditions formu lated in explicit form. The incidence of a longitudinal wave on a quarry was simulated with computing the complete velocity of the medium and the stress tensor. These results can be used to estimate the stability of the quarry walls. ACKNOWLEDGMENTS This work was supported by the Ministry of Education and Science of the Russian Federation (contract no. 14.132.21.1809) and by the Russian Foundation for Basic Research (project no. 120731028). REFERENCES 1. V. V. Gaskin and V. I. Sobolev, “Simulation of seismic processes in extended building strictures,” Sovrem. Tekh nol. Sist. Anal. Model., No. 5, 26–33 (2005). 2. T. V. Mikheeva, “Overview of simulation software tools for the study of operation and control mechanisms of manufacturing systems,” Izv. Altaisk. Gos. Univ., No. 1, 87–90 (2009). 3. A. N. Semchukov, Candidate’s Dissertation in Mathematics and Physics (2004). 4. O. E. Khvostova and A. A. Kurkin, “Mathematical modeling of subsidence tsunami by the particle method,” Vestn. Belgorod. Gos. Tekhnol. Univ., No. 4, 96–100 (2009). 5. Z. A. Nazarova, “Accuracy and stability testing of software programs (DIMAS, ARC) for computing the coor dinates of an earthquake’s hypocenter in the case of the Kamchatka seismic network,” in Abstracts of the Scien tific Conference on Complex Geophysical Monitoring of the Russian Far East, November 11–18, 2007 (Petropav lovskKamchatskii), pp. 191–195. 6. I. B. Petrov and N. I. Khokhlov, “Modeling of seismic phenomena by the gridcharacteristic method,” Trudy Mosk. Fiz.Tekh. Inst. 3 (3), 159–167 (2011). 7. V. I. Golubev, I. E. Kvasov, and I. B. Petrov, “Influence of natural disasters on ground facilities,” Math. Model. Comput. Simul. 4 (2), 129–134 (2012). 8. L. M. Balakina, A. V. Vvedenskaya, N. V. Golubeva, L. A. Misharina, and E. I. Shirokova, Elastic Stress Field of the Earth and Focal Mechanism of Earthquakes (Nauka, Moscow, 1972) [in Russian]. 9. I. B. Petrov and A. S. Kholodov, “Numerical study of some dynamic problems of the mechanics of a deformable rigid body by the meshcharacteristic method,” USSR Comput. Math. Math. Phys. 24 (3), 61–73 (1984). 10. K. M. Magomedov and A. S. Kholodov, GridCharacteristic Numerical Methods (Nauka, Moscow, 1988) [in Russian]. 11. K. M. Magomedov and A. S. Kholodov, “The construction of difference schemes for hyperbolic equations based on characteristic ratios,” Comput. Math. Math. Phys. 9 (2), 158–176 (1969). 12. A. Harten, “High resolution schemes for hyperbolic conservation laws,” J. Comput. Phys. 135 (2), 260–278 (1997). 13. I. B. Petrov and N. I. Khokhlov, “Comparison of TVD limiters for the numerical solution of solid mechanics equations produced by the gridcharacteristic method,” in Mathematical Models and Control Problems: Collec tion of Scientific Works of Moscow Institute of Physics and Technology (Moscow, 2011), pp. 104–111 [in Russian]. 14. P. L. Roe, “Characteristicbased schemes for the Euler equations,” Annu. Rev. Fluid Mech. 18, 337–365 (1986). 15. B. Van Leer, “Towards the ultimate conservative difference scheme III: Upstreamcentered finitedifference schemes for ideal compressible flow,” J. Comp. Phys. 23, 263–275 (1977). 16. A. S. Kholodov and Ya. A. Kholodov, “Monotonicity criteria for difference schemes designed for hyperbolic equations,” Comput. Math. Math. Phys. 46, 1560–1588 (2006). 17. A. S. Kholodov, “The construction of difference schemes of increased order of accuracy for equations of hyper bolic type,” USSR Comput. Math. Math. Phys. 20 (6), 234–253 (1980). 18. F. Shahzad, “Software Development for Fault Plane Solution and Isoseismal Map,” M. Sc. Thesis QuaidiAzam Univ. Islamabad (2006).
Translated by I. Ruzanova COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 53
No. 10
2013