EXPLICIT-IMPLICIT SCHEMES FOR CALCULATING THE DYNAMICS OF LAYERED MEDIA WITH NONLINEAR CONDITIONS AT CONTACT BOUNDARIES Текст научной статьи по специальности «Физика»

DOI: 10.17516/1997-1397-2021-14-6-768-778 УДК 519.63

Explicit-Implicit Schemes for Calculating the Dynamics of Layered Media with Nonlinear Conditions at Contact Boundaries

Ilia S. Nikitin*

Institute of Computer Aided Design RAS Moscow, Russian Federation

Vasily I. GolubeV

Moscow Institute of Physics and Technology (National Research University)

Moscow Region, Russian Federation

Received 10.06.2021, received in revised form 27.07.2021, accepted 15.09.2021 Abstract. In this paper we consider the problem of dynamic loading of a deformable solid medium containing slip planes with nonlinear slip conditions on them. An explicit-implicit scheme was constructed for the numerical solution of the constitutive system of equations, which exactly reduces to correcting the stress tensor values after performing the elastic step. An implicit approximation of the constitutive relations containing a small parameter in the denominator of the nonlinear free term was used with the second order of the approximation. The correction procedure is applicable for those cases when the viscosity parameter of interlayers providing the sliding mode of the contact boundaries is not small. The solution of the problem of the seismic waves propagation in an inhomogeneous fractured geological massif in a two-dimensional case was obtained numerically.

Keywords: layered media, fractured media, numerical simulation, explicit-implicit method, grid-characteristic method, slip conditions.

Citation: I.S. Nikitin, V.I. Golubev, Explicit-Implicit Schemes for Calculating the Dynamics of Layered Media with Nonlinear Conditions at Contact Boundaries, J. Sib. Fed. Univ. Math. Phys., 2021, 14(6), 768-778. DOI: 10.17516/1997-1397-2021-14-6-768-778.

Seismic exploration is a standard method of prospecting for hydrocarbon deposits and is actively used by many service companies in the oil and gas complex. It is based on the propagation of seismic waves in heterogeneous media, their reflection from contrasting boundaries, during the registration of which the internal structure of the subsurface space can be restored. Computational methods that allow calculating wavefields in a given geological model is a crucially important. With the development of modern computer systems, they are being actively improved and developed.

Various mathematical models are used to describe dynamic processes in geological media: acoustic, isotropic elastic, anisotropic elastic, elastic-viscoplastic. Unfortunately, the vast majority of applied problems cannot be solved analytically, in view of which various numerical methods are actively used [1-4]. Well prepared reviews are presented in [5,6]. Among recent works, one can highlight the work [7], devoted to the construction of a hybrid numerical method;

* i_nikitin@list.ru https://orcid.org/0000-0003-3499-6910

tw.golubev@mail.ru, golubev.vi@mipt.ru https://orcid.org/0000-0003-3113-7299 © Siberian Federal University. All rights reserved

and [8,9], in which a numerical method was constructed for calculating wave processes in a block medium with interacting structural elements. In paper [10], the grid-characteristic method was successfully generalized for inverse migration problems.

Fractured media are of the practical interest, since they either contain oil or are confined to the places of occurrence of its large reserves. Earlier in works [11,12], continuous models of deformable solid media with a discrete set of slip planes (layered, block media) and with nonlinear slip conditions at contact boundaries of structural elements were constructed. In all these cases, the constitutive relations of the system include equations with a nonlinear free term, the form of which is determined by the selected contact condition at the interlayer boundary.

For a stable numerical solution of the system of differential equations, explicit-implicit schemes were developed with an explicit approximation of motion equations and an implicit approximation of constitutive relations containing a small parameter in the denominator of the nonlinear free term [13,14]. Various effective formulas for correcting stress components after an "elastic" time step were analytically derived from implicit nonlinear finite-difference approximations. To calculate the "elastic" step, the grid-characteristic method on structured grids was used. For monotonization, a grid-characteristic criterion of the monotonicity was used, which is based on the characteristic property of the exact solution of the linear transport equation.

However, while an explicit elastic step allows the construction of schemes of a higher order of approximation [15,16], the correction of stresses in previously published works was carried out by means of an implicit first-order approximation [13,14]. In this paper, to match the orders of approximation of the explicit elastic and implicit correction steps, an implicit second-order approximation was constructed for the constitutive equations of a layered medium with friction at contact boundaries and refined correction formulas were obtained after the elastic calculation step. The use of this method with the consistent approximation for explicit elastic and implicit correction steps made it possible to increase the accuracy of calculations and to carry out a numerical simulation of non-stationary problems of scattering of elastic waves and generation of a response from fractured clusters, both in a homogeneous and in an inhomogeneous geological massif.

