Содержание к диссертации
Введение
1 Механико-математическая модель среды 10
1.1 Линейно-упругая среда 10
1.2 Линейно-акустическая среда 13
1.3 Упруго-пластическая среда 16
2 Разрывный метод Галёркина 18
2.1 Одномерная задача Римана распада произвольного разрыва для линейного уравнения упругости с переменными коэффициентами 18
2.2 Численный метод решения двумерных динамических задач 23
2.3 Численный метод решения трехмерных динамических задач 30
2.4 Численный метод расчёта граничных ячеек
2.4.1 Закрепленная граница 33
2.4.2 Свободная граница 34
2.4.3 Заданная внешняя скорость или сила 34
2.4.4 Условие поглощения 34
2.5 Условия на контактной границе 35
2.5.1 Полное слипание 35
2.5.2 Свободное скольжение 35
2.5.3 Сухое трение 36
2.5.4 Численное решение задачи контактного разрыва 2.6 Задание точечных источников и объёмных сил 41
2.7 Алгоритм выбора кратного шага по времени для многостадийных интеграторов 2.7.1 Описание алгоритма 43
2.7.2 Особенности реализации алгоритма 43
2.8 Регуляризация численного метода 46
2.8.1 Выбор системы базисных полиномов 47
2.9 Учёт разрушения материалов 48
2.10 Особенности программной реализации 50
2.11 Реализация метода для высокопроизводительных многопроцессорных вычислительных систем 52
2.12 Верификационные расчёты 56
2.12.1 Задача Лэмба с точечным источником 56
2.12.1.1 2D задача Лэмба 56
2.12.1.2 3D задача Лэмба с точечным источником
2.12.2 Расчёт волновых процессов в упругом образце с жёстким включением 61
2.12.3 Распространение упругих волн в геторогенных средах: контакт воздух-сталь 62
2.12.4 Распространение волн в многослойной геологической среде. Сравнение результатов различных численных методов. 64
2.12.4.1 Постановка задачи
2.12.4.2 Сравнение результатов расчётов
2.12.4.3 Заключение
2.12.5 Столкновение двух пластин
2.12.6 Столкновение шара и пластины
2.12.7 Моделирование плавучести сферического тела
2.12.8 Решение квази-статических задач
2.12.9 Расчёт с кратными шагами по времени
3 Совместное решение систем уравнений теории упругости и акустики
3.1 Численное решение задач сейсморазведки в шельфовой зоне
3.1.1 Распространение волн в системе упругая среда-флюид
3.1.2 Сейсморазведка при наличии ледяного покрова
3.1.3 Волновой отклик от геологической структуры 3.1.4 Численное моделирование возмущений от подводных объектов
3.2 Оценка влияния степени раскрытости трещины на сейсмический отклик .
3.2.1 Введение
3.2.2 Методика численного моделирования
3.2.3 Верификация метода
3.2.4 Характеристики моделей сплошных сред
3.2.5 Одномерный анализ применимости модели бесконечно тонкой трещины (БТТ)
3.2.6 Описание параметров численного метода и расчётных сеток
3.2.7 Двумерный анализ - сравнение сейсмограмм полученных от БТТ с сейсмограммами от трещин конечной раскрытости
3.2.8 Сопоставление сейсмограмм волнового отклика от кластера трещин
3.2.9 Анализ фазовых переходов у обменной дифрагированной волны на записи горизонтальной компоненты
3.2.10 Выводы 66 66
Взаимодействие массивных ледовых образований со стойкой нефтедобываю щей платформы 104
4.1 Реология льда 105
4.2 Особенности численного решения задачи 108
4.3 Результаты численных экспериментов 111
5 Мониторинг состояния рельсов 113
5.1 Ультразвуковая дефектоскопия рельсов 113
5.1.1 Распространение упругих волн в рельсе. 2D приближение. 114
5.1.2 Численное моделирование волновой картины в поврежденном рельсе 115
5.1.3 Распространение ультразвуковых волн в профиле рельса с горизонтальным расслоением головки 115
5.1.4 Ультразвуковая диагностика рельса на наличие дефекта в головке на начальном этапе формирования 116
5.2 Определение остаточной прочности рельса 117
Заключение 125
Список иллюстраций
- Линейно-акустическая среда
- Заданная внешняя скорость или сила
- Распространение волн в системе упругая среда-флюид
- Особенности численного решения задачи
Линейно-акустическая среда
Существует несколько основных классов численных методов для решения систем уравнений в частных производных: конечно-разностные методы (finite-difference methods, FDM); конечно-элементные методы (finite-element methods, FEM); конечно-объёмные методы (finite-volume methods, FVM).
Для решения гиперболических систем уравнений наиболее подходят конечно-разностные и конечно-объёмные методы. В 1971 году был предложен новый класс конечно-элементных методов [70], а именно разрывный метод Галёркина (discontinuous Galerkin finite element method, DG-FEM), сочетающий в себе ряд положительных свойств для решения гиперболических систем [71]: консервативность схемы; устойчивость схемы при больших скачках решения, которые присутствуют и в аналитических решениях, т.е. не являются по своей природе нефизичными осцилляциями. гибкость в выборе базиса, по которому раскладывается решение в ячейках. -адаптивность (локально измельчаем сетку, сохраняя порядок полиномов (), локально увеличиваем порядок полиномов ()).
Метод уже успел зарекомендовать себя, и применялся также для решения нелинейных гиперболических задач, параболических и эллиптических. Получены оценки погрешности, исследована сходимость метода как теоретически, так и на практике [72].
В книге [73] представлена таблица (см. таблицу 2.1), наглядно сопоставляющая основные численные методы решения систем уравнений в частных производных, из которой можно сделать вывод, что разрывный метод Галёркина по совокупности свойств наиболее подходит для решения, в частности, гиперболических задач, исследуемых в данной работе.
При построении численного метода нам предстоит построить функцию численного потока, которая тесно связана с решением задачи Римана [68]. Таблица 2.1: Сводная таблица свойств основных методов дискретизации систем уравнений в частных производных. /- успех, X - определенные трудности, (/) - необходим ряд дополнительных модификаций. Сложная Высокий порядок, Явная форма К онсервативност геометрия hp-адаптивность записи FDM X
Рассмотрим квази-одномерную задачу Римана о распаде произвольного начального разрыва или, для краткости просто задачу Римана для линейной гиперболической системы уравнений упругости, которая является по сути задачей Коши:
Стоит отметить, что мы рассматриваем общий случай, когда матрица зависит от координаты. Матрица зависит от параметров среды и для двумерного случая определяется формулой (1.8), а для трехмерного - формулой (1.14). Мы назвали задачу квази-одномерной, т.к. нормаль к границе раздела сред n и направление распространения волн сонаправлены с координатной осью и поэтому наличие оставшихся слагаемых системы уравнений упругости мы не учитываем. Далее приведём подробно выкладки для двумерного случая, а для трёхмерного - только конечный результат ввиду относительной громоздкости и аналогичности рассуждений. коэффициентами
Вычислительная область (см. рис. 2.1) разбита на участки, решение в которых постоянно. Границами этих участков являются характеристики системы уравнений 2.1: = , которые выпущены из начала координат. На каждой характеристике происходит скачок решения. Величину скачка можно определить как из физических соображений, так и из полностью математических. С точки зрения физики величины скачков на разрыве можно найти, применяя законы сохранения [60,74], которые формулируются аналогично условиям Ранкина-Гюгонио для газодинамики. Для линейной системы дифференциальных уравнений это условие сводится к тому что при переходе через характеристику решение изменяется на величину, пропорциональную соответствующему правому собственному вектору (столбцу).
Покажем, как это условие получается формально с математической точки зрения. Мы рассматриваем гиперболическую систему уравнений 2.1, следовательно имеющую полный набор действительных собственных значений и соответствующих собственных векторов (правых и левых). Пока будем считать, что А = А+ = А. Позже это ограничение будет снято. Тогда для матрицы А справедливо: где L - матрица из левых собственных векторов (строк), R - матрица из правых собственных векторов (столбцов), Л - диагональная матрица из собственных значений Si. Без ограничения общности можно считать, что собственные значения Si упорядочены по неубыванию, а собственные вектора нормированы. Из 2.4 и 2.5 видно, что L = ВТ1.
Матрица Л - диагональная, поэтому 2.6 распадается на систему независимых уравнений, каждое из которых является уравнением переноса. Решением этих уравнений является система разрывов Wi, распространяющихся со скоростями Si от точки х = 0. Следовательно начальный скачок решений и+ - и" ввиду 2.8 раскладывается по правым собственным векторам: где а - столбец неизвестных коэффициентов. Т.к. собственные вектора линейно независимы (т.е. R l существует), то система уравнений 2.9 разрешима относительно неизвестных ац и имеет единственное решение. Уравнение 2.9 можно переписать в виде:
Теперь вернёмся к общему виду уравнения 2.1, когда параметры сред при х 0 и при х 0 различны. Для этого достаточно составить матрицу R так: т.е. собственные вектора, которым соответствуют неположительные собственные значения Si 0 берём из А , а положительные Si 0 из А+. Нулевому собственному значению соответствует собственный вектор, не зависящий от параметров среды (см. 1.12), что избавляет от неопределенности выбора параметров среды (— или +) для этого вектора (если бы мы оперировали с матрицей L из левых собственных векторов, то присутствовала бы неопределённость для нулевого собственного значения). Стоит отметить, что 2.1 записано в неконсервативной форме, те. производная по координате берётся не от А(х)и(х), а только от и(х). В случае же консервативной записи уравнения 2.1 запись его в переменных Римана имела бы следующий вид:
Здесь учли, что L(x)ur = (L(x)u) — Lr(x)u, где Lr(x) = а и Віх) = Lr(x)R(x) (недиаго-нальная матрица в общем случае). Т.е. уравнения переноса для Wi теперь не являются независимыми. В [58] как раз показывается, что используя приём с составлением комбинированной матрицы R (см. формулу 2.11) для случая кусочно-постоянных параметров среды (см. условие 2.3), можем ограничиться решением задачи Римана для гетерогенного пространства для х = 0 и не учитывать дополнительное слагаемое в правой части, т.к. Ьх(х) = 0 при \х\ 0. Важно понимать, что это не значит, мы решаем задачу распада Римана и для консервативной формы уравнения
Заданная внешняя скорость или сила
Здесь предполагается суммирование по повторяющимся индексам. Такой способ представления решения соответствует всем конечно-элементным методам, но в нашем случае на границах элементов допускаются разрывы в численном решении. Фактически область можно было бы разбить на произвольные симплексы (четырёхугольники, например) при условии задания соответствующего базиса на них, т.е. ограничений на форму ячейки нет. На самом деле нет и строгого ограничения на тип используемых базисных функций. Например, в [77,78] предлагается раскладывать решение в сумму плоских волн. Выбор типа базиса определяется видом решения конкретной задачи, и может меняться как от ячейки к ячейке, так и во времени. Такая гибкость достигается благодаря возможности решения претерпевать разрывы на границе ячеек, в отличие от традиционных конечно-элементных методов [79]. В работе [6,76] предложена модификация численной схемы, когда внутри каждого треугольника Т якобианы Ард и Врд являются функциями координаты и также раскладываются по системе базисных полиномов Ф{Г\х,у). Целесообразность такого приёма несколько сомнительна, т.к. она значительно увеличивает затраты на вычислительные ресурсы. Также разложение параметров среды по полиномам степени выше 1 может вызвать появление нефизичных осцилляций, т.е. привести, например, к отрицательным значениям физических параметров сред А,//,р. Также ввиду того, что используются неструктурированные расчётные сетки можно всегда построить расчётную сетку в соответствии с геометрией неоднородностей. Но всё же такая модификация может быть существенна в задачах с непрерывно изменяющимися параметрами среды, градиент которых велик, при счёте на относительно грубой расчётной сетке.
Подставив в систему уравнений 1.7 численную аппроксимацию (и)р получим: где в правой части стоит не 0, некоторая функция невязки Д(т)(:г,г/,). Большое количество численных методов строятся из условия, что функция невязки должна быть ортогональна некоторому заданному набору функций. В зависимости от этого набора выделяют метод кол-локаций, метод Галёркина, метод подобластей и другие. Если в качестве такого набора функций выбрать полиномиальный базис Фг(т)(ж,у), по которому раскладывали решение в каждой ячейке, получим семейсво методов Галёркина.
Выпишем явно условие ортогональности невязки базисным функциям, т.е. умножим 2.31 на базисную функцию Ф т) и проинтегрируем по треугольнику Т :
Существует общая формула интегрирования по частям для n-мерного пространства, но, в частности, для 2D случая она эквивалентна формуле Грина, а в 3D - формуле Гаусса-Остроградского [80].
Ввиду разрывности матриц Ард, Врд и решения Uh на границе ячейки в общем случае напрямую произвести интегрирование по частям не удаётся. Проблема решается с помощью введения функции численного потока F через j-е ребро треугольника в глобальной системе координат OXY. Известно большое число возможных вариантов задания численного потока [81-86]. В нашем случае наиболее естественно [81,82,85] выбрать набегающий поток (upwind flux), который отражает волновую природу решаемой системы уравнений. Нахождение набегающего потока сводится к решению задачи распада произвольного разрыва Римана [68]. В случае рассматриваемого нами изотропного пространства многомерную задачу Римана можно приближённо свести к одномерной (т.н. x-split Riemann problem [83]). В случае линейных систем гиперболических уравнений одномерная задача Римана имеет аналитическое решение [81,82]. Для случая систем уравнений упругости и акустики задача Римана решена в [62] и рассмотрена нами подробно в предыдущем разделе 2.1. соответствует системе координат OX Y . Тензора 2-го ранга, и тензор напряжения, в частности, преобразуется в матричной нотации по правилу о = ( хх ху ) = Т2апТ [87]. В тензорной же нотации преобразование можно записать так: Тч- = tutjmT(m, где числа ч- -элементы матрицы поворота Т2. Символьно перемножив три матрицы, получаем выражение для элементов матрицы Тг такой, что (сгхх,ауу, оху)т = Ti(a xx, а , о ху)т.
Ещё раз подчеркнём, что пространство изотропно, поэтому мы можем приближённо вычислить поток в системе координат OX Y , т.е. решить ж-split Riemann problem, а затем повернуть систему координат обратно.
В предыдущем разделе 2.1 мы обозначили, что возможно составить много выражений для численного потока, используя аналитическое решение задачи Римана, но остановились именно на 2.27. Только выражение 2.27 удовлетворяет условию согласованности численного потока. Допустим мы знаем точное решение uexact исходной линейной системы уравнений упругости 1.5. Тогда подставляя в выражение 2.29 и" = и+ = uexact (контактные условия полного слипания 2.14 уже в нём учтены) получаем F = (m)uexact. Далее проводим фактически сворачивание формулы интегрирования по частям для 2.33 и получаем тождество. Выполнение условия согласованности потока крайне важно. Если использовать вместо 2.27, например, А и} или А т +2А 3 гг1-, то схема оказывается неустойчивой. Это особенно проявляется при контакте сред с сильно отличными импедансами (см. выражение 2.16) (на 2 и более порядков).
При вычислении интеграла f kF dS вдоль ребра, т.е. при подстановке в 2.35 од V j них и тех же координат точек на этом ребре, выражения для и и и+ перейдут в пределы решения с внутренней и внешней сторон j-го ребра треугольника. Т.е. на поток влияет не всё решение в ячейках T(m) иTW, а только значения вдоль соответствующего ребра. Объединяя выражения 2.29, 2.35, 2.36 получим:
Распространение волн в системе упругая среда-флюид
Неотражающее граничное условие - некая абстракция, подразумевающая, что никакие волны не входят снаружи в вычислительную область, а выходящие волны не отражаются от таких границ. Создание качественных неотражающих условий - отдельное направление исследования, поэтому в работе используется простой, естественно вписывающийся в численный метод подход, дающий удовлетворительный результат. Выражение для численного потока 2.37 состоит из двух слагаемых: вытекающий и втекающий поток. Если мы занулим втекающий поток, то получим некий вариант неотражающего граничного условия. Такой метод неотражающего условия хорошо работает только при нормальном падении волн на границу. При касательно падении или падении на углы расчётной области, граничное условие ”фонит”. Более сложным, но и более качественным неотражающим условием является PML (Perfectly Matched Layers) [8,90,91]. Поток для условия поглощения можно также записать в виде 2.63, задав rabsorb = diag{0,0 ... ,0
Учёт всех контактных условий одинаков для всех типов условий: при решении квазиодномерной задачи Римана (см. раздел 2.1) составляется линейная система уравнений, связывающая скачки решений при переходе через ненулевые характеристики, которая для нулевой характеристи дополняется физическими соотношениями на контакте. Решение задачи Римана используется уже для вычисления численного потока (см. раздел 2.1, 2.2).
В разделе 2.1 рассуждения как раз проводились для случая полного слипания, когда на контакте сред в точке = 0 должны выполняться соотношения 2.14, выражающие непрерывность скорости v и плотности сил ( n) на контакте.
В случае свободного 2D скольжения соотношения на скачках дополняются условиями непрерывности нормальной составляющей скорости и отсутствием сдвиговых напряжений на контакте:
Для случая полного слипания, когда выполняется условие непрерывности всех компонент скорости и тензора напряжения на контакте, коэффициент 0 совпадает со случаем скольжения, т.е. продольные волны распространяются в этих случаях одинаково. Для лучшего понимания физического смысла условия скольжения можно считать, что коэффициент 1 для скольжения (см. формулу 2.70) получается из соответствующего коэффициента для слипания (см. формулу 2.19) путём зануления продольной скорости и сдвигового напряжения в граничащей ячейке, т.е. + = + = 0, что равносильно контакту с акустической средой, в которой нет поперечных волн.
Хотелось бы записать конечные выражения для соответствующих потоков для всех типов контактов единообразно в матричном виде, как это получилось сделать для граничных условий. К сожалению, компактно и аккуратно численные потоки выписываются только для полного слипания (см. выражение 2.29), либо для скольжения при контакте сред с одинаковыми параметрами, поэтому здеcь приводиться не будут.
Скольжение начинается сдвиговое напряжение превышает некоторый порог, который пропорционален нормальному напряжению. Во время скольжения скорость касающихся берегов и сдвиговое напряжение имеют противоположные знаки. Все эти явления можно описать математически:
В данной работе предполагаем, что коэффициент трения р/ - константа для заданной пары контактирующих материалов. Известны более сложные модели зависимости коэффициента трения от локальных параметров среды, в частности, для геофизики [92,93]. В [94] приводятся варианты законов трения для случаев контакта резины и полимеров, бетонных конструкций и почвы и т.д. - проверка контакта двух берегов (напомним, что положительное нормальное напряжении соответствует растяжению, отрицательное - сжатию). При положительном значении TS 0 контакта между берегами нет, и получаем случай двух граничных ячеек. Важно подчеркнуть, что ”корректировке” в случае проскальзывания подлежат также величины и и в отличие от приблежения твёрдого недеформируемого тела, где этого делать не нужно. Отсутствие корректировки приведёт к неверной волновой картине в случае деформируемых тел. Алгоритм 1 приведён для 3D случая. Для 2D случая необходимо просто положить все величины, отвечающие за z-компоненты нулями.
Учёт наличия статического решения и0 зависит от решаемой задачи: в некоторых постановках задачи устойчивости строительных конструкций или сейсмического отклика от неод-нородностей, находящихся на значительной глубине, при амплитудах динамического воздействия сравнимых со статическими нагрузками учёт и0 необходим.
В разрывном методе Галёркина нам необходимо вычислять интегралы по ребрам (граням) от численного потока F (см. выражение 2.37) через ребра (грани). В случае линейных условий контакта (полное слипание или свободное скольжение) интеграл от потока можно вычислить, используя предрасчитанные интегралы от базисных функций и/или их производных. В случае нелинейного контактного условия сухого трения интегралы от потока F J необходимо вычислять численно, используя, например, квадратурные формулы интегрирования Гаусса-Лежандра-Лобатто.
После введения нелинейного контактного условия система уравнений упругости перестаёт быть линейным, и возможно появление нефизичных осцилляций. В работах [10,12,95] отмечается, хоть метод и становится нелинейным, но нефизичных осцилляций в случае сухого трения не наблюдается, даже при более сложной зависимости /if от локальных параметров. И объясняется тем, что используется аналитическое решение задачи Римана при вычислении потока.
Отдельно остановимся на физическом смысле условия т& 0, отвечающего на вопрос, находятся ли в контакте две среды. Предположим, что статическое решение ахх = 0, тогда условие контакта сведётся к существует аналитическое решение, связывающую давление и амплитуду скорости частиц бегущей волны: р = — а = cpv. Волна сжатия ( т 0, р 0 - давление) с таким соотношением между скоростью и давлением пройдёт через такой контакт, не отразившись. Случай прохождения плоской волны можно считать пограничным, т.к. пластины ещё ”не давят” друг на друга, но контакт между пластинами есть и волны сжатия проходят из одной среды в другую. Тогда выражение 2.84 при равенстве его нулю можно интерпретировать следующим образом: а%.х - это суперпозиция двух волн, набегающих на контакт с разных сторон. Если, например, а х или ахх уменьшить, то это будет соответствовать тому, что давление с данной стороны на другую увеличилось. В этом случае контакт между средами также будет выполняться. Следовательно для контакта должно выполняться соответствующее неравенство. Условие 2.84 симметричное, поэтому если есть контакт 1-го тела со 2-м, то и 2-го с 1-м тоже будет.
Особенности численного решения задачи
Анализу сейсмических эффектов, связанных с вторичными волнами от заполненных флюидом трещин конечной толщины (раскрытости), в последние годы посвящены ряд работ. Наиболее полное исследование на базе численного моделирования волн типа guided wave, распространяющихся вдоль трещин с конечной раскрытостью, их отражения и рассеяния от концов трещин приведено в [168]. Образование при падении волнового фронта на трещиноватый пласт насыщенный флюидом волн Крауклиса (их также иногда называют Stoneley guided waves), колеблющихся внутри трещины и при определенных условиях рассеивающихся наружу, рассмотрено в работах [63,169,170]. Однако все эти эффекты, связанные с трещинами конечной раскрытости, рассматривались на очень ограниченных фрагментах среды, не касались основных волн и не представляли полной картины волновых откликов на дневной поверхности. Работ с оценкой влияния раскрытости субвертикальной трещины на отклик от системы таких трещин авторы не встречали вовсе.
Разрывный метод Галёркина, сочетающий в себе высокий порядок сходимости по пространственным координатам и времени, а также возможность задания различных граничных условий на поверхности трещин, позволяют решать широкий круг задач сейсморазведки. Важной моделью трещеноватости является абстрактная модель Бесконечно Тонкой флюидо-насыщенной Трещины (БТТ) [171–173] с аспектным отношением (отношение толщины трещины к её высоте, ) в пределе равным 0. Для широкого применения модели БТТ необходимо было доказать правомерность её использования на основе сопоставления с трещинами конечной раскрытости, характерными для реальных сред.
В реальных условиях геологоразведочных работ на нефть и газ сведения о параметрах трещин в продуктивных породах получают по результатам анализа керна, шлифов и данных сканеров FMI, FMS , дипольного акустического каротажа X-MAK., а также изучения обнажений. По этим данным, приведенным в ряде обобщающих работ [174–176], известно, что толщина (раскрытость) макротрещин (трещиноватых коридоров), сопоставимых с длиной сейсмической волны, колеблется в пределах от первых миллиметров до их десятых долей. В среднем их раскрытость может быть принята равной 1мм. Тогда для макротрещин аспектное отношение будет 1 : 100000. Мезотрещины - массовые трещины, которыми пронизаны про-луктивные пласты, имеют высоту в первые метры и раскрытость от долей миллиметров до десятков микрон. Их аспектное отношение в среднем может быть принято равным 1 : 10000 и в отдельных случаях 1 : 1000. В какой мере отклик от модели БТТ соответствует откликам от трещин с указанными выше параметрами раскрытости, определяло первую задачу данной работы.
Для оценки применимости модели БТТ, в частности, для расчета откликов волн адекватных откликам от трещин с реальной конечной раскрытостью проведены расчеты, с использованием программно-вычислительного комплекса на основе Разрывного метода Галер-кина [35]. Метод позволяет рассчитывать волновые поля как для трещин с конечной раскры-тостью, так и с нулевой (БТТ).
Для изучения характера распространения упругих волн, вызванных прохождением падающего (продольного) волнового фронта через субвертикальную макро- или мегатрещину, был использован программно-вычислительный комплекс [35], основанный на разрывном методе Галёркина на треугольных неравномерных сетках [1]. Метод имеет ряд положительных особенностей [81], основной из которых для авторов является возможность использования высокого порядка сходимости по пространственным координатам и времени из-за волновой специфики решаемой задачи. В работе в качестве системы базисных полиномов использовались полиномы Лагранжа и Дубинера 4-6 порядков. Интегрирование по времени проводилось с помощью интегратора Дорманда-Принца 5-го порядка, либо Рунге-Кутта-Фельберга 4-го порядка с адаптивным шагом по времени.
Для среды использовалась модель идеального изотропного линейно-упругого материала [57] (см. раздел 1.1). Для описания поведения флюида в случае конечной раскрытости использовалась двумерная модель линейно-акустического материала [59]:
Рассмотрены условия контакта массивной породы с заполненными жидкостью трещинами. В модели бесконечно тонкой флюидонасыщенной трещины с нулевой раскрытостью контакт с жидкостью задавался условием свободного скольжения на берегах трещины. В случае ненулевой раскрытости решалась совместная система линейной упругости и акустики [3, 31,35,49, 62].
Для верификации расчётной программы, основанной на разрывном методе Галёркина [35] был проведён расчёт, описанный в статье [177]. В ней с использованием той же системы уравнений и конечно-разностного численной схемы с неравномерной сеткой проведен расчет отклика от вертикальной трещины в модели тонкого вертикального пласта (TLM, thin-layer model) и модели Шоенберга [178], инициированного точечным источником p-волны (монополем).
На рис. 3.7 представлено сопоставление полученных сейсмограмм: красной сплошной линией показано расчетное решение, пунктирной – полученное в [177]. Небольшое несоответствие порядка нескольких процентов можно отчасти объяснить погрешностью оцифровки растровой картины сейсмограммы для сопоставления. На рис. 3.8 приведена волновая картина горизонтальной и вертикальной составляющей скорости смещения частиц в момент времени 2.8 мс, которая визуально полностью совпадает с [177]. Расчётная сетка состояла из 1.5 105 ячеек, размер ребра которых менялся от 2 мм у трещины до 36 см на краях расчётной области. На ширину трещины приходилось 2 ячейки, зато использовались полиномы Дубинера 6-го порядка, что эквивалентно 11 точкам на трещину для конечно-разностного метода (не 12, т.к. точки фактически дублируются на границе ячеек из-за особенностей численного метода), использованного в [177], где авторы используют всего 5 точек на толщину трещины. Расчётная сетка изотропна, в отличие от оригинального расчёта, где вертикальный от минимального горизонтального шагов сетки различались в 60 раз.