CC BY
7
DOI: 10.17516/1997-1397-2020-13-5-644-654 УДК 519.63
Finite Difference Schemes for Modelling the Propagation
of Axisymmetric Elastic Longitudinal Waves
Vladimir M. Sadovskii* Oxana V. Sadovskaya^ Evgenii A. Efimov*
Institute of Computational Modeling SB RAS Krasnoyarsk, Russian Federation
Received 16.06.2020, received in revised form 11.07.2020, accepted 16.08.2020 Abstract. An efficient finite difference shock-capturing scheme for the solution of direct seismic problems is constructed. Problem formulation is based on equations of the dynamics of elastic medium with axial symmetry. When implementating the scheme on multiprocessor computing systems, the two-cyclic splitting method with respect to spatial variables is used. One-dimensional systems of equations that arise in the context of splitting procedure are represented as subsystems for longitudinal, transverse and torsional waves. The case of longitudinal waves is considered in this paper. The results of simulations with the use of explicit grid-characteristic schemes and implicit schemes of the "predictor-corrector" type with controllable dissipation of energy are compared with exact solutions that describe propagation of monochromatic waves.
Keywords: elastic medium, cylindrical waves, splitting method, finite difference scheme, monotonicity, dissipativity, parallel computing.
Citation: V. M. Sadovskii, O.V. Sadovskaya, E. A. Efimov, Finite Difference Schemes for Modelling the Propagation of Axisymmetric Elastic Longitudinal Waves, J. Sib. Fed. Univ. Math. Phys., 2020, 13(5), 644-654. DOI: 10.17516/1997-1397-2020-13-5-644-654.
Introduction
The key problem in simulation of axially symmetric wave propagation consists in selection of appropriate approximation of lowest terms in equations of dynamic elasticity written in cylindrical system of coordinates. These terms cause the degeneracy of the equations at the axis of symmetry. Our aim is to find such appropriate approximations while remaining within the framework of conservative Godunov's scheme which is widely used in numerical solution of two-and three-dimensional problems [1]. The scheme can be modified to simulate wave propagation in granular and porous materials with different resistances to compression and tension [2-4], wave propagation and fracturing in blocky media [5-7] and other nonlinear processes.
Many methods were developed to solve axially symmetric equations of dynamic elasticity. The method of characteristics was used for the analysis of one-dimensional cylindrical wave [8,9]. To solve two- and three-dimensional equations finite difference schemes based on the method of characteristics were proposed [10,11]. These schemes allow one to compute the discontinuities of velocities and stresses. Such methods were applied for the analysis of wave processes in linear elastic, viscoelastic and elastic-plastic media [12,13]. Numerical analysis of the dynamics of plates and shells of revolution was implemented with the methods based on the axially symmetric equations [14-16].
* sadov@icm.krasn.ru https://orcid.org/0000-0001-9695-0032 to_sadov@icm.krasn.ru https://orcid.org/0000-0001-7222-8088 teugene6467@mail.ru © Siberian Federal University. All rights reserved
1. Equations of axially symmetric motion
Equations of the dynamic theory of elasticity with axially symmetric fields of stresses and velocities are written as
•) d (r&rz) dvv d (r arv) d (r avz)
pr dvr ~dt _ d(r a, dr
pr dvz dt _ d(r a, dr
1 dap v d
E dt E dt
1 darp dvp
- dt dr
I y * " / j r f _ ^ f ' | ^ f ' | j
dz p ' dt dr dz rp '
) d(raz) 1 dar v d , n dvr
~ + ~e ~st — Edt[(Jp + _ (1)
\ vr 1 daz v d , N dvz
Tz + _ V Ed- — Em{ar + ^ _ ~dz>
vp 1 darz dvz dvr 1 dapz dvp
r 1 ¡i dt dr dz 1 ¡i dt dz
Here E = 2 p (1 + v) is Young's modulus, p is the shear modulus, v is the Poisson ratio; the axes r and z of the cylindrical coordinate system are directed along the radius and the axis of symmetry, respectively. This form of equations is convenient for obtaining the equation of energy balance. Let us multiply the equations of motions by vr, vv, vz and the constitutive equations by rar, rav, raz, rarv, rarz, r avz. Then the left-hand and right-hand sides of the obtained equations are summed up. The result is
d ( vl + vp + v2
dt\pr-2-+
d f \ d
(2)
dr^vr ar + r vp Jrp + rvz Jrz j + dz^rvr Jrz + r vp Jpz + rvz Jz^ ,
where W is the elastic potential which is a quadratic form with respect to stresses: W _ (j2 + ap + j2 + 2 alp + 2 a^ + 2 Jpz - 1+v J + aP + Jz №
System (1) is a system of hyperbolic partial differential equations. It can be splitted into two independent subsystems. The first subsystem (equations 1,3-6 and 8) describes motion in plane of symmetry, and the second subsystem (equations 2, 7 and 9) describes torsional motion. Motion in rz-plane is represented by superposition of longitudinal and transverse waves with velocities
Lp
2 - 1 — v fp
Cs _ 4 -
p 1 — 2 v v p
respectively. The torsional wave has velocity cs.
For seismic analysis, system of equations (1) is solved using the method of two-cyclic splitting with respect to spatial variables with solution of one-dimensional problems in parallel mode at different stages. Contrary to the conventional splitting, the method of two-cyclic splitting maintains the second-order accuracy if second-order finite difference schemes are used to solve one-dimensional systems [17]. Numerical implementation of splitting in z direction presents no difficulties. We have the system with constant coefficients after all derivatives with respect to r and terms containing r are canceled. Then obtained system is splitted into subsystems of plane P-and S-waves. These subsystems can be solved using Godunov's scheme [1] or grid-characteristic finite difference scheme with limiting reconstruction of Riemann invariants [18].
One-dimensional system of equations in the direction of radial axis r is splitted into three subsystems of longitudinal, transverse and torsional waves. The system contains terms without derivatives. Thus direct application of standard finite difference schemes for plane problem of
the elasticity theory may lead to unwanted effects such as asymptotic instability, accumulation of rounding errors in simulations with many time steps or imbalance in momentum and energy. In consequence of these features the correctness of numerical results may be doubtful.
A common approach to suppress such effects is to use fully conservative finite difference schemes [19, 20] in combination with the method of artificial viscosity [21] that smoothes off oscillations of numerical solution in calculation of discontinuities due to artificial dissipation of energy. In the case of equations of the dynamic theory of elasticity, an approach to construct schemes with controllable dissipation of energy was developed [22,23]. This approach is applied to the subsystem of equations for one-dimensional cylindrical longitudinal waves.
2. Finite difference schemes for propagation of cylindrical longitudinal wave
The equations for longitudinal waves can be written in an equivalent form in terms of the Lame parameters A = 2jv/(1 — 2v) and j as follows
dvr d (r ar) do~r \dvr . vr
pr-dt =~dT — a*, Ht = (A + 2j) + Ar,
dav = ^ + (a + 2 j) — dOz = a(— + —^
dt dr r ' dt \ dr r J
Equation (2) of the energy balance for this subsystem takes the form
d (prVl + rW ^ = d(rvr ar)
dt\ 2 J dr
Integrating equations (3) over the rectangular space-time grid leads to discrete equations of the "corrector" step
0 vr — vr r+a+ — r pr 0 ar — Or M , 0 , v+ — v- , , v°r
pr =-h--av, —— = (A + 2 + Aro-
ov—oV = Av_^r + (a + 2 j) 4 , O. = A + A 4 .
t h r0 t h r0
(4)
In what follows, the values with circumflex belong to the upper time layer, and the values without circumflex belong to the lower time layer. The values with superscripts "±" belong to the right-hand and left-hand boundaries of a cell, respectively and r0 = (r+ + rr)/2. The input velocity v° and stress , alongside with a± and v±, are determined at the "predictor" step.
Multiplying equations (4) respectively by (vr + vr)/2, r0(ar + ar )/2, r0(aV + aV)/2 and r°(vz + az )/2, the difference analogue of the energy balance equation (2) for longitudinal waves is obtained. It has the form
0 v2 — v2 0 WW — W r+v+a+ — r~vrarr ^
pr0 -!——r + r0-= -r—!—-— D ,
2 t t h
D = r+a+ — rro- fv+ + v- vr + vA + v+ — v- /r+a+ + r-a- r0 Or + a/, +
22
0 vr + vr 0 aV + O V + 2 Vr 2 ' The idea to control the dissipation of energy consists in setting expression for D explicitly in the form of a positive definite quadratic form. This form can be identically equal to zero. Then we have a dissipation-free (totally conservative) scheme.
Let the quadratic form be D = 7 (v+ — vr ) /h? with a free parameter 7 = O(h) > 0, with an assumption that
0 _ vr + vr _ v+ + v r 0 _ av + av _ + V-Vr _ 2 _ 2 ' _ 2 _ 2 r+a+ + r-a- 0 &r + v+
(5)
' r T ' <-> r 0 u r T Vr
--r - _ Y -
2 2 Y
In this case, artificial dissipation of energy in the scheme is nonnegative. This automatically ensures stability of calculations. Further, it decreases with refinement of the grid and only depends on the rate of deformation of the medium. When 7 = 0 the scheme is totally conservative, and the law of energy conservation is preserved on discrete level. However, in practice this is inapplicable in calculations of discontinuities and solutions with high gradients because it results in nonmonotonic solutions.
Considering equations (4), the closure equations in the scheme with controllable energy dissipation lead to the system
r+a+ — r-a- = aj-i/2 v+ + bj-1/2v- + fj-1/2 ,
r+a+ + r-a- = Cj-i/2 v+ + dj-1/2 v- + gj-1/2 ,
where coefficients a,j_1/2, bj-1/2, Cj-1/2, dj-1/2, fj-1/2 and gj-1/2 depend on the cell number j = 1,2,... ,n (fractional indices belong to the centers of cells). They are calculated as
phr0 . rh t phr0 rh r
aj-i/2 = + (A + 2N4-^0 + A 2, bj-i/2 = + (A + 2— A 2 ,
„0 r.
t t r 2 y t t r 2 y
cj-1/2 _ a 2 + (A + 2 p) + h, dj-i/2 _ a 2 - (a + 2 - h"
„0 pvr _ O „,0,
fj-1/2 _ hav - 2 hr- -, 9j-1/2
t
Hence
2 r+a+ _ (aj-1/2 + Cj-1/2) v+ + (bj-1/2 + dj-1/2) v- + fj-1/2 + gj-1/2 , 2 r-V- _ icj-1/2 - aj-1/2) v+ + {ß'j-1/2 - bj-1/2) v- + gj-1/2 - f j-1/2 ■
(6)
Equating these expressions and changing index j, we obtain system of linear equations with three-diagonal matrix to determine velocities v+ = vj and v- = vj-1 at the cell boundaries:
Aj vj+1 + Cj vj + Bj vj-1 = Fj (7)
with coefficients
Aj = Cj+1/2 — aj+1/2 , Cj = dj+1/2 — bj + 1/2 — aj-1/2 — Cj-1/2 ,
Bj = — bj-1/2 — dj-1/2 , Fj = fj + 1/2 + fj-1/2 — gj+1/2 + gj-1/2 ■
0
This system is supplemented with the boundary condition v0 =0 at the axis of symmetry in
n r
the first cell of the grid and with the condition vn _ v in the last cell of the grid if the particle velocity v is set at r _ R or with the condition
(an-1/2 + cn-1/2) < + (bn-1/2 + dn-1/2) vn 1 + fn-1/2 + 9n-1/2 = 2
which follows from (6) if the external stress a is set at the boundary. In either case, the system of equations with boundary conditions is solved with the use of the three-point sweep method.
In this manner, the algorithm of transition to the next time level starts with computation of v± and a± using equations (6), (7) at the "predictor" stage. Then, v° and are determined from (5). The final computations of vr, ar, aV and az are performed using equations (4) of the "corrector" stage.
For comparison, we consider three versions of explicit finite difference schemes based on the solution of the Riemann problem. The schemes are constructed by means of approximation of the first equation in (3) written in equivalent (nonconservative) form
dvr dar
P dt dr + r The "predictor-corrector" scheme with explicit approximation of lowest terms
Vr - Vr a+ - a- ar - av ¿r - ,, , 0 v+ - v- Vr
P-= ---+-^— , -= (A + 2 p)---+ A —
t h r0 t h r0
= + (a + 2p) V0 , *±-£± = A f V^V- +
t h ru t \ h r0
+ a+ ar -, a- , ¿r
V+--= Vr--, V- +--= Vz +--,
pcp pcp pcp pcp
is applicable if the Courant number Kp = cp T/h is in the range from 0 to 0.8. When the Courant number exceeds 0.8 the parasitic oscillations appear in the vicinity of the axis of symmetry. The amplitude of oscillations unlimitedly increases with increasing Kp from 0.9 to 1. This results in distortion of the solution. When the value of Kp is low, the viscosity of the scheme considerably smoothes off the solution. For these two reasons, it is inadvisable to use this scheme in simulations.
The scheme with implicit approximation of lowest terms is derived by replacing the stresses ar, av and the velocity Vr in the lowest terms with ar, av and Vr. This scheme is stable and monotone for 0 < Kp < 1 but it also smoothes off the solution at sufficiently low values of the Courant number.
The scheme with implicit approximation of lowest terms by the Crank-Nicolson method includes
ar + ar + av Vr + Vr
2 ' 2 ' 2 . In regard to the accuracy of numerical solution, this scheme has certain advantages over the explicit and implicit schemes.
3. Results of computations
Computational schemes were verified by comparing the results of simulations with the exact solution obtained for the monochromatic wave with frequency w with the use of the method of separation of variables. The exact solution has the form
a0 • . T If N a0 „„„ Y/A , o \ T If \ 2(A + P)
Vr =- sin wt.Ji(£p), ar = cos wt[ (A + 2 p) .h(Çp)--r-Ji(Çp) ,
Pcp A + 2 P \ Çp J
a = a0 cos wt\AJ2(Çp) - 2 (A+ p) Ji (Çp M , az = t4^ cos wt( J2Ç) - 2 Ji(Çp )
* A + 2 j ^ < a + 2 j V ' S
where Sp = ur/cp is the dimensionless variable, Jk (x) is the Bessel function of an integer order k.
Tabs. 1-4 present relative errors of the schemes for various frequencies as a function of the Courant number. The dimensionless frequency u = uR/cs, where R is the radius of the
computational domain, was varied between 10 and 50. At such frequencies the number of half-waves in the computational domain varies between one and a half and seven and a half (Fig. 1).
Table 1. Relative errors for the dissipation-free scheme (7 = 0)
^^^ CJ 10 20 30 40 50
0.5 0.00023 0.00110 0.00837 0.01408 0.03500
0.75 0.00009 0.00045 0.00405 0.00701 0.01818
1 0.00019 0.00081 0.00201 0.00349 0.00563
1.25 0.00049 0.00210 0.00977 0.01621 0.03527
1.5 0.00085 0.00372 0.01920 0.03209 0.07125
Table 2. Relative errors for the scheme with explicit approximation of lowest terms
^^^ CJ 10 20 30 40 50
0.5 0.04334 0.09600 0.37513 0.40646 0.63148
0.75 0.02148 0.05053 0.19970 0.25226 0.39363
Table 3. Relative errors for the scheme with implicit approximation of lowest terms
^^^ CJ 10 20 30 40 50
0.5 0.04802 0.10547 0.39297 0.41076 0.65137
0.75 0.02860 0.06545 0.23341 0.26373 0.44424
1 0.01025 0.02666 0.03320 0.04613 0.05591
Table 4. Relative errors forn the scheme with Crank-Nicolson approximation of lowest terms
^^^ CJ 10 20 30 40 50
0.5 0.03824 0.09015 0.36576 0.40030 0.62509
0.75 0.01350 0.04288 0.18149 0.23903 0.37536
0.97 0.00912 0.02940 0.02788 0.05726 0.04975
1 0.01220 0.03379 0.05646 0.06663 0.12026
To calculate the error of numerical solution a discrete equivalent of the norm of the space L^(0,T; L2(0,R)) was used:
(v
= sup
0<t<T
p ^ + W )dr2
The time T was set so that the cylindrical longitudinal wave in the interval (0, T) travels a distance 2 R with single reflection from the axis of symmetry.
The finite difference grid has 200 cells. Analysis of the data in the tables shows that numerical solution obtained with implicit approximation of lowest terms and with approximation by
r
Fig. 1. Exact solution for cylindrical longitudinal waves; dimensionless velocity versus dimensionless distance
the Crank-Nicolson method is inaccurate if one half-wave contains less than 60-70 cells. The dissipation-free scheme gives more accurate results at all frequencies within the specified range.
Figs. 2 and 3 show the velocity profiles behind the front of a strong discontinuity when sudden constant stress is applied at the boundary of the domain. These results are obtained using the scheme with Crank-Nicolson approximation (Fig. 2) and the dissipation-free scheme (Fig. 3).
Fig. 2. Velocity profiles behind the front of discontinuity: scheme with approximation of lowest terms by the Crank-Nicolson method
Fig. 3. Velocity profiles behind the front of discontinuity: dissipation-free scheme
In the case of the dissipation-free scheme, the profiles of velocities and stresses are only monotonous at Kp = l, while at Kp = 0.9 parasitic oscillations appear ahead of the wave front, and they grow as the Courant number decreases. The oscillations are smoothed off setting j > 0 in doing so the artificial dissipation of energy is introduced or replacing the sudden application of stress at the boundary by a monotone increasing stress that is changed from zero to a constant value during not less than 10 time steps.
To implement the scheme with controllable energy dissipation on computational clusters it is possible to use the iterative process which demonstrates rarely high rate of convergence of approximate solutions in umerical experiments. It appears that errors presented in Tab. 1 are already obtained with one or two iterations.
The problem consists in parallel computing at the "predictor" stage of the finite difference scheme. System of equations (7) is changed at the junction points of the neighbor processors and corresponding three-point equations of the system are replaced by the relations of the Godunov scheme:
. vrj+l/2 + vrj-l/2 arj+l/2 — ar j — 1/2
vr =---+
2 2 pcp
where fractional indices mark the velocities and stresses that belong to the boundary of grid cell of the neighbor processors. This procedure allows one to implement the three-point sweep method in parallel mode on a cluster and obtain solution in the first approximation.
Referring to figures given above, the nonconservative finite difference scheme with approximation of lowest terms by the Crank-Nicolson method provides much more reliable results for solutions with discontinuities over the whole range of the Courant number Kp < 1 (this is stability condition of the scheme).
It is worth to mention that the analogous schemes with conservative equations (3) inadequately distort the pattern of wave reflection from the axis of symmetry even in the case of smooth solutions, and this may finally result in the total loss of accuracy.
Conclusions
The finite difference scheme with controllable dissipation of energy and typical grid-characteristic schemes of the "predictor-corrector" type were considered. The comparison of the results of computations shows that the scheme with controllable dissipation has undisputable advantages over the other approaches in the case of smooth solutions. As for solutions with discontinuities, the grid-characteristics schemes are preferable due to their monotonicity. In the case of solutions with discontinuities the scheme with controllable energy dissipation produces parasitic oscillations. To smooth off these oscillations one should introduce artificial energy dissipation.
This work is supported by the Krasnoyarsk Mathematical Center and financed by the Ministry of Science and Higher Education of the Russian Federation in the framework of the establishment and development of regional Centers for Mathematics Research and Education (Agreement No. 075-02-2020-1631). Computations were performed using computational cluster at the Institute of Computational Modeling SB RAS, Krasnoyarsk, and at the Joint Supercomputer Center of the Russian Academy of Sciences, Moscow.
References
[1] S.K.Godunov, A.V.Zabrodin, M.Ya.Ivanov, A.N.Kraiko, G.P.Prokopov, Numerical Solution of Multidimensional Problems of Gas Dynamics, Moscow, Nauka, 1976 (in Russian).
[2] O.Sadovskaya, V.Sadovskii, Mathematical Modeling in Mechanics of Granular Materials, Ser. Advanced Structured Materials, vol. 21, Heidelberg - New York - Dordrecht - London, Springer, 2012. DOI: 10.1007/978-3-642-29053-4
[3] O.V.Sadovskaya, V.M.Sadovskii, Numerical implementation of mathematical model of the dynamics of a porous medium on supercomputers of cluster architecture, AIP Conf. Proc., vol. 1684, 2015, 070005-1-070005-9. DOI: 10.1063/1.4934306
[4] V.M.Sadovskii, O.V.Sadovskaya, Analyzing the deformation of a porous medium with account for the collapse of pores, J. Appl. Mech. Techn. Phys., 57(2016), no. 5, 808-818. DOI: 10.1134/S0021894416050072
[5] V.M.Sadovskii, O.V.Sadovskaya, Modeling of elastic waves in a blocky medium based on equations of the Cosserat continuum, Wave Motion, 52(2015), 138-150.
DOI: 10.1016/j.wavemoti.2014.09.008
[6] V.M.Sadovskii, O.V.Sadovskaya, A.A.Lukyanov, Modeling of wave processes in blocky media with porous and fluid-saturated interlayers, J. Comput. Phys., 345(2017), 834-855. DOI: 10.1016/j.jcp.2017.06.001
[7] V.M.Sadovskii, O.V.Sadovskaya, Supercomputer Modeling of Wave Propagation in Blocky Media Accounting Fractures of Interlayers, In: Nonlinear Wave Dynamics of Materials and Structures (Eds. H.Altenbach, V.A.Eremeyev, I.S.Pavlov, A.V.Porubov), Ser. Advanced Structured Materials, vol. 122, chapt. 22, Cham, Springer, 2020, 379-398.
DOI: 10.1007/978-3-030-38708-2_22
[8] P.F.Sabodash, R.A.Cherednichenko, Application of the method of spatial characteristics to the solution of axisymmetric problems on the propagation of elastic waves, Prikl. Mekh. Tehn. Fiz., 12(1971), no. 4, 101-109 (in Russian).
[9] K.M.Magomedov, A.S.Kholodov, Grid-Characteristic Numerical Methods, Moscow, Nauka, 1988 (in Russian).
[10] V.N.Kukudzhanov, Difference Methods for the Solution of Problems of Mechanics of De-formable Media, Moscow, MFTI, 1992 (in Russian).
[11] V.N.Kukudzhanov, Numerical Continuum Mechanics, Ser. De Gruyter Studies in Mathematical Physics, vol. 15, Berlin - Boston, De Gruyter, 2013. DOI: 10.1515/9783110273380
[12] N.G.Burago, A.B.Zhuravlev, I.S.Nikitin, Continuum model and method of calculating for dynamics of inelastic layered medium, Math. Models Comput. Simul., 11(2019), no. 3, 488-498. DOI: 10.1134/S2070048219030098
[13] I.S.Nikitin, Theory of Inelastic Layered and Blocky Media, Moscow, Fizmatlit, 2019 (in Russian).
[14] N.G.Burago, V.N.Kukudzhanov, Buckling and Supercritical Deformations of Elastic-Plastic Shells under Axial Symmetry, In: Collection of Numerical Methods in the Mechanics of Deformable Solids, Moscow, Computing Center of the USSR Academy of Sciences, 1978, 47-66 (in Russian).
[15] V.G.Bazhenov, E.V.Igonicheva, Nonlinear Processes of Shock Buckling of Elastic Structural Elements in the Form of Orthotropic Shells of Rotation, N. Novgorod, UNN, 1991 (in Russian).
[16] V.G.Bazhenov, D.T.Chekmarev, Solution of the Problems of Non-Stationary Dynamics of Plates and Shells by the Variational-Difference Method, N. Novgorod, UNN, 2000 (in Russian).
[17] G.I.Marchuk, Splitting and Alternating Direction Methods, In: Handbook of Numerical Analysis, vol. 1 (Eds. P.G. Ciarlet, J.-L. Lions), North-Holland, Elsevier, 1990, 197-462. DOI: 10.1016/S1570-8659(05)80035-3
[18] A.G.Kulikovskii, N.V.Pogorelov, A.Yu.Semenov, Mathematical Aspects of Numerical Solution of Hyperbolic Systems, Ser. Monographs and Surveys in Pure and Applied Mathematics, vol. 118, Boca Raton - London - New York - Wasington, Chapman & Hall/CRC, 2001. DOI: 10.1201/9781482273991
[19] Yu.P.Popov, A.A.Samarskii, Completely conservative difference schemes, USSR Comput. Math. Math. Phys, 9(1969), no. 4, 296-305. DOI: 10.1016/0041-5553(69)90049-4
[20] A.A.Samarskii, The Theory of Difference Schemes, Ser. Pure and Applied Mathematics, New York - Basel, CRC Press, 2001. DOI: 10.1201/9780203908518
[21] B.L.Rozdestvenskii, N.N.Janenko, Systems of Quasilinear Equations and Their Applications to Gas Dynamics, Ser. Translations of Mathematical Monographs, vol. 55, Providence, American Mathematical Society, 1983.
[22] G.V.Ivanov, V.D.Kurguzov, Schemes for solving one-dimensional problems of the dynamics of inhomogeneous elastic bodies on the basis of approximation by linear polynomials, Dynam. Contin. Medium, 49(1981), 27-44 (in Russian).
[23] G.V.Ivanov, Yu.M.Volchkov, I.O.Bogulskii, S.A.Anisimov, V.D.Kurguzov, Numerical Solution of Dynamic Elastic-Plastic Problems of Deformable Solids, Novosibirsk, Sib. Univ. Izd., 2002 (in Russian).
Разностные схемы для анализа продольных волн на основе осесимметричных уравнений динамической теории упругости
Владимир М. Садовский Оксана В. Садовская Евгений А. Ефимов
Институт вычислительного моделирования СО РАН Красноярск, Российская Федерация
Аннотация. Цель исследования состоит в построении экономичной разностной схемы сквозного счета для решения прямых задач сейсмики на основе уравнений динамики упругой среды в осесим-метричной постановке. При численной реализации схемы на многопроцессорных вычислительных системах применяется метод двуциклического расщепления по пространственным переменным. Одномерные системы уравнений на этапах расщепления распадаются на подсистемы продольных, поперечных и крутильных волн. В данной работе рассматривается случай продольных волн. Проводится сравнение явных сеточно-характеристических схем и неявных схем типа "предиктор-корректор" с контролируемой диссипацией энергии на точных решениях, описывающих бегущие монохроматические волны.
Ключевые слова: упругая среда, цилиндрические волны, метод расщепления, разностная схема, монотонность, диссипативность, параллельная реализация.