1. Mathematical model and numerical method

Let us formulate nonlinear conditions for the interaction of contact boundaries of structural elements. In a Cartesian rectangular coordinate system xi (i = 1,2,3), we consider an infinite elastic medium with an oriented system of periodically repeating parallel slip planes. The orientation of this system is set by the unit normal n. The distance between sliding planes is constant and equal e. The density of the material p, as well as Lame moduli of elasticity A and n are assumed to be given constants. The stress state is described by the stress tensor a. The vector of the shear stress on the slip plane is t = a ■ n — (n ■ a ■ n), the normal stress is an = n ■ a ■ n. We assume that layers are in a compressed state and an < 0.

Let us introduce the vector of the shear velocity 7, determined by the jump of the tangential

velocity [VT] at the contact boundary: 7 = i.

e

It is assumed that physically between elastic layers there are thin layers with the thickness d C e, however, we neglect the thickness of these layers and replace them with sliding conditions at compressed boundaries of layers. The conditions of the contact interaction, depending on physicomechanical properties of interlayers, are taken in the form the of local slip condition for Coulomb friction with a small viscous additive.

At an < 0 (compressed contact boundaries):

t = q\<n\(y/\y\ + nY), (i)

or expressing the slip velocity y in terms of shear stresses t

Y = (2)

where n is the viscosity coefficient, q is the Coulomb friction coefficient, (F(y)) = F(y)H(y), H(y) is the Heaviside function, H(y) =0 if y < 0 and H(y) = 1 if y > 0. The contact plane with indicated interaction conditions is called the slip plane.

In order to deal with the continuous model of a medium containing a system of such slip planes, we will consider y as a continuous function of coordinates and time, which has the meaning of distributed slip velocities. These relations allow us to take into account the contribution of slip velocities to the tensor of velocities of an inelastic deformation eY:

eY = (n ® y + Y ® n) /2, Y • n = 0. (3)

The total strain velocity tensor e is obtained by summing all elastic and inelastic components and is equal to

e = ee + eY, e = (Vv + VvT) /2. (4)

Here v is the "macroscopic" velocity of medium particles, ee is the elastic strain velocity tensor, which is related to the stress tensor by Hooke's law:

a = X (ee : I) I + 2^ee. (5)

The system is closed by motion equations:

pV = V • a. (6)

The system of equations for the model of a layered medium.

In a layered medium consisting of elastic layers, there is a single system of slip-delamination planes with a normal n. If the normal n to interlayer boundaries is oriented along the x3 coordinate axis, then its components satisfy the relation nj = S3.

The conditions for y, corresponding to the local contact condition, have the form:

* = 13« - 0 • \T\ = ^k. (7)

Taking into account the choice of the orientation of the normal in a component wise form, the system of equations of the model will take the form:

pVi = Oijj ,&u = Xvk,k +2^vi,i, <733 = Xvk,k + 2^v3,3, (8)

i=3

<ij . = M (vi,j + vj,i), i= j, <3j = M (v3,j + vj,3) - Vlj. (9)

i,j=3 j=3

The system of equations for the model of a block medium.

The block medium is formed by uniformly stacked elastic cubes (parallelepipeds) with three possible slip-delamination planes oriented by mutually perpendicular unit normals n(s),

s = 1,2, 3. In this case, the tensor of the velocity of an inelastic deformation will take the form:

3

eY = ^2 (n(s) ® Y(s) + Y(s) ® /2, Y(s) • n(s) = 0. (10)

If we orient three normals to slip-delamination planes along coordinate axes of the Cartesian

where Ss

coordinate system, then the following relation will be true for normals: nj = Ss, where Ss is

the Kronecker symbol.

The through conditions for Y(i), corresponding to the local contact conditions, have the form:

j=1A (£n1) •■=j Jzr*- (11)

As in the case of a layered medium, we can write the system of equations for a block medium in a component wise form:

pvi = , &jj = Avk,k + 2^vj,j, (12)

¿ij = M (vi,j + vj,i) — Mljl) — MYj, i = j- (13)

1.1. Implicit approximation of constitutive relations

