CC BY
10
УДК 519.63
Д. Н. Ворощук, В. А. Миряха, И. Б. Петров
Московский физико-технический институт (государственный университет)
Исследование 3В-сейсмического отклика от вертикального геологического разлома разрывным методом Галеркина
При помощи математического моделирования исследуется сейсмический отклик от вертикального геологического разлома. Используется разрывный метод Галеркина на частично структурированных тетраэдральных сетках, реализованный для расчета пространственных процессов в гетерогенных средах. Результаты демонстрируют эффективность метода для решения прямых задач сейсморазведки.
Ключевые слова: разрывный метод Галеркина, частично структурированные тетраэдральные сетки, сейсморазведка, численное моделирование.
D. N. Voroshchuk, V. A. Miryaha, I. B. Petrov
Moscow Institute of Physics and Technology (State University)
3D-Seismic response investigation from a vertical crack by Discontinuous Galerkin numerical method
A seismic response from a vertical geological crack is investigated by mathematical modeling. The Discontinuous Galerkin method for a semistructured tetrahedral mesh is implemented for the calculation of space dynamic processes in heterogeneous media. The results obtained show the efficiency of this method for direct seismic problems.
Key words: DGM, semistructured tetrahedral mesh, exploration seismology, numerical modeling.
1. Введение
Несмотря на то, что сейсмическая разведка является самым востребованным методом
поиска полезных ископаемых и более 95% геофизических работ на планете проводится
методом сейсмической разведки, она по-прежнему является самым дорогим геофизическим методом. Последний факт стимулирует использование дополнительных методов, в
частности, математического моделирования. В статье приводится реализация одного из
численных методов на примере прямой задачи сейсморазведки - исследовании сейсмического отклика от вертикального разлома. Разрывный метод Галеркина относится к классу
конечно-элементных методов. Метод обладает рядом положительных свойств, выделяющих его перед другими при решении гиперболических систем уравнений. К ним относится консервативность схемы, гибкость в выборе базиса, Ир-адаптивность. Применение частично структурированных тетраэдральных сеток, оставляя возможность задания криволинейных границы, позволяет существенно ускорить процесс численного расчёта. Благодаря этому расчёты с количеством степеней свободы ~ 6 ■ 107 удалось провести на персональном компьютере, что свидетельствует об эффективности реализованного алгоритма.
2. Постановка задачи
Проводилась серия расчетов. В первом случае рассматривалась однородная упругая среда. В последующих в центральной части области на расстоянии 1500 м от верхней границы располагался вертикальный разлом с геометрическими размерами (глубина х ширина х высота), равными 200 х 20 х 200 м. Общие размеры области интегрирования во всех случаях равнялись (Г х Ш х В) 1800 х 2000 х 2000 м. В данной работе в качестве сей-смосигнала использовалась р-волна с плоским фронтом, направленная вертикально вниз. Середина волнового фронта располагается на глубине 400 м. На рис. 1 приведена общая схема численного эксперимента: И = 2000 м, Е = 1800 м, В = 20 м, А = 200 м, С = 100 м.
О
А
-г-й В
! Е
■1с
Рис. 1. Геометрические параметры модели
Среда описывалась тройкой параметров р,ср,сь - плотность, продольная и поперечная скорости волн. В рассмотренных экспериментах значения параметров были следующими: (р, ср, ) = (2450, 3200м, 1780м) - для вмещающей среды и на 20-40% параметры уменьшались для среды, заполняющей разлом.
3. Описание метода
Подробное описание метода приведено в [1-3]. Остановимся лишь на том частном случае, который использовался для проведения расчетов, рассмотренных в статье. Распространение волновых процессов рассматривается в рамках линейной теории упругости. Для записи основных уравнений используются второй закон Ньютона и следствие из закона Гука:
( дю. Р
да,
г]
01 дХг
_ 1 Ыз = ^
/ + диЛ
\дxj дХг ) '
(1)
где р - плотность среды, щ - смещение точки (вектор), С^ы - тензор упругих напряжений, е^ - тензор деформаций. При рассмотрении линейно-упругой среды тензор упругих напряжений зависит от параметров Ляме. Учитывая факт симметричности тензора деформаций, можно переписать систему уравнений
да,
Р
г]
дЪ да,
дх;
г]
дИ
= Л
\дхк V гз \дх^ дХг) ,
(2)
где - параметры Лямэ в общем случае являются функциями координат, но в рамках текущей статьи будут рассматриваться только изотропные среды, для которых параметры
являются константами. Для трехмерного случая система выглядит следующим образом:
( д , д Л д Л д —*хх - (Л + 2ц) —и - \—V - Х—'ш = 0,
д
д д д Ж*уу -Л^и - (Л + 2ц) ^= °,
д
9 д
д д д д т*" -ЛдХи -Л%(Л + 2ц)= °,
д
рти д рЖу д
рЖш
д
- + —и) = 0,
д_
'дх'
д_ ду' д
д д дЬ** - Ц(Тг" + ду' = °, д д д ^-ц( д~ги + дхш)=°, д д д
*хх Т7~*ху Т7~*хх — 0
д х д д
д
д
д
м *ху м *уу м *ух — 0,
д х д д
д д д
д х д д
(3)
Приведенная система является гиперболической и записана в переменных скорость-напряжение. Приведем матричное представление:
дяР + А дЯ± + в дЯ± ^
ря дх ря ду ря дх
д
+ ^ = 0,
(4)
где Q = (ахх, ауу, охх, аху, ауг, ахх, ьх,и, V, ш),
Ард =
0 0 0 0 0 0 -( Л + 2ц) 0 0 \
0 0 0 0 0 0 - Л 0 0
0 0 0 0 0 0 - Л 0 0
0 0 0 0 0 0 0 ц 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 - ц
1/р 0 0 0 0 0 0 0 0
0 0 0 -1/р 0 0 0 0 0
0 0 0 0 0 -1/р 0 0 0 /
(5)
Врд =
0 0 0 0 0 0 0 - Л 0 \
0 0 0 0 0 0 0 -( Л + 2ц) 0
0 0 0 0 0 0 0 - Л 0
0 0 0 0 0 0 - ц 0 0
0 0 0 0 0 0 0 0 - ц
0 0 0 0 0 0 0 0 0
0 0 0 -1/ р 0 0 0 0 0
0 -1/ р 0 0 0 0 0 0 0
0 0 0 0 -1/р 0 0 0 0 /
(6)
( 0 0 0 0 0 0 0 0 —X
0 0 0 0 0 0 0 0 —Л
0 0 0 0 0 0 0 0 —(Л + 2^)
0 0 0 0 0 0 0 0 0
Срд = 0 0 0 0 0 0 0 —V 0
0 0 0 0 0 0 —V 0 0
0 0 0 0 0 —1/р 0 0 0
0 0 0 0 — 1/р 0 0 0 0
0 0 —1/р 0 0 0 0 0 0
Скорости распространения упругих волн определяются собственными
бианов указанных матриц и равны:
81 = Ср , «2 = "С8, вз = —Сз, 84 = 0, =0,
¿6 = = 0, в7 = с8, = С3, ¿9 = Ср.
(7)
Выражения для скоростей распространения р и волн имеют вид:
Ср
/
Л + 2^ Р '
(9)
Запишем явным образом матрицу, составленную из собственных векторов, соответствующих собственным значениям «1 — зд:
ЯА =
/А + 2^ 0 0 0 0 0 0 0 X + 2^\
Л 0 0 0 1 0 0 0 X
Л 0 0 0 0 1 0 0 X
0 V 0 0 0 0 0 V 0
0 0 0 1 0 0 0 0 0
0 0 V 0 0 0 V 0 0
Ср 0 0 0 0 0 0 0 Ср
0 Сз 0 0 0 0 0 —Сз 0
0 0 Сз 0 0 0 —Сз 0 0
(10)
Для построения численной схемы введем несколько ограничений, в рамках которых будет рассматриваться задача. Представленные далее вычисления проводились на тетраэдральной сетке. Все коэффициенты приведенных ранее матриц являются постоянными в рамках элемента расчетной сетки.
С каждым элементом сетки свяжем систему координат (£, ), в которой оси координат совпадают с тремя ребрами тетраэдра. Зададим формулы перехода их глобальной системы координат в локальную (так будем называть систему, связанную с элементом) следующими формулами:
С = Щ [Х1(у4г3 — УзХь) + Хз(у1Х4 — У4Х1) + Х4(уз21 — У1?з) +
+(У1(^3 — Z4) + Уз(г4 — г{) + У4(г1 — гз))х+ + (Х1 (хА — Хз) + Хз(Х1 — хА) + хА(хз — Х1))у+ +(^1 (уз — ш) + Хз(у4 — У1) + ХА(У1 — Уз)) г],
V = Щ [У1(Х4Х2 — х2хА) + У2(Х1Х4 — ^4^1) + 21 — ^2) +
+ (Ш (¿4 — х2) + У2(%1 — 24) + ^(^2 — ^1))ж+ + (Х1(Х2 — £4) + #2(24 — 21) + £4(21 — ¿2))У+ + Ы(У4 — Ы + Ж2(^1 — Ш) + Ж4 (У2 — Ш)>],
(11)
С
в
где Хг, у г, Zi
( = щ[г\{хзу2 - Х2уз) + г2(х!уз - Хзу{) + 2з(х2у\ - Х1у2)+ +(2/1 (- 2з) + у2(2з - 21) + Уз(^ - 22))х+
+(Х1(2з - 22) + Х2(21 - 2з) + Хз(г2 - Zl))y+
+ (Х1(У2 - Уз) + Х2(уз - ш) + Хз(У1 - ЫЖ | = Х1(у2(24 - Zз) + Уз(г2 - 24) + ш(¿з - ¿2))+ +Х2(2/1 (2з - 24) + 2/з(-г4 - -21) + 2/4(21 - 2з)) + +Хз(2/1 (24 - 22) + у2(г 1 - 24) + 2/4(22 - 21))+ +Х4(2/1(22 - 2з) + 2/2(2з - 21) + 2/з(21 - 22)),
+ Х = Х1 + ( Х - Х1) + ( Хз - Х1) + ( Х4 - Х1) ,
+2/ = У1 + (У2 - + (уз - У1 )л + (ш - ш)С,
+2 = 21 + ( 22 - 21){ + (2з - г1)г] + (^ - 21)(, координаты вершин ячеек в глобальной системе.
(13)
(14)
(15)
Численное решение Ор в рамках одного элемента расчетной сетки будем искать в виде разложения по пространственно-зависимому базису полиномов Фг(£() степени не выше
К, коэффициентов, зависимых от времени О^, то есть
(я{Г]ш ,ч,с, *) = (£ £ ]тт о.
(16)
В качестве базиса Ф^ (£ ,г], () выбиралась ортогональная система полиномов, образованная при помощи алгоритма ортогонализации Грамма-Шмидта из системы
1 2 2 2
1, Х, у, г, Х , у , 2 , ...
В результате подстановки решения и интегрирования по объему каждой ячейки получим следующую запись для матричного представления системы:
/г+ /гФ (^Ж + ■^ + Ж) = 0 (17)
На границах ячеек численное решение может претерпевать разрыв. Для описания процессов взаимодействия между ячейками вводится функция потока Расчет потока удобнее производить в системе координат, связанной с плоскостью грани и нормали к ней. Для построения указанной системы понадобиться:
• ввести вектор внешней нормали к грани ячейки п = (пх,пу,пг)Т,
ввести два вектора, лежащих в плоскости грани, сонаправленных с ребрами грани и
5 x, $ у, $ х ) и t — x, ~ку, ^ х )
исходящих из одной вершины. Обозначим их как 8 = ( 8х, 8У, 8Х)т и t = (£х, 1У, )т.
На основе этих векторов и будем строить оси координат. Матрица перевода Тря вектора Ор из глобальной в систему координат, связанную с гранью (обозначим вектор в этой системе как Од), будет иметь следующий вид:
Т =
Трд =
( пх ^х 2шх <§ х 2 & X 2п^ t ^ 0 0 0
п1 г> у Ьу 2Пу 5 у 2 8 yt у 2пу1/ у 0 0 0
п\ ь х 1X 2пг& г 2 $ х^ х 2п^£ ^ 0 0 0
пу пх 8 у8 х t yt X пувх + пх^у Зу^х + $хЛу ПуЪх + пт1у 0 0 0
пг пу 8 ¿8 у t ^ у пг&у + пу&г ^ у + х Пх1у + пу^х 0 0 0
8 ¿8 х х пг$х + пх$г 8 г^х + Пг1х + пх^г 0 0 0
0 0 0 0 0 0 пх 8 х t X
0 0 0 0 0 0 ПУ £у
\ 0 0 0 0 0 0 )
Перевод записывается так:
Яр = (19)
С учетом последних двух формул можно записать выражение 10 следующим образом: , С ъ.^а Г (дФ*л г, , дФ*. Л , дфк
]Г (т) ^ ^ + ]дГ (т) — ]Г (тд ^ + Ж ^ + ^) ^ = '0.
(20)
Обозначим через Я^)Фг(т) аппроксимацию решения элемента т на плоскость грани, а
через Я^Ф(тз) - аппроксимацию решения соседнего элемента т по грани ]. Обозначим абсолютное значение якобиана матрицы | А(т) | , с учетом приведенных ранее значений собственных чисел и векторов можно получить следующее выражение:
|А£ | = < 1ЛР31(К1 )-1, (21)
где |Лр5| - диагональная матрица с ненулевыми элементами (| «11, 1821 18 з |, | 18 51, | «б|, | | «в|, | «д|).
Численный поток через одну грань будем задавать следующим образом:
Рр" = \тр, (аМ + |АМ|) (Тга)-1Я!гЧм +
р — 0^РЯ ~ N \-LrsJ
1 У У (22) +2тря {а^ — а^) (т„г1^)ф(т).
Останавливаться на выводе не будем, заметим только, что приведенная формула строго верна только для случая равенства парамеров среды по обе стороны от грани ячейки. В рассматриваемом случае ее можно считать приближенной.
Подставим значение численного потока в выражение (20), получим следующее уравнение:
д+ят [ Ф*+ £ Н КТ + а^) (ПГ^) [ Ф<Т)Ф((Т)^+
сп ^ J Г(т) 2 V / J(QГ(т)).
3=1
4 1 Г
+ £ Н (■АТ — А^) (Т^ЯЧ Ф^Ф^— (23)
3=1 2 V У J(дг (ш)).
—^) /г^Ф.ЛЧ — ^/Г— ^Т) /ГдФт= 0.
Последнее записано в глобальной системе координат. Для оптимизации вычислений имеет смысл переписать уравнение в систему координат ( ), в которой три ребра тетраэдров, которым представлены элементы ячейки, совпадают с осями координат.
Лх Лу Лх = ^(24)
_д_ д1 дУ
дх
/91 дл. дс \
дх Эх Эх
д£ дц дС,
ду ду ду
д| дц д^
\ дх дх дх /
( ТтЛ
1 дц
(25)
Далее будем обозначать через Г е тетраэдр в указанной системе (в отличие от Г(т) в глобальной системе). С учетом приведенных соотношений для локальной и глобальной систем координат можно переписать (23) в следующем виде:
д г 4 1
д^^'77Гв Ф*Ф^^ + Е^1 {а{Т + (т^Г^^^-^
+ Етз4 {АТ — |А^Т)|) (ПГ1^^^ — л^И / ^— (26)
РЯ 2 \ дг | V1 ^ К^гв/ ^31 1-л'Ы -'РЯ~^д1 1-1 Jг д£
—Кя&Т^31/, дФ*Ф.^Л^ — СмЯф^^ дФ*Ф^Л^ = 0,
где
А* = А 9% + В 9% + С 9% (27) = Ард 9Х + Врд ду + Срд 9? (27)
9т] 9т] 9г1
Вр(< = Ар" 9Х + Вм Оу + См м (28)
Ср<1 = Ат 9Х + Врд 9; + Срд ^. (29)
Приведем значения величин, которые возможно рассчитать аналитически единожды и тем самым существенно увеличить производительность метода:
мк1 = I Фк Ф1<%<к]<К, (30)
•¡Ге
Кк1 = -яГ (31)
•> Ге
К1 = 1Е 9фГФ^1]^ (32)
Кк1 = I Ф1<%Кт]КС (33)
Ге 9>
Рассмотрим вычисление потоков Р—''' ,Р+1':''г'Н, упомянутых в уравнении (26). Для вычисления значения потока через грань тетраэдра производится параметризиция грани. Вводятся система координат (\т), связанная с плоскостью грани, и аналогичная система координат для смежного элемента (хт).
= I Ф
1д(ГЕ )
>к (^(х, г)) Ф1 (£Л(х, т)) 1х1т, VI <з< 4,
Р+М* = [ Фк №(х, г)) Ф1 №г)(х{11)(х, г), т(н)(х, Г))) 1х1т, (34)
Ге), V У V )
VI < х < 4, VI < г < 3,
где £ - вектор пространственных координат, индексы +,-Д,Ь используются для обозначения порядковой нумерации граней и взаимной ориентации со смежными элементами.
Далее уравнение (26) численно интегрируется методом Рунге-Кутты. В рассматриваемой задаче для ячеек, находящихся на границе расчетной области, задавались следующие потоки: Граничное условие поглощения:
РрМзогЬ = \ТР, (А^ + |АМ|)(ТГ, Г^Ф^. (35)
Условие задания свободной границы:
РрАЬзогЬ = \тт (А^ + |А(™)|)(Тг,)"1(С;Г)Фг(т) +
1 2 (36)
+-Тт(А^ - |А^|)Г„(Т^О^Ф^,
где ГГ5 = Над(-1,1,1, -1,1, -1,1,1,1).
4. Результаты численных экспериментов
Результаты представлены в виде 2Б- и 3Б-векторных полей скоростей, последовательности волновых картин. На всех изображениях показана разница между экспериментом с вертикальным разломом и без него. Таким образом, удалось избежать значительной части помех, вызванных неидельностью задания граничных условий (поглощения, свободной границы) и конечным порядком точности (имеется в виду порядок аппроксимации решения).
На рис. 2 изображен процесс распространения отклика после падения р-волны на разлом. Хорошо виден процесс появления последовательности кратных волн.
Рис. 2. Распространение отклика от вертикального разлома
Изображение читается слева направо, сверху вниз. Срез пространственной волновой картины представлен на рис. 3.
Рис. 3. Волновая картина
Рис. 4. Векторное поле скоростей
Рис. 5. Векторное поле скоростей на плоскости
На обоих изображениях отчетливо видны кратные волны, образующие на удалении эллипсоиды. Полученные результаты соответствуют аналогичным 2Б-расчетам и полностью согласуются с ожиданиями. Пространственный срез волновой картины дает хорошее представление о структуре волновых фронтов. Понимание направления распространения формируется после изучения диаграммы векторных полей. Пример визуализации векторных полей скоростей приведен на рис. 4, где слева изображена вся расчетная область, а справа приведены проекции на плоскости смежных боковых граней. Распределение скоростей, сооветствующее 2Б-срезу, изображенному на волновой картине на рис. 2, показано на рис. 5.
5. Заключение
Проведенные расчеты показывают возможность применения рассмотренного численного метода для моделирования волновых процессов, происходящих в геологических породах для получения волновых картин и их исследования. Реализованный метод на тетраэдральных сетках позволяет корректно задавать условия на контактных границах и границах области интегрирования. Это позволяет использовать его в практических задачах, для которых требование к сложности геометрии является одним из наиболее значимых. Расчетным путем показано, что метод позволяет получать отраженные от контактных границ волны и волновые картины сейсмических процессов.
Работа выполнена при финансовой поддержке РФФИ в рамках научного проекта 16-29-15097 офи-м.
Литература
1. Kaser M., Dumbser M. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes - I. The two-dimensional isotropic case with external source terms // Geophysical Journal International. 2006. V. 166, N 2. P. 855-877.
2. Dumbser M. and Kaser M. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes - II. The three-dimensional isotropic case // Geophysical Journal International. 2006. V. 167, N 1. P. 319-336.
3. Dumbser M., Kaser M. An arbitrary high-order Discontinuous Galerkin method for elastic waves on unstructured meshes - V. Local time stepping and p-adaptivity // Geophysical Journal International. 2007. V. 171, N 2. P. 695-717.
4. Kvasov I.E., Petrov, I.B. Numerical study of the anisotropy of wave responses from a fractured reservoir using the grid-characteristic method // Mathematical Models and Computer Simulations. — 2012. V. 4, N 3. P. 336-343.
5. Левянт В.Б., Петров И.Б., Муратов М.В. Численное моделирование волновых откликов от системы (кластера) субвертикальных макротрещин // Технологии сейсморазведки. 2012. № 1. С. 5-21.
6. Муратов М.В., Петров И.Б. Расчет волновых откликов от систем субвертикальных макротрещин с использованием сеточно-характеристического метода // Математическое моделирование. 2013. T. 25, № 3. C. 89-104.
7. Левянт В.Б., Петров И.Б., Голубев В.И., Муратов М.В. Численное 3D моделирование объемного волнового отклика от систем вертикальных макротрещин // Технологии сейсморазведки. 2014. № 2. С. 5-21.
8. Shewchuk J.R. A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator // https://www.cs.cmu.edu/ quake/triangle.html
9. Karypis G. METIS, Serial graph partitioninge // http://www.cs.umn.edu/ karypis/metis. 1998.
10. Миряха В.А., Санников А.В., Петров И.Б. Численное моделирование динамических процессов в твердых деформируемых телах разрывным методом Галеркина // Математическое моделирование. 2015. T. 27, № 3. С. 96-108.
11. Petrov I.B., Favorskaya A.V., Khokhlov N.I., Miryakha V.A., Sannikov A.V., Golubev V.I. Monitoring the state of the moving train by use of high performance systems and modern computation methods // Mathematical Models and Computer Simulations. 2015. V. 7, N 1. P. 51-61.
12. Hesthaven J.S., Warburton T. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications // Springer, 2008.
13. Ладонкина М.Е., Неклюдова О.А., Тишкин В.Ф. Использование разрывного метода Галеркина при решении задач газовой динамики // Математическое моделирование. 2014. T. 26, № 1. С. 17-32.
14. LeVeque, R.J. Finite Volume Methods for Hyperbolic Problems // Cambridge University Press, 2002.
15. Wilcox, Lucas C. and Stadler, Georg and Burstedde, Carsten and Ghattas, Omar A highorder discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media // Journal of Computational Physics. 2010. V. 229, N 24. P. 9373-9396.
16. Saenger, Erik H. and Shapiro, Serge a. Effective velocities in fractured media: A numerical study using the rotated staggered finite-difference grid // Geophysical Prospecting. 2002. V. 50, N 2. P. 183-194.
17. Saenger, Erik H. and Kriiger, Oliver S. and Shapiro, Serge a. Effective elastic properties of randomly fractured soils: 3D numerical experiments // Geophysical Prospecting. 2004. V. 52, N 3. P. 183-195.
References
1. Kaser M., Dumbser M. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes - I. The two-dimensional isotropic case with external source terms. Geophysical Journal International. 2006. V. 166, N 2. P. 855-877.
2. Dumbser M., Kaser M. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes - II. The three-dimensional isotropic case. Geophysical Journal International. 2006. V. 167, N 1. P. 319-336.
3. Dumbser M., Kaser M. An arbitrary high-order Discontinuous Galerkin method for elastic waves on unstructured meshes - V. Local time stepping and p-adaptivity. Geophysical Journal International. 2007. V. 171, N 2. P. 695-717.
4. Kvasov, I.E. and Petrov, I.B. Numerical study of the anisotropy of wave responses from a fractured reservoir using the grid-characteristic method. Mathematical Models and Computer Simulations. 2012. V. 4, N 3. P. 336-343.
5. Leviant V.B., Petrov I.B. and Muratov M.V. Numerical simulation of wave responses from subvertical macrofractures system. Seismic Technologies. 2012. N 1. P. 5-21. (in Russian).
6. Muratov M.V., Petrov I.B. Simulation of wave responses from subvertical macrofracture systems using grid-characteristic method. Matem. Mod. 2013. V. 25, N 3. P. 89-104. (in Russian).
7. Leviant V.B., Petrov I.B., Golubev V.I., Muratov M.V. 3D modelling of seismic responses from large vertical fractures. Seismic Technologies. 2012. N 2. P. 5-21. (in Russian).
8. Jonathan R. Shewchuk A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator. https://www.cs.cmu.edu/ quake/triangle.html
9. George Karypis METIS, Serial graph partitioninge. http://www.cs.umn.edu/ karypis/metis . 1998.
10. Miryaha V.A., Sannikov A.V., Petrov I.B. Discontinuous Galerkin method for numerical simulation of dynamic processes in solids. Matem. Mod. 2015. V. 27, N 3. P. 96-108. (in Russian).
11. Petrov I.B., Favorskaya A.V., Khokhlov N.I., Miryakha V.A., Sannikov A.V., Golubev V.I. Monitoring the state of the moving train by use of high performance systems and modern computation methods. Mathematical Models and Computer Simulations. 2015. V. 7, N 1. P. 51-61.
12. Hesthaven, J.S. and Warburton, T. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. Springer, 2008.
13. M. E. Ladonkina, O. A. Neklyudova and V. F. Tishkin Application of the RKDG method for gas dynamics problems.Matem. Mod. 2014. V. 26, N 1. P. 17-32. (in Russian).
14. LeVeque, R.J. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press. 2002. (in Russian).
15. Wilcox L.C., Stadler G., Burstedde C., Ghattas O. A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media. Journal of Computational Physics. 2010. V. 229, N 24. P. 9373-9396.
16. Saenger, Erik H. and Shapiro, Serge a. Effective velocities in fractured media: A numerical study using the rotated staggered finite-difference grid. Geophysical Prospecting. 2002. V. 50, N 2. P. 183-194.
17. Saenger E.H., Kruger O.S., Shapiro S.A. Effective elastic properties of randomly fractured soils: 3D numerical experiments. Geophysical Prospecting. 2004. V. 52, N 3. P. 183-195.
Поступила в редакцию 10.06.2016
