1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 63, № 10, 2023

Журнал вычислительной математики и математической физики, 2023, T. 63, № 10, стр. 1706-1720

Моделирование распространения динамических возмущений, в пористых средах сеточно-характеристическим методом с явным выделением неоднородностей

И. А. Митьковец 1, Н. И. Хохлов 1, *

1 МФТИ (НИУ)
141701 М.о., Долгопрудный, Институтский пер., 9, Россия

* E-mail: khokhlov.ni@mipt.ru

Поступила в редакцию 12.05.2023
После доработки 12.05.2023
Принята к публикации 26.06.2023

Полный текст (PDF)

Аннотация

Рассматривается вопрос численного моделирования распространения волновых возмущений в гетерогенных средах с наличием пористых включений, а также вопрос явного выделения пористых неоднородностей. В качестве подхода для явного выделения неоднородностей предложен метод наложенных сеток. Для численного решения возникающих систем дифференциальных уравнений в частных производных применяется сеточно-характеристический метод. Рассмотрены особенности предложенного метода, проведена верификация предложенных алгоритмов, приводится серия тестовых расчетов. Библ. 48. Фиг. 9.

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

1. ВВЕДЕНИЕ

Одним из самых старых и распространенных методов поиска и разведки месторождений углеводородов на текущий момент, по-видимому, является сейсмическая разведка. В то время, как основной ее задачей является определение структуры подповерхностного пространства по данным наземных (и, возможно, скважинных) наблюдений, важной составляющей решения данной задачи является описание процесса распространения сейсмических волн от источника вглубь геологической среды. Решение прямых задач сейсморазведки может быть использовано, например, для уточнения запасов углеводородов в районах с большим числом эксплуатируемых скважин, построения синтетических сейсмограмм для нужд региональной сейсморазведки, уточнения присутствия углеводородов и их характеристик в известной геологической структуре. Значительное развитие в данном направлении возможно с использованием численного моделирования распространения сейсмических волн в реалистичных моделях геологического массива, поскольку реализуется возможность задания произвольной внутренней структуры среды с последующей оценкой синтетического сигнала-отклика.

Настоящая работа посвящена исследованию пространственных динамических процессов, протекающих в геологических средах, содержащих пористые включения, в процессе сейсмической разведки. Геологические породы, содержащие пористые и трещиноватые включения, являются одним из источников залежей углеводородов. Более того, такие структуры с низкой пористостью и проницаемостью используются для контролируемого хранения и миграции газа. Относительно высокая пористость и развитая сеть микротрещин имеют важное значение для добычи природного газа в плотных газоносных песчаниковых коллекторах, являющихся одним из основных факторов обогащения газа (см. [1]). В то же время развитие природных микротрещин способствует формированию сети пор-трещин при гидроразрыве пласта, что является одним из определяющих факторов добычи углеводородов. При решении обратной задачи, опираясь на сейсмические данные и имея таким образом модель геологического устройства под поверхностью, можно определять оптимальные места бурения скважин. Вследствие чего для нефте- и газодобычи необходимы модели, позволяющие максимально точно описывать физические свойства подобных структур. В частности, важно создание моделей, описывающих сейсмические характеристики пор и микротрещин (см. [1]).

При распространении упругих волн через среду, имеющую такие неоднородности, как поры или трещины, будет происходит ослабление возмущения. Помимо таких факторов, как потеря плотности энергии, вызванных расширением фронта волны и трения зернистых материалов, это связано с частичным отражением упругой волны от неоднородностей (см. [2]). Величина последнего эффекта зависит от пористости среды, формы пор, трещин, их взаимного расположения и материала, которым они заполнены. Все эти параметры могут варьироваться в реальных геологических породах, поэтому они должны учитываться в реалистичной модели среды (см. [3]). Важно, чтобы разработанные математические модели и численные методы учитывали эти факторы.

Существуют теоретические модели подобных сред. Классической является модель Гасcмана (см. [4, 5]), в которой устанавливается связь между упругими параметрами насыщенной жидкости или газом пористой среды. В модели Гассмана предполагается, что мы можем описывать упругие свойства пористой среды аналогично однородному изотропному материалу. Для этого в модели предполагается, что порода состоит из твердого каркаса (скелета), по которому равномерно распределен флюид в жидкой или газообразной фазах. Основываясь на этом предположении, производится расчет “эффективных” параметров среды, описывающих упругие свойства материала. Предположения, заложенные в модель Гассмана, накладывают определенные ограничения на ее применение. Основное ограничение связано с частотой распространяющегося сигнала. При достаточно высоких частотах модель Гассмана плохо описывает упругие характеристики двухфазных пород из-за отсутствия учета влияния движения флюида относительно твердого каркаса. Физическое влияние пор учитывается в твердо-жидкостной модели Био, описывающей пористые среды, в которых жидкая и твердая фазы имеют схожую плотность (см. [6]). Как и в модели Гассмана, в модели Био считается, что поры равномерно распределены по всему объему среды, но при этом учитывается частота проходящей волны и в движение вязкой жидкости в порах. На основе этой модели были разработаны другие, учитывающие свойства жидкости, проницаемости и углов падения, такие как модель Squirt и BISQ (см. [7–10]), частично насыщенная модель (см. [11, 12]). Была разработана обобщенная модель Био для анизотропных вязко упругих сред, поскольку в природе существуют породы, обладающие анизотропией по отношению к направлению упругой волны (см. [13]). Также в ряде исследований проводилось изучение затухания продольных волн в частично насыщенных породах (см. [14–16]). Микроструктурный подход к аналитическому изучению затухания и дисперсии волн, проходящих через пористую структуру, проводился в работах [17–19]. Для описания среды в них использовалось представление пор в виде эллипса с заданным соотношением полуосей и пористость (доля объема занимаемого полостями пор).

На основе нелинейной теории фильтрации (см. [20]) была предложена линеаризованная модель Доровского (см. [21]). В этой модели используются все скорости распространения упругих волн и требуется всего три параметра для описания пористой среды, против четырех в модели Био. Однако она не позволяет точно описывать волны Стоунли при различных частотах (см. [22]).