The obtained system belongs to the class of semilinear hyperbolic systems; their numerical solution can be constructed using various explicit schemes. However, in the slip mode, a nonlinear free term with a small viscosity parameter in the denominator is included in the governing equations for shear stresses. The system becomes rigid and common explicit schemes will be unstable. To overcome these difficulties, we propose to use an explicit-implicit method. Implicit approximation is used only for those equations that contain a small parameter in the denominator of the free term, the rest of equations are approximated explicitly.

Let us describe this method using the example of the equation for ¿3j in the compressed contact boundary mode ¿33 < 0 for a layered medium:

¿3j = M (v3,j + vj,3) — M¿3j {\t\/(q\^ss\) — 1 / (n\T I), \t \ = V¿i3¿i3, i = 1, 2- (14)

3=3

Earlier, as a result of the analytical solution of an algebraic equation for an implicit first-order approximation, the following correction formulas were obtained after an elastic step for the case of compressed layers [13, 14]:

When \Tn+'\ > q\an+e\

¿in3+1 = q\¿nз+e1\ (¿in3+7\Ten+1\) (1 + SlTn+'l) / (1 + Squill) , (15)

Yi = (¿nl1 — ¿23+1) / (mA), ¿i3 = ¿13 + MAt (vnl1 + v3n+*), i = 1,2. (16)

When \Tn+1\ < qWn+I

¿n1 = ¿Hi\ Y?+1 =0,S = v/ (MAt), \Tn+1\ = yl¿n!t1¿inз+1, i = 1, 2. (17)

Here, indices n +1 and n correspond to the values on the upper and lower time layers, At is the time step, the subscript e marks values of the stress tensor after the elastic time step. It is assumed that the values of vrn+1 and ¿^l1 have already been determined from an explicit approximation of motion equations and equations for the normal component of the stress tensor that do not contain small viscosity parameters in the denominator of the free term. To agree approximation orders of the elastic step and the corrective numerical procedure, we construct an implicit second-order scheme for governing equations of a layered medium with nonlinear contact conditions at interlayer boundaries. We do not write out further the terms related to

the approximation of spatial derivatives, since it can be carried out by various schemes. A specific implementation based on the use of the inverse method of characteristics, which has a second order of approximation in 2D and 3D cases, is given in the next section. An implicit approximation in time of the second order for the equation for tangential stresses is:

n+1 Ji3

'i3

At

( n+1 + n+1N W/im+l A ^"T1 ,/ |Tn| A ^3 1

+ }" n\V nF+M +\ V nR/

or

i3

i3

u/im+i_ A ^3+1+ /

2n\q|^3n3+1| / Itn+1| + \

It nl ,qK31

-1

/ It"U

This nonlinear system of algebraic equations for ^n^1 can be rewritten as:

^+2

[/ IT"+1| A / J^L A

Wql^l / |t"+1I + \ql^l / lTnij

Let us collapse this equation sequentially with anj+1, ^¡3, af3 and introduce notations: for unknown convolutions at the calculation step

X2 =

for known quantities

a3+1_a3+1

°i3 °i3 ,

Y2

^n _n+1 r/2 a"3 • , Z

e n+1 ai3 • °i3 ,

1 = ai3 • ai3> S = ai3 • ai3> S = ai3 • ai3-

We get three algebraic equations for X, Y and Z:

1

^ 2 + 2 sz2+2

X 1 -1 X2 + / T 1 -1 \ Y2

^?3+1| X \q|o33l ) t

X -1) \ Y2 + / T -1) T2

K3+1! ' X WK3I ) T

X -1) )Z2 + / T -1) S2

K3+1! ' X \q|°33| ' T .

SZ2

^ = SS2

^ = SS2.

We will also denote for brevity:

AT = 1

(— - l) 2\ ql^l /

1

- 1) /T, AX = S +-

2( q|o"3+1l 1

- 1 /X.

In this case, the algebraic system can be written in the form:

X 2AX + Y2 AT = 5Z2, Y2 AX + T2 AT = 5S2, Z2 AX + S2AT = 5T?. From second and third equations, we can obtain:

5S2 - T2 AT

Y2

AX

Z2

SS2 - S2 AT

AX .

(18)

(19)

(20)

(21) (22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

IL

}

From the first equation we obtain:

*2 AX = ^^^ — ATSS—^ (31)

or

X2 AX2 = AP2, AP2 = S (SS2 — S2AT) — AT (SS2 — T2AT). (32)

The system for defining the convolution X will take the form:

XAX = AP, (33)

1 ( X

AX = S +-( — A /X, (34)

2\ q^l /1 ()

From this system

SX + -1 ^ — A = AP, (35)

2\ qWn3+1\ / , ()

and when X > \ ¿3+11 we get

X = 1 + 2AP+1 ). (36)

2S +1/(qWn++1\) ( )

So, for this convolution we have the expression:

X = ^T^l • AX = %. (37)

The difference equation for ¿n^1 can be written in the form

¿n+1AX + ¿n AT = S¿iз (38)

or

n+1 S¿i3 — ¿¡3AT YS_at1~anAT

¿i3 =-AX-= X-AP-• (39)

The final correction formula of the second order of approximation will take the form:

¿n+1 ^¿n+1\ 1 + 2AP ¿¿^ — AT¿3з (40)

¿i3 = q¿3 11 + 2SqK3+1\-AP-, (40)

AP = VS (SS2 — S2 AT) — AT (SS2 — AT ■ T2), (41)

AT=*(w—') sm- S=^ (42)

S = T = ^¿¿3, i =1,2• (43)

The correction formula for shear stresses can be written in a different, more convenient form. For this purpose, we introduce the notation ¿i3+1 for additional components with corrections to the elastic calculated values containing, at first glance, a "bad" small parameter in the denominator:

A 1 / \ Ti

¿"3+1 = ¿i3 — -¿i3 = ¿e3— 2S (q® — V T3, T = ^3- (44) We also introduce a notation for the convolution:

= y/n+1

'„i3 „i3

(45)

It is easy to show that AP = SE"+1. So, in the square root expression in the formula for AP the value is always positive. In these notations, the correction formula for shear stresses on the n + 1 time layer will take the form:

1 + 26 V n+1 V n+1 „n+1 = q|„n+1| 1 + 26 ^_Vi3

ffi3 = ql„33 11 + 26q\„n+'1\ Vn+1

(46)

Let us show that the small parameter in the denominator of the expression for E2+1 is canceled after some transformations:

X =

= / n+1 n+1 = I n+11 1 + 26V n+1 ' i3 „ i3 = q\„33 L 0J- | n+11 1 + 26q\„33+ \

, T = v„„3 = q\„n3

1 + 26 Vn 1 + 26q\„ 33 \

^^ ' ¿(qR3T 1) = 26(

26

(sn - q\„J3\ (sn - q\„?3\)

1 + 26q| 3n3 |

1 + 26q| 3n3 |

and, finally,

-n+1 = 7i3 = '

i3

Vn - q| 3n3|

i3

(47)

(48)

(49)

1 + 2Sq\a^l •

It can be seen from this formula that differs from the elastic stress on the upper time layer by a relatively small value, and the formula for a"+1 does not contain a small viscosity parameter in the denominator.

Note that the obtained nonlinear algebraic system from the implicit approximation was solved exactly, without taking into account the smallness of the parameter S. Therefore, the proposed correction procedure is also applicable for those cases when the viscosity parameter of interlayers providing the sliding mode of contact boundaries is not small.

1.2. Explicit approximation of motion equations

Let us consider in more details the numerical method used to solve the elastic part of the problem. In the three-dimensional case, the vector of unknowns contains nine components

u = (v1,v2,v3, a 11, a12, a13, a22, a23, a33) . Note that the original system of equations can be written in the canonical form

du „ du „ du „ du

UT = + + ,

dt dx1 dx2 dx3

(50)

(51)

where A1, A2, A3 are 9 x 9 matrices and (x1,x2,x3) is an orthonormal coordinate system. After splitting in coordinate directions, the original system splits into three independent one-dimensional systems

du = Aj = 1,2,3. (52)

dt

Each of them is a hyperbolic system of equations, which means that its matrix can be diagonalized

du

% = ^dXj, j = 1,2,3,

(53)

n

n

e

where Qj is composed of eigenvectors and Aj is a diagonal matrix composed of eigenvalues

Aj = diag (ci, —ci, C2, — C2,C2, — C2, 0,0,0), j = 1, 2, 3, (54)

and c\ = 4 j + and c2 = 4 . After passing to the Riemann invariants v = Qj u, each of

VP VP

the systems splits into a set of linear transport equations

dv dv ,

m + A âT = 0' j = 1'23 (55)

The transport equation is solved using an one-dimensional interpolation procedure with a given order of approximation. Then the solution is restored on a new time layer. Note that to calculate the dynamic behavior of a layered (or block) medium with nonlinear slip at the contact boundaries, it is enough to apply the correcting relations given above after each elastic time step.

2. Simulation results

On the basis of the proposed computational algorithm consistent in the order of approximation, the numerical simulation of multidimensional problems of reflection and scattering of elastic waves from a near-surface source on fractured buried geostructures was carried out.

At the first stage, a comparison was made of the calculation results obtained using the first and the second order of approximation. For this, a two-dimensional problem of isotropic geological massif loading with sizes of 2500 x 3000 m, a velocity of propagation of longitudinal waves of 4500 m/s, velocity of propagation of transverse waves of 2250 m/s, and density of 2500 kg/m3 was considered. At a depth of 2000 m, a fractured cluster with a length of 1500 x 100 m was located in it. The dynamic behavior of the medium in this region was described within the framework of a continual model with parameters q = 0.1, ó = 0.3, n = ^^, . It was illuminated by a longitudinal plane wave with a harmonic time dependence with a frequency of 30 Hz and an amplitude of 1 cm/s. A rectangular computational grid was built with a spatial step of 5 m, covering the computational domain and containing 300 000 nodes. The time step was chosen from the Courant condition and was 1 ms. In total, 1 s of the loading process was simulated. For the numerical solution of one-dimensional linear transport equations, the Rusanov scheme of the third order of accuracy was used. In this 2D case, due to the use of coordinate splitting technique, the total order of the elastic step is reduced to the second.

Fig. 1 shows the difference in seismic signals obtained using both computational algorithms. In general, when the P-wave interacts with inclined cracks, reflected longitudinal and transverse waves are generated — from the horizontal part of the inclusion, as well as two fronts of spherical waves — from the side boundaries of the inhomogeneity. Analysis of the results shows that the difference in seismic responses reaches 0.8%, and it is the most significant for the P-S response.

At the second stage, the seismic wave propagation was calculated in the Marmousi2 model [14], complicated by the presence of a cluster of vertically oriented fractures in the volume [4000 — 4500] x [1000 — 1500] m. The simulation domain was covered with a rectangular computational grid with a spatial step of 5 m, in each node of which the elastic parameters of the medium were stored. A spherically symmetric near-surface explosion with a time dependence in the form of a Ricker pulse with a main frequency of 30 Hz was considered as a signal source. It was located at the point (8500,5). The fractured cluster was described by the following parameters: q = 0.1, ó = 0.3, n = (1,0). Fig. 2 (top) shows the wavefield superimposed on the

Fig. 1. Seismic response (left) and absolute value of velocity modulus (right) for first order and second order schemes usage

density model of the medium. One can see the formation of numerous waves reflected from curvilinear contact boundaries. At the same time, it is not possible to identify the position of the heterogeneity itself. Fig. 2 (bottom) shows the difference wavefield obtained by subtracting the simulation results for a fully elastic model. The generation of spherically diverging response waves is clearly visible. Note that the level of the response amplitude is also at the level of 1%.

Conclusions

In this paper, we consider the problem of dynamic loading of a deformable elastic medium with a discrete set of slip planes and nonlinear slip conditions at the contact boundaries of structural elements. A continual model was used to describe its behavior. For a stable numerical solution of the governing system of differential equations, an explicit-implicit scheme with an explicit approximation of the motion equations and an implicit approximation of the constitutive relations containing a small parameter in the denominator of nonlinear free terms is used. To match the orders of approximation of the explicit elastic and implicit correction steps, an implicit second-order approximation is constructed for the constitutive equations of a layered medium with friction at the contact boundaries, and refined correction formulas are obtained after the elastic calculation step. To calculate the "elastic" step, the grid-characteristic method on rectangular grids was used. The use of the method of consistent approximation made it possible to increase the accuracy of calculations and to carry out a numerical simulation of non-stationary problems of the scattering of elastic waves and the generation of a response from an oriented fractured cluster, both in a homogeneous and inhomogeneous geological massif.

This work was carried out with the financial support of the Russian Science Foundation, project no. 19-71-10060.

Fig. 2. Wave field in Marmousi2 fractured model (top) and difference with pure elastic model (bottom)

References

[1] M.Dumbser, M.Kaser, J.De La Puente, Arbitrary high-order finite volume schemes for seismic wave propagation on unstructured meshes in 2D and 3D, Geophysical Journal International 171(2007), 665-694. DOI: 10.1111/j.1365-246X.2007.03421.x

[2] J.T.Etgen, M.J.O'Brien, Computational methods for large-scale 3D acoustic finite-difference modeling: a tutorial, Geophysics, 72(2007), SM223-SM230. DOI: 10.1190/1.2753753

[3] S.Hestholm, Acoustic VTI modeling using high-order finite differences, Geophysics, 74(2009), T67-T73. DOI: 10.1190/1.3157242

[4] J.W.D.Hobro, C.H.Chapman, J.O.A.Robertsson, A method for correcting acoustic finite-difference amplitudes for elastic effects, Geophysics, 79(2014), T243-T255.

DOI: 10.1190/geo2013-0335.1

[5] J.Virieux, H.Calandra, R.E.Plessix, A Review of the Spectral, Pseudo-Spectral, Finite-Difference and Finite-Element Modelling Techniques for Geophysical Imaging, Geophysical Prospecting, 59(2011), no. 5, 794-813. DOI:10.1111/j.1365-2478.2011.00967.x

[6] J.M.Carcione, C.Herman Gerard, P.E. Kroode, Y2K Review Article: Seismic Modeling, Rev. Lit. Arts Amer., 67(2002), no. 4, 1304-1325.

[7] V.Lisitsa, V.Tcheverda, C.Botter, Combination of the Discontinuous Galerkin Method with Finite Differences for Simulation of Seismic Wave Propagation, J. Comput. Phys., 311(2016), 142-157. DOI: 10.1016/j.jcp.2016.02.005

[8] V.M.Sadovskii, O.V.Sadovskaya, Supercomputer Modeling of Wave Propagation in Blocky Media Accounting Fractures of Interlayers, Advanced Structured Materials, 122(2020), 379-398.

[9] V.M.Sadovskii, O.V.Sadovskaya, I.E.Petrakov, On the theory of constitutive equations for composites with different resistance in compression and tension, Composite Structures, 268(2021), paper 113921.

[10] V.I.Golubev, The Usage of Grid-Characteristic Method in Seismic Migration Problems, Smart Innovation, Systems and Technologies, 133(2019), 143-155.

[11] I.S.Nikitin, Dynamic models of layered and block media with slip, friction, and separation, Mechanics of Solids, 43(2008), no. 4, 652-661.

[12] N.G.Burago, I.S.Nikitin, Improved model of a layered medium with slip on the contact boundaries, Journal of Applied Mathematics and Mechanics, 80(2016), no. 2, 164-172. DOI: 10.1016/j.jappmathmech.2016.06.010

[13] V.Golubev, I.Nikitin, Y.Golubeva, I.Petrov, Numerical simulation of the dynamic loading process of initially damaged media, AIP Conference Proceedings, 2309(2020), paper 0033949.

[14] V.Golubev, I.Nikitin, A.Ekimenko, Simulation of seismic responses from fractured MAR-MOUSI2 model, AIP Conference Proceedings, 2312(2020), paper 050006.

[15] V.I.Golubev, N.I.Khokhlov, I.S.Nikitin, M.A.Churyakov, Application of compact grid-characteristic schemes for acoustic problems, Journal of Physics: Conference Series, 1479(2020), no. 1, paper 012058.

[16] V.Golubev, A.Shevchenko, I.Petrov, Simulation of Seismic Wave Propagation in a Multi-component Oil Deposit Model, International Journal of Applied Mechanics, 12(8)(2020), paper 2050084.

Явно-неявные схемы для расчета динамики слоистых сред с нелинейными условиями на контактных границах

Илья С. Никитин

Институт автоматизации проектирования РАН Москва, Российская Федерация

Василий И. Голубев

Московский физико-технический институт Московская область, Российская Федерация

Аннотация. В настоящей работе рассматривается задача динамического нагружения деформируемой твердой среды, содержащей плоскости скольжения с нелинейными условиями проскальзывания на них. Построена явно-неявная схема для численного решения определяющей системы уравнений, в точности сводящаяся к корректировке значений тензора напряжений после выполнения упругого шага. Неявная аппроксимация определяющих соотношений, содержащих малый параметр в знаменателе нелинейных свободных членов, выполнена со вторым порядком. Процедура корректировки применима и для тех случаев, когда параметр вязкости у прослоек, обеспечивающих режим скольжения контактных границ, не мал. Численно получено решение задачи о распространении сейсмических волн в неоднородном трещиноватом геологическом массиве в двумерной постановке.

Ключевые слова: слоистые среды, трещиноватые среды, компьютерное моделирование, явно-неявный метод, сеточно-характеристический метод, условия скольжения.