Тем не менее неоднородная структура пор и микротрещин может приводить к сложным зависимостям между скоростью прохождения упругой волны и пористостью, которые не поддаются классическим формулам (см. [23]). Степенная зависимость скорости прохождения упругой волны от соотношения полуосей эллипса, описывающего поры, была продемонстрирована в [24]. При этом для поверхностных волн такая зависимость не была выявлена. Чувствительность продольной скорости упругой волны к соотношению сторон исследуется в [25]. Отмечается различие скорости для породы с иглообразными и сферическими порами при одинаковых макропараметрах. Также отмечается влияние насыщающего флюида (газа или жидкости). Важность учета распределения пор по размерам и толщины их стенок была продемонстрирована с помощью моделирования на микроуровне в работе [26].

Из вышеизложенного следует, что развитие более точных математических моделей пористых сред, в том числе с явным выделением неоднородностей, является актуальной задачей. В настоящей работе предлагается новый подход для явного выделения пористых неоднородностей, основанный на наложенных или химерных сетках. Ранее данный подход уже применялся для выделения геологических трещиноватых неоднородностей. Так, в работе [27] рассматривался вопрос применения повернутых наложенных сеток для выделения отдельных трещин. В [28] применялся подход на искривленных наложенных сетках, который позволяет более точно выделять трещиноватые неоднородности. Также в работе [29] описано применение метода наложенных сеток для явного выделения топографии земной поверхности.

В качестве численного метода для решения возникающей системы дифференциальных уравнений в настоящей работе используется сеточно-характеристический метод (см. [30]), который изначально разрабатывался для решения задач газовой динамики. Одним из основоположников данного метода был А.С. Холодов (см. [31, 32]). Затем сеточно-характеристический метод был обобщен для решения задач упругости. Данное направление начал развивать И.Б. Петров. Первые работы на эту тему можно найти в [33, 34]. И.Б. Петров внес большой вклад в развитие сеточно-характеристического метода для данного направления. Позже он был обобщен для решения задач сейсмики (см. [35]), явного выделения геологических неоднородностей (см. [36]), задач сейсмостойкости (см. [37]) и других динамических задач механики деформируемого твердого тела. В настоящей работе также применяется сеточно-характеристический метод.

2. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ И ЧИСЛЕННЫЙ МЕТОД

2.1. Математическая модель

Для решения обобщенной задачи моделирования распространения сейсмической волны в грунте будем использовать линейно-упругую и изотропную модель среды (см. [38–41]). В тензорной форме она имеет вид (см. [30])

(1)
$\rho {{{\mathbf{v}}}_{t}}{\mathbf{ = }}{{\left( {\nabla \cdot {\mathbf{T}}} \right)}^{{\text{T}}}},$
(2)
${\mathbf{T}} = \lambda \left( {\nabla \cdot {\mathbf{v}}} \right){\mathbf{I}} + \mu \left( {\nabla \otimes {\mathbf{v}} + {{{\left( {\nabla \otimes {\mathbf{v}}} \right)}}^{{\text{T}}}}} \right),$
где $\rho $ – плотность среды, ${\mathbf{v}}$ – скорость перемещения среды, ${\mathbf{T}}$ – тензор напряжений Коши, $\lambda $ и $\mu $ – параметры Ламе, характеристики упругих деформаций, ${\mathbf{I}}$ – единичный тензор, $ \otimes $ – оператор тензорного произведения ${{\left( {\nabla \otimes {\mathbf{v}}} \right)}_{{i,j}}} = {{\nabla }_{i}}{{{\mathbf{v}}}_{j}}$.

Для задания пористых неоднородностей предлагается их явное выделение. Пористое включение представляется в виде пустой или заполненной поры. В случае, если пора пустая, на границе области задается условие свободной границы. Если пора заполнена другим материалом, задается контактное условие на границе раздела поверхностей. При описании свободных граничных условий свободной поверхности мы полагаем проекцию тензора напряженности на нормаль к границе, равной нулю в граничных узлах, таким образом, должно выполняться

(3)
${\mathbf{T}} \cdot {\mathbf{n}} = 0,$
где n – вектор единичной длины, перпендикулярный границе и ориентированный против области моделирования.

В случае заполненной поры необходимо контактное условие на границе раздела сред. Для описания передачи возмущения между сетками с различными характеристиками среды нами был реализован механизм контакта сеток. Реализация данного механизма предполагает, что граничные узлы таких сеток в месте их контакта геометрически совпадают. В таком случае мы полагаем равенство нормальных составляющих тензоров напряженности в совпадающих узлах таких сеток и равенство скорости контактных сред. Эти условия можно записать в виде соотношений

(4)
${{{\mathbf{v}}}_{0}} = {{{\mathbf{v}}}_{1}},$
(5)
${{{\mathbf{f}}}_{0}} = - {{{\mathbf{f}}}_{1}},$
(6)
${{{\mathbf{f}}}_{0}} = {{{\mathbf{T}}}_{0}} \cdot {\mathbf{n}},\quad {{{\mathbf{f}}}_{1}} = {{{\mathbf{T}}}_{1}} \cdot ( - {\mathbf{n}}),$
где индексы $0$ и $1$ обозначают различные контактные сетки; fi – сила, действующая со стороны области, описываемой $i$-й сеткой; ${{{\vec {v}}}_{i}}$ – скорость среды на рассматриваемой границе. Более подробное описание данных граничных и контактных условий приводится в [30].

2.2. Численный метод

Сеточно-характеристический метод использует характерные свойства систем гиперболических уравнений, описывающих распространение упругих волн (см. [42, 43]). Математические принципы сеточно-характеристического метода основаны на представлении уравнений движения линейно-упругой среды в матричной форме. Для двумерного случая она имеет вид

(7)
$х{{{\mathbf{q}}}_{t}} + {{{\mathbf{A}}}_{1}}{{{\mathbf{q}}}_{x}} + {{{\mathbf{A}}}_{{\text{2}}}}{{{\mathbf{q}}}_{y}} = 0,$
где q – это вектор, состоящий из пяти компонент:

(8)
${\mathbf{q}} = \left[ {\begin{array}{*{20}{c}} {\mathbf{v}} \\ {\mathbf{T}} \end{array}} \right] = {{\left[ {\begin{array}{*{20}{c}} {{{{v}}_{1}}}&{{{{v}}_{2}}}&{\begin{array}{*{20}{c}} {{{T}_{{11}}}}&{{{T}_{{22}}}}&{{{T}_{{12}}}} \end{array}} \end{array}} \right]}^{{\text{T}}}}.$

Матрицы ${{{\mathbf{A}}}_{k}}$, $k = 1,2$, являются квадратными матрицами размера $5 \times 5$. Произведение матриц ${{{\mathbf{A}}}_{k}}$ и вектора q может быть вычислено следующим образом:

(9)
${{{\mathbf{A}}}_{k}}\left[ {\begin{array}{*{20}{c}} {\mathbf{v}} \\ {\mathbf{T}} \end{array}} \right] = - \left[ {\begin{array}{*{20}{c}} {{{\rho }^{{ - 1}}}\left( {{\mathbf{T}} \cdot {\mathbf{n}}} \right)} \\ {\lambda \left( {{\mathbf{v}} \cdot {\mathbf{n}}} \right){\mathbf{I}} + \mu \left( {{\mathbf{n}} \otimes {\mathbf{v}} + {\mathbf{v}} \otimes {\mathbf{n}}} \right)} \end{array}} \right],$
где ${\mathbf{n}}$ является единичным вектором, направленным вдоль направлений $x$, $y$ для матриц ${{{\mathbf{A}}}_{1}}$, ${{{\mathbf{A}}}_{2}}$ соответственно. Чтобы представить эту формулу явного интегрирования по времени, мы используем спектральное разложение матрицы ${{{\mathbf{A}}}_{k}}$. Например, для матрицы ${{{\mathbf{A}}}_{1}}$
(10)
${{{\mathbf{A}}}_{{\text{1}}}} = {{\left( {{{\Omega }_{{\text{1}}}}} \right)}^{{ - {\text{1}}}}}{{\Lambda }_{{\text{1}}}}{{\Omega }_{{\text{1}}}},$
где ${{\Lambda }_{{\text{1}}}}$ – диагональная матрица $5 \times 5$, образованная собственными значениями матрицы ${{{\mathbf{A}}}_{1}}$; $\Omega _{{\text{1}}}^{{ - {\text{1}}}}$ – матрица $5 \times 5$, образованная соответствующими собственными векторами. Обратим внимание, что матрицы ${{{\mathbf{A}}}_{1}}$ и ${{{\mathbf{A}}}_{2}}$ имеют одинаковый набор собственных значений:
(11)
$\left\{ {{{c}_{p}}, - {{c}_{p}}, - {{c}_{s}},{{c}_{s}},0,0,} \right\}.$
Здесь ${{{\text{c}}}_{p}}$ – продольная скорость волны, равная
(12)
$\sqrt {\left( {{{\rho }^{{ - 1}}}\left( {\lambda + 2\mu } \right)} \right)} ,$
и ${{c}_{s}}$ – поперечная скорость волны, равная

(13)
$\sqrt {\left( {{{\rho }^{{ - 1}}}\mu } \right)} .$

В [30, 44] показано, что решение уравнения (7) для вектора q в координатах $xиy$ для следующего шага по времени может быть записано следующим образом:

(14)
${\mathbf{q}}\left( {t + \tau ,x,y} \right){\text{ = }}\sum\limits_{j{\text{ = 1}}}^J \,{{{\mathbf{X}}}_{{{\text{1,}}j}}}{\mathbf{q}}\left( {t,x - {{\Lambda }_{{{\text{1,}}j}}}\tau ,y} \right),$
(15)
${\mathbf{q}}\left( {t + \tau ,x,y} \right){\text{ = }}\sum\limits_{j{\text{ = 1}}}^J \,{{{\mathbf{X}}}_{{{\text{2,}}j}}}{\mathbf{q}}\left( {t,x,y - {{\Lambda }_{{{\text{2,}}j}}}\tau } \right).$
Здесь $\tau $ – шаг интегрирования по времени, ${{{\mathbf{X}}}_{{1,j}}}$ и ${{{\mathbf{X}}}_{{2,j}}}$ – характеристические матрицы, которые можно выразить с использованием компонентов матриц ${{{\mathbf{A}}}_{1}}$, ${{{\mathbf{A}}}_{2}}$ и их собственных значений следующим образом:
(16)
${{{\mathbf{X}}}_{{i{\text{,}}j}}} = {{\varpi }_{{\text{*}}}}_{{i{\text{,}}j}}{{\varpi }_{{i,j}}},\quad i = 1,2,$
где ${{\varpi }_{{\text{*}}}}_{{i{\text{,}}j}}$ есть $j$-я колонка матрицы ${{\left( {{{\Omega }_{{{\kern 1pt} i}}}} \right)}^{{ - 1}}}$, ${{\varpi }_{{i{\text{,}}j}}}$ есть $j$-я строчка матрицы $\Omega {{{\kern 1pt} }_{i}}$. Скалярные компоненты вектора ${{\varpi }_{{\text{*}}}}_{{i{\text{,}}j}}$ определяются следующим образом:

(17)
${{\omega }_{{1,2}}} = {{\left( {{{\Omega }_{{\text{1}}}}\left[ {\begin{array}{*{20}{c}} {\mathbf{v}} \\ {\mathbf{T}} \end{array}} \right]} \right)}_{{{\text{1,2}}}}} = {\mathbf{n}} \cdot {\mathbf{v}} \mp {{\left( {{{c}_{p}}\rho } \right)}^{{ - {\text{1}}}}}\left( {{{{\mathbf{N}}}_{{00}}} * {\mathbf{T}}} \right),$
(18)
${{\omega }_{{3,4}}} = {{{\mathbf{n}}}_{{\text{1}}}} \cdot {\mathbf{v}} \mp {{\left( {{{c}_{s}}\rho } \right)}^{{ - {\text{1}}}}}\left( {{{{\mathbf{N}}}_{{{\text{01}}}}} * {\mathbf{T}}} \right),$
(19)
${{\omega }_{{5,6}}} = {{{\mathbf{n}}}_{{\text{2}}}} \cdot {\mathbf{v}} \mp {{\left( {{{c}_{s}}\rho } \right)}^{{ - {\text{1}}}}}\left( {{{{\mathbf{N}}}_{{{\text{02}}}}} * {\mathbf{T}}} \right),$
(20)
${{\omega }_{7}} = {{{\mathbf{N}}}_{{{\text{12}}}}} * {\mathbf{T}},$
(21)
${{\omega }_{8}} = \left( {{{{\mathbf{N}}}_{{{\text{11}}}}} - {{{\mathbf{N}}}_{{{\text{22}}}}}} \right) * {\mathbf{T}},$
(22)
${{\omega }_{9}} = \left( {{{{\mathbf{N}}}_{{{\text{11}}}}} + {{{\mathbf{N}}}_{{{\text{22}}}}} - \frac{{{\text{2}}\lambda }}{{\lambda + {\text{2}}\mu }}{{{\mathbf{N}}}_{{00}}}} \right) * {\mathbf{T}}.$

В уравнениях (17)–(22) звездочка $ * $ обозначает свертку двух тензоров ранга $2$. Выражения (16)–(17) могут быть использованы для нахождения решения, вектора q, в любой момент времени, $t + \tau $, при заданных начальных условиях. Таким образом, получаем прямой алгоритм пошагового выполнения – численное моделирование распространения упругих волн в неоднородных средах. Для численного интегрирования в работе применяется сеточно-характеристическая схема третьего порядка точности (см. [30]).

2.3. Метод наложенных сеток

Основная цель работы – разработка нового подхода, основанного на учете различных неоднородностей с использованием наложенных сеток. В качестве неоднородностей рассматриваются единичные поры, представляющие собой замкнутое включение какой-либо формы. При использовании метода наложенных сеток, он же метод химерных сеток, используются основная, как правило, прямолинейная регулярная сетка, и, как минимум, одна наложенная сетка. На фиг. 1а представлен пример реализации таких сеток: роль основной сетки выполняет прямолинейная регулярная сетка, ячейки который являются квадратами с заданной длиной ребер, роль наложенной сетки – криволинейная сетка, образующая кольцо. В процессе моделирования распространения упругой волны, после расчета значений напряженности на следующем временном шаге для всех узлов основной сетки, происходит интерполяция этих значений в соответствующие узлы наложенной сетки и уже после этого рассчитывается распространение возмущения в рамках наложенной сетки. Вычисление следующего временного шага завершается после того, как происходит обратная интерполяция из наложенной сетки в основную. Заметим, что для корректности процедуры интерполяции между вычислительными сетками необходимо, чтобы основная и наложенная сетки имели идентичные характеристики среды. Достаточно подробно метод наложенных сеток исследован в [27], там же приводится исследование его точности при использовании сеточно-характеристического метода.

Фиг. 1.

Расположение наложенных сеток для описания пор: (а) – пример основной регулярной прямолинейной сетки (черный цвет) и наложенной на нее криволинейной (зеленый цвет) для случая моделирования полого отверстия, красным цветом отмечены замкнутые границы, к которым применяются периодические граничные условия; (б) – пример основной регулярной прямолинейной сетки (белый цвет), ребра внешней наложенной сетки (зеленый), внутренней наложенной сетки (оранжевый); серым цветом обозначена область основной сетки, имеющая характеристики “заполнителя”.

2.4. Учет единичных отверстий

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

В случае отверстий, заполненных материалом, отличным от материала среды, используются две наложенные сетки в форме кольца, вложенные друг в друга (см. фиг. 1б). Замыкание сеток самих на себя происходит аналогичным образом с применением периодических граничных условий. Основная сетка в такой постановке не является однородной, узлы, находящиеся в центре описываемого отверстия, имеют параметры материала “заполнителя” поры, отличные от параметров материала среды “каркаса”. При этом узлы основной сетки, образующие собой “лестничную” границу заполнителя и каркаса, должны геометрически находиться под внутренней наложенной сеткой. Эта внутренняя наложенная сетка имеет параметры материала заполнителя и не имеет граничных условия на внутренней и внешней сторонах образуемого ею кольца. Внешняя наложенная сетка, напротив, имеет параметры среды каркаса и также не имеет граничных условий, помимо периодических на замыкаемых границах. Между внутренней стороной внешней наложенной сетки и внешней стороной внутренней наложенной сетки реализован контакт двух различных сред (см. (4)). Таким образом, возмущение, распространяющееся в основной сетке в направлении центра заполненного отверстия после интерполяции, попадает на внешнюю стороной внешней наложенной сетки, распространяется по этой сетке, проходит через контакт различных сред и продолжает распространение во внутренней сетке, после которого в результате интерполяции попадает обратно в основную сетку.

Для демонстрации работы предлагаемого метода нами было проведено моделирование распространения упругой волны в присутствии отверстий обоих типов. На фиг. 2 представлена волновая картина распространения волн в образце, имеющем круглое полое отверстие в различные моменты времени. В данном расчете использовалась основная равномерная квадратная сетка размером $101 \times 101$ узлов с шагом ячеек 2 м и криволинейная наложенная сетка размером 85 × 5 узлов с характерным шагом ячеек 2 м, плотность используемого материала 2500 кг/м3, продольная скорость упругой волны $3000$ м/с, поперечная скорость $1500$ м/с. Таким образом, рассматривается пора радиусом $22$ м. При моделировании использовался шаг интегрирования $dt = 2 \times {{10}^{{ - 4}}}$ с, всего было произведено $400$ шагов. Визуализация результата моделирования распространения плоской волны в присутствии заполненного отверстия приведена на фиг. 3. В данном расчете использовались аналогичные параметры основной сетки, материала “каркаса” и шага интегрирования по времени. Параметры материала поры следующие: плотность 2500 кг/м3, продольная скорость упругой волны $5400$ м/с, поперечная скорость $2700$ м/с. Размеры внешней и внутренней наложенных сеток совпадают и составляют $83 \times 5$ узлов. Характерный размер ячеек внутренней наложенной сетки также составляет $2$ м, а внешней 2.3 м. Плоская волна с частотой $100$ Гц, заданная функцией Рикера, движется параллельно оси $Oy$.

Фиг. 2.

Распространение упругой волны в присутствии полого отверстия.

Фиг. 3.

Распространение упругой волны в присутствии заполненного отверстия.

3. ВЕРИФИКАЦИЯ МЕТОДА

Рассмотрим вопросы верификации предложенного метода. В первую очередь было проведено сравнение решения задачи Лэмба сеточно-характеристическим методом с использованием метода наложенных сеток с аналогичными постановками известных решений. Целью данного шага верификации является демонстрация применимости сеточно-характеристического метода с использованием метода наложенных сеток для корректного моделирования распространения эластичных волн линейно в упругой среде с учетом граничных условий свободной поверхности. За один из эталонных результатов было взято решение, полученное сеточно-характеристическим методом с регулярной прямолинейной сеткой и граничными условиями свободной поверхности (см. [30]). Другим решением, с которым мы проводили сравнение своего метода, было решение, полученное программным комплексом specfem2d (см. [45]), основанное на методе спектральных элементов с использованием неструктурных сеток. Подробная методика и результаты данных сравнений были описаны в [46]. Кроме того, в [29] проведено сравнение результатов моделирования распространения упругой волны под синусоидальной свободной поверхностью с аналогичной постановкой в программном комплексе specfem2d. В результате этих исследований показано, что применение метода наложенных сеток при использовании сеточно-характеристического метода позволяет достигать результаты, незначительно отличающиеся от результатов работы эталонных моделей, приведенных выше.

Для верификации реализации нашего метода, позволяющего учитывать единичные отверстия при моделировании распространения упругой волны в неоднородной среде, была проведена серия расчетов, проверяющих симметрию реализации системы “наложенная сетка–основная сетка”, в случае периодического граничного условия на замыкаемых границах наложенной сетки, которая образует кольцо, описывающее пору. Верификация проводилась для обоих случаев: полого отверстия с граничными условиями свободной поверхности на внутренней поверхности и отверстия, заполненного упругим материалом, отличным от материала окружающей среды.

Описание эксперимента. В постановке задачи верификации единичного отверстия рассматривается распространение упругих волн в квадратной области в присутствии круглого отверстия, параметры моделирования описаны в п. 2.4. В дополнение к описанной выше постановке вокруг наложенной сетки на расстоянии 5 м от внешнего края сетки располагается $36$ приемников. На описанное отверстие по основной сетке набегают 2 плоские волны, частотой $100$ Гц, симметрично относительно прямой, проходящей через центр отверстия под переменным углом, постановка численного эксперимента приведена на фиг. 4. Сравнивалась сумма ошибок принятого приемниками сигнала по норме ${{L}_{2}}$ в течение периода моделирования. Что выражается формулами

(23)
${{L}_{2}}_{{r,x}} = \frac{{\sqrt {\sum\limits_{t = 0}^{{\text{steps}}} \,{{V}_{{x,t,r}}} - V_{{x,t,r + {\text{cnt}}/2}}^{2}} }}{{\sqrt {\sum\limits_{t = 0}^{{\text{steps}}} \,V_{{x,t,r}}^{2}} }},$
(24)
${{L}_{2}}_{x} = \sum\limits_{r = 1}^{{\text{cnt}}/2} \,{{L}_{2}}_{{r,x}},$
(25)
${{L}_{2}}_{{r,y}} = \frac{{\sqrt {\sum\limits_{t = 0}^{{\text{steps}}} \,{{V}_{{y,t,r}}} - V_{{y,t,r + {\text{cnt}}/2}}^{2}} }}{{\sqrt {\sum\limits_{t = 0}^{{\text{steps}}} \,V_{{y,t,r}}^{2}} }},$
(26)
${{L}_{2}}_{y} = \sum\limits_{r = 1}^{{\text{cnt/2}}} \,{{L}_{2}}_{{r,y}},$
где ${{V}_{{y,t,r}}}$ – возмущение, сонаправленное с осью $Ox$, на $t$-м шаге вычислений в точке расположения приемника под номером $r$; steps – общее количество итераций по времени в текущем расчете; cnt – суммарное количество приемников в эксперименте.

Фиг. 4.

Визуализация постановки задачи верификации симметричности отверстия: белым обозначены границы наложенной сетки; желтым и зеленым – попарно сравниваемые группы приемников; черная линия демонстрирует принцип, по которому сравниваются приемники.

Была проведена серия экспериментов при различных углах наклона фронта плоской волны, относительно осей равномерной регулярной основной сетки, для изучения анизотропии использованного метода. Диапазон измерений от 0° до 90°, более широкий диапазон не потребовался ввиду симметричности расположения источников волн и равномерности основной сетки. Результаты проведенных серий экспериментов для полого и заполненного отверстий представлены на фиг. 5а и 5б соответственно. На приведенных графиках представлены экспериментально вычисленные зависимости (24) и (26) от угла между фронтами плоских волн, распространяющихся в направлении отверстия. Полученные графики на фиг. 5а и 5б демонстрируют, что ошибка по норме ${{L}_{2}}$ отдельно взятой компоненты возмущения максимальна в случае перпендикулярности направления распространения упругой волны и компоненты возмущения. Из абсолютных величин ошибки по норме ${{L}_{2}}$ для полого и заполненного отверстий можно сделать вывод, что интерполяция между основной и наложенной сетками не вносит существенной ошибки в данных экспериментах (менее 1%).

Фиг. 5.

Суммарная ошибка по норме ${{L}_{2}}$ при различных углах между осями основной сетки и фронтом волны: (а) – полое отверстие, (б) – заполненное отверстие.

4. ПРИМЕРЫ РАСЧЕТОВ

В данном разделе рассмотрим примеры, демонстрирующие возможности разработанного подхода для численного решения задач распространения волновых возмущений в пористых средах.

4.1. Прохождение фронта волны через пористый слой

Рассмотрим вопрос прохождения фронта плоской волны через пористый слой. На основе проведенного численного эксперимента могут быть определены эффективные упругие параметры пористой среды. Пример показывает возможности разработанного подхода для проведения такого рода исследований. Для моделирования тонкого пористого слоя каждая пора выделяется и рассчитывается отдельно. Рассматривается тонкий пористый слой, состоящий из полых отверстий, представленных кругами одинакового радиуса. Слой находится между двумя полупространствами. Упругие параметры полупространств и каркаса пористого слоя совпадают. На фиг. 6а представлена схема расположения основных сеток. “Верхняя” и “нижняя” сетки (выделены черным) – равномерные регулярные сетки, не имеющие интерполяционного взаимодействия с другими сетками. Они служат для описания среды, в которой находятся источник и приемник соответственно. Между ними располагаются сетки-“сегменты”, смещенные на произвольное количество узлов друг относительно друга. Каждая сетка-“сегмент” необходима для описания отдельного полого отверстия в среде, центр которого совпадает с геометрическим центром “сегмента”. Метод описания поры с использованием наложенной сетки описан выше. Смещение “сегментов” друг относительно друга необходимо для исследования различных конфигураций расположения отверстий и исключения эффектов, вызванных их структурой, таких как интерференция. Для наглядной визуализации проведенного эксперимента на фиг. 7 приведена картина распространения линейно-упругой волны в тонком пористом слое в различные моменты времени. При моделировании описанной выше задачи использовалось $50$ сегментов сеток, каждый из которых содержал наложенную сетку, описывающую единичную пору радиуса от 11 до 41 м. Варьируя радиус пор между экспериментами (в каждом отдельном эксперименте радиусы всех пор одинаковы), эксперимент позволяет измерять свойства такого пористого слоя в зависимости от пористости среды в диапазоне от 4 до $52\% $. В предлагаемом методе при увеличении диапазона пористости квадратично растет вычислительная сложность моделирования. Каждый из этих сегментов имеет размер $51 \times 51$ узлов с длиной ребра $2$ м. Наложенные сетки также имеют характерные размеры ячеек $2$ м и имеют размеры от $44 \times 6$ до $162 \times 6$ узлов в зависимости от радиуса пор в эксперименте. “Верхняя сетка” и “нижняя сетка” , представленные на фиг. 6а, имеют размеры $51 \times 501$ и $51 \times 251$ соответственно. Размер ячеек в этих сетках также составляет $2$ м. В данном эксперименте используется источник Рикера с длиной волны $6000$ м. Приемник, регистрирующий напряженность, расположен на высоте $200$ м от нижнего края “нижней сетки”. Физические параметры среды во всех сетках идентичны: плотность материала $2500$ кг/м3, продольная скорость упругой волны $3000$ м/с, поперечная скорость $1500$ м/с.

Фиг. 6.

Распространение упругой волны в пористом слое. Белыми точками обозначены поры. Градиентом показана амплитуда упругой волны.

Фиг. 7.

Расположение вычислительных сеток при моделировании тонкого пористого слоя: (а) – взаимное расположение основных сеток; (б) – периодически граничные условия (крайние узлы отмечены красным) и контакт между сетками-“сегментами” (крайние узлы отмечены зеленым).

Существует ряд аналитических теорий, позволяющих оценить значение модуля упругости материла в зависимости от его степени пористости $p = {{S}_{{{\text{pore}}}}}{\text{/}}({{S}_{0}} + {{S}_{{{\text{pore}}}}})$, где ${{S}_{{{\text{pore}}}}}$ – площадь пор, ${{S}_{0}}$ – площадь, занимаемая средой каркаса. Наиболее известные из них модели – Reuss (27), Voigt (28) и Hashin–Shtrikman (29) (см. [47]), которые в нашем случае среды, состоящей из двух материалов, принимают вид

(27)
$\frac{1}{K} = \frac{p}{{{{K}_{{{\text{pore}}}}}}} + \frac{{1 - p}}{{{{K}_{0}}}},$
(28)
$K = p{{K}_{{{\text{pore}}}}} + (1 - p){{K}_{0}},$
(29)
$K = {{K}_{0}} + \frac{p}{{{{{({{K}_{{{\text{pore}}}}} - {{K}_{0}})}}^{{ - 1}}} + (1 - p)({{K}_{0}} + 4{{\mu }_{0}}{\text{/}}{{{3)}}^{{ - 1}}}}},$
где ${{K}_{0}},\;{{K}_{{{\text{pore}}}}},\;K$ – модули упругости материала каркаса, заполнителя и получившегося образца соответственно; ${{\mu }_{0}}$ – модуль сдвига материала каркаса. В приведенном эксперименте приемник расположен в нижней части моделируемой области, что позволило измерить зависимость величины возмущения, вызванного распространением от времени. Данная зависимость представлена на фиг. 8а. Используя зависимости, представленные на этом графике, учитывая расположение приемника, величины участков среды, в которых отсутствуют вкрапления пор, и параметры используемого источника Рикера, вычислены значения продольной скорости упругой волны для каждого значения пористости в эксперименте. Поскольку в рассматриваемой постановке измерение поперечной скорости упругой волны не представляется возможным, для ее оценки будем использовать приближенное значение модуля сдвига Hashin–Shtrikman (30) (см. [47])
(30)
$\mu = {{\mu }_{0}} + \frac{p}{{{{{({{\mu }_{{{\text{pore}}}}} - {{\mu }_{0}})}}^{{ - 1}}} + \frac{{2p({{K}_{0}} + 2{{\mu }_{0}})}}{{5{{\mu }_{0}}({{K}_{0}} + 4{{\mu }_{0}}{\text{/}}3)}}}},$
где ${{\mu }_{{{\text{pore}}}}},\;\mu $ – модули сдвига материала заполнителя и получившегося образца соответственно. Получившаяся зависимость модуля сдвига от пористости материала представлена на фиг. 8б и отмечена красными крестами. Для сравнения на фиг. 8б представлены теоретические оценки зависимости модуля упругости, согласно моделям Reuss (27), Voigt (28) и Hashin–Shtrikman (29). Как следует из графика, предлагаемый метод дает результаты, укладывающиеся в теоретические оценки. Заметим, что вычислительная сложность данного метода растет квадратично с приближением пористости к предельным значениям $0$ и $1$, чем обусловливается диапазон измерений на фиг. 8.

Фиг. 8.

Результаты моделирования тонкого пористого слоя при различных значениях пористости: (а) – зависимость значения возмущения на приемнике от времени для различных значений пористости; (б) – графики сравнения зависимости объемного модуля упругости от пористости среды, полученные в результате моделирования предлагаемым методом с теоретическими границами Reuss, Voigt и моделью Hashin–Shtrikman.

4.2. Моделирование сегмента реалистичной породы

Предлагаемый в статье метод позволяет при моделировании распространения упругих волн использовать различные геометрические формы отверстий. Для демонстрации потенциальной области применения данного метода рассмотрим пример, который воспроизводит постановку, используемую в ([48], Fig. 1), где все поры представлялись эллипсами. Вследствие высокой трудоемкости получения точных координат пор, для моделирования рассмотрим только верхнюю половину представленного образца. Следует отметить, что данное ограничение связано с трудностью оцифровки вручную постановки эксперимента из работы [48] и не связано с ограничениями метода. В данном эксперименте используются одна основная регулярная прямолинейная сетка и 120 наложенных криволинейных сеток различных размеров. Количество пор в эксперименте равно соответственно 120. Регулярная сетка имеет размер $850 \times 1300$ узлов, размером 1 × 1 м. Шаг интегрирования по времени составил $1.5 \times {{10}^{{ - 5}}}$ с. Наложенные сетки, описывающие поры, имеют различные размеры в диапазоне от $39 \times 5$ до $304 \times 5$ узлов, что, в свою очередь, позволяет моделировать полые отверстия в форме эллипса с размерами полуосей от $7.1,5.0$ м до $59.1,36.1$ м. Физические параметры среды во всех сетках идентичны и составляют: плотность материала $2500$ кг/м3, продольная скорость упругой волны $3000$ м/с, поперечная скорость $1500$ м/с. Использовался источник с длиной волны $300$ м. Волновые картины данного расчета в различные моменты времени представлены на фиг. 9.

Фиг. 9.

Результаты моделирования реалистичной пористой среды, волновая картина в различные моменты времени.

Данный эксперимент демонстрирует возможности разработанного алгоритма. Предложенный метод позволяет задавать сложные конфигурации пористых структур в достаточно большом количестве и проводить расчет за разумное время.

5. ЗАКЛЮЧЕНИЕ

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

В качестве примеров использования разработанного алгоритма приведена серия тестовых расчетов. Один из них – способ, позволяющий производить моделирование пористого слоя с различными параметрами пористости и толщины. Для частного случая такой конфигурации приводится эксперимент по численному исследованию упругих параметров пористого слоя в зависимости от пористости среды. В заключение работы продемонстрирована возможность моделирования реалистичной пористой среды, состоящей из пор различных форм и размеров. В результате можно сделать заключение, что рассмотренный метод наложенных сеток имеет потенциал для изучения свойств пористых сред различного состава и конфигураций.

В качестве дальнейшего развития данного направления стоит отметить перенос алгоритмов на трехмерный случай и распараллеливания в системах с общей и распределенной памятью.

Список литературы

  1. Qi Yingkai, Chen Xuehua, Zhao Qingwei, Luo Xin, Feng Chunqiang. Seismic wave modeling of fluid-saturated fractured porous rock: Including fluid pressure diffusion effects of discrete distributed large-scale fractures // EGUsphere. 2023. № 1. P. 1–26.

  2. Liu Jiong, Wei Xiu Cheng, Ji Yu Xin, Chen Tian Sheng, Liu Chun Yuan, Zhang Chun Tao, Dai Ming Gang. An analysis of seismic scattering attenuation in a random elastic medium // Appl. Geophys. 2011. V. 8. № 12. P. 344–354.

  3. Wei Yijun, Ba Jing, Carcione J.M. Stress effects on wave velocities of rocks: Contribution of crack closure, squirt flow and acoustoelasticity // J. Geophys. Res.: Solid Earth. 2022. V. 127. № 10. P. 2022JB025253.

  4. Gassmann F. On elasticity of porous media // Classics of Elastic Wave Theory. 2007. № 1. P. 389–408.

  5. Berryman J.G. Origin of gassmann’s equations // Geophysics. 1999. V. 64. P. 1627–1629.

  6. Biot M.A. Theory of propagation of elastic waves in a fluid-saturated porous solid. ii. higher frequency range // J. Acoustical Soc. Am. 1956. V. 28. № 6. P. 179.

  7. Dvorkin J., Nur A. Dynamic poroelasticity: a unified model with the squirt and the biot mechanisms // Geophysics. 1993. V. 58. P. 524–533.

  8. Dvorkin J., Nolen-Hoeksema R., Nur A. The squirt-flow mechanism: macroscopic description // Geophysics. 1994. V. 59. P. 428–438.

  9. Dvorkin J., Mavko G., Nur A. Squirt flow in fully saturated rocks // Geophysics. 1995. V. 60. P. 97–107.

  10. Yang Dinghui, Zhang Zhongjie. Effects of the biot and the squirt-flow coupling interaction on anisotropic elastic waves // Chinese Sci. Bull. 2000. V. 45. P. 2130–2138.

  11. Pride S.R., Berryman J.G., Harris J.M. Seismic attenuation due to wave-induced flow // J. Geophys. Res.: Solid Earth. 2004. № 1. P. 109.

  12. Müller T.M., Toms-Stewart J., Wenzlau F. Velocity-saturation relation for partially saturated rocks with fractal pore fluid distribution // Geophys. Res. Lett. 2008. V. 35. № 5. P. 9306.

  13. Huang Xingguo, Greenhalgh Stewart, Han Li, Liu Xu. Generalized effective biot theory and seismic wave propagation in anisotropic, poroviscoelastic media // J. Geophys. Res.: Solid Earth. 2022. V. 127. № 3. P. 2021JB023590.

  14. Jing B.A., Carcione J.M., Hong Cao, Qi-Zhen Du, Zhen-Yu Yuan, Ming-Hui Lu. Velocity dispersion and attenuation of p waves in partially-saturated rocks: Wave propagation equations in double-porosity medium // Chinese J. Geophys. 2012. V. 55. № 1. P. 219–231.

  15. Amalokwu K., Best A.I., Sothcott J., Chapman M., Minshull T., Li X.Y. Water saturation effects on elastic wave attenuation in porous rocks with aligned fractures // Geophys. J. Inter. 2014.V. 197. № 5. P. 943–947.

  16. Sun Weitao, Ba Jing, Müller T.M., Carcione J.M., Cao Hong. Comparison of p-wave attenuation models of wave-induced flow // Geophys. Prospect. 2015. V. 63. № 3. P. 378–390.

  17. Kachanov M. Elastic solids with many cracks and related problems // Adv. Appl. Mech. 1993. V. 30. P. 259–445.

  18. Gu’eguen Y., Sarout J. Crack-induced anisotropy in crustal rocks: Predicted dry and fluid-saturated thomsen’s parameters // Physics of the Earth and Planetary Interiors. 2009. V. 172. P. 116–124.

  19. Gu’eguen Y., Sarout J. Characteristics of anisotropy and dispersion in cracked medium // Tectonophysics. 2011. V. 503. № 4. P. 165–172.

  20. Dorovsky N.V. Continual theory of filtration // Sov. Geology and Geophysics. 1989. P. 34–39.

  21. Blokhin A.M., Dorovskii V.N. Mathematical modelling in the theory of multivelocity continuum. 1995. P. 183.

  22. Dorovsky V.N., Perepechko Yu.V., Fedorov A.I. Stoneley waves in the biot–johnson and continuum filtration theories // Russian Geology and Geophys. 2012. V. 53. № 5. P. 475–483.

  23. Guo Zhiqi, Qin Xiaoying, Zhang Yiming, Niu Cong, Wang Di, Ling Yun. Numerical investigation of the effect of heterogeneous pore structures on elastic properties of tight gas sandstones // Frontiers in Earth Science. 2021. V. 9. № 4. P. 219.

  24. Li Tianyang, Wang Zizhen, Yu Nian, Wang Ruihe, Wang Yuzhong. Numerical study of pore structure effects on acoustic logging data in the borehole environment. 2020. V. 28. № 5. https://doi.org/10.1142/S0218348X20500498

  25. Ozotta O., Saberi M.R., Kolawole O., Malki M.L., Rasouli V., Pu Hui. Pore morphology effect on elastic and fluid flow properties in bakken formation using rock physics modeling // Geomechanics and Geophysics for Geo-Energy and Geo-Resources. 2022. V. 8. № 12. P. 1–19.

  26. Aney Sh., Rege A.The effect of pore sizes on the elastic behaviour of open-porous cellular materials // Math. and Mech. of Solids. 2022. № 10.

  27. Khokhlov N., Favorskaya A., Stetsyuk V., Mitskovets I. Grid-characteristic method using Chimera meshes for simulation of elastic waves scattering on geological fractured zones // J. Comput. Phys. 2021. V. 446. P. 110637.

  28. Khokhlov N.I., Favorskaya A., Furgailo V. Grid-characteristic method on overlapping curvilinear meshes for modeling elastic waves scattering on geological fractures // Minerals. 2022. V. 12. № 12. P. 1597.

  29. Mitskovets I., Stetsyuk V., Khokhlov N. Novel approach for modeling curved topography using overset grids and grid-characteristic method // European Association of Geoscientists Engineers 2020. № 12. P. 1–5.

  30. Favorskaya A.V., Zhdanov M.S., Khokhlov N.I., Petrov I.B. Modelling the wave phenomena in acoustic and elastic media with sharp variations of physical properties using the grid-characteristic method // Geophys. Prospect. 2018. V. 66. № 10. P. 1485–1502.

  31. Magomedov K.M., Kholodov A.S. The construction of difference schemes for hyperbolic equations based on characteristic relations // USSR Comput. Math. and Math. Phys. 1969. V. 9. № 2. P. 158–176.

  32. Korotin P.N., Petrov I.B., Pirogov V.B., Kholodov A.S. On a numerical solution of related problems of supersonic flow over deformable shells of finite thickness // USSR Comput. Math. and Math. Phys. 1987. V. 27. № 4. P. 181–188.

  33. Petrov I.E., Kholodov A.S. Numerical study of some dynamic problems of the mechanics of a deformable rigid body by the mesh-characteristic method // USSR Comput. Math. and Math. Phys. 1984. V. 24. № 3. P. 61–73.

  34. Petrov I.B., Tormasov A.G., Kholodov A.S. On the use of hybrid grid-characteristic schemes for the numerical solution of three-dimensional problems in the dynamics of a deformable solid // USSR Comput. Math. and Math. Phys. 1990. V. 30. № 4. P. 191–196.

  35. Kvasov I.E., Pankratov S.A., Petrov I.B. Numerical simulation of seismic responses in multilayer geologic media by the grid-characteristic method // Math. Model. Comput. Simulat. 2011. V. 3. № 2. P. 196–204.

  36. Muratov M.V., Petrov I.B. Estimation of wave responses from subvertical macrofracture systems using a grid characteristic method // Math. Model. Comput. Simulat. 2013. V. 5. № 5. P. 479–491.

  37. Petrov I.B., Khokhlov N.I. Modeling 3D seismic problems using high-performance computing systems // Math. Model. Comput. Simulat. 2014. V. 6. № 4. P. 342–350.

  38. Aki Keiiti, Richards P.G. Quantitative seismology, 2nd ed. // Quse. 2022. V. 68. P. 1546–1546.

  39. LeVeque R.J. Finite volume methods for hyperbolic problems // Finite Volume Methods for Hyperbolic Problems. 2002. V. 8.

  40. Zhdanov M.S. Geophysical inverse theory and regularization problems. 2002. P. 609.

  41. Zhdanov M.S. Inverse theory and applications in geophysics // Inverse Theory and Appl. Geophys. 2015. V. 9. P. 1–704.

  42. Petrov I.B., Favorskaya A.V., Sannikov A.V., Kvasov I.E. Grid-characteristic method using high-order interpolation on tetrahedral hierarchical meshes with a multiple time step // Math. Model. Comput. Simulat. 2013. V. 5. № 9. P. 409–415.

  43. Golubev V.I., Petrov I.B., Khokhlov N.I. Numerical simulation of seismic activity by the grid-characteristic method // Comput. Math. and Math. Phys. 2013. V. 53. № 10. P. 1523–1533.

  44. Khokhlov N.I., Golubev V.I. On the class of compact grid-characteristic schemes // Smart Innovation, Systems and Technolog. 2019. V. 133. P. 64–77.

  45. Komatitsch D., Tromp J. Introduction to the spectral element method for three-dimensional seismic wave propagation // Geophys. J. Inter. 1999. V. 139. № 12. P. 806–822.

  46. Khokhlov N.I., Stetsyuk V.O., Mitskovets I.A. Overset grids approach for topography modeling in elastic-wave modeling using the grid-characteristic method // Компьют. иссле д. и моделирование. 2019. V. 11. P. 1049–1059.

  47. Mavko G., Mukerji T., Dvorkin J. Effective elastic media: bounds and mixing laws // The Rock Physics Handbook. 2009. № 3. P. 169–228.

  48. Wang Z., Wang R., Li T., Qiu Hao, Wang F. Pore-scale modeling of pore structure effects on p-wave scattering attenuation in dry rocks // PLoS ONE. 2015. № 5. P. 10.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал