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

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

Граничные и контактные условия повышенного порядка аппроксимации для сеточно-характеристических схем в задачах акустики

А. В. Шевченко 1, 2, В. И. Голубев 1, 2, *, **

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

2 ИАП РАН
123056 Москва, ул. 2-я Брестcкая, 19/18, Россия

* E-mail: w.golubev@mail.ru
** E-mail: golubev.vi@mipt.ru

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

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

Аннотация

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

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

1. ВВЕДЕНИЕ

Математическое моделирование является важным практическим инструментом в ряде научных и прикладных областей: изготовление новых материалов, авиапромышленность, оптимизация дорожного трафика и др. Одним из возможных подходов является формулирование определяющей системы уравнений в полных или частных производных, возможно высшего порядка, описывающей поведение рассматриваемой системы. Для ее численного решения в дальнейшем могут быть применены различные численные методы: конечно-разностные схемы [1], конечно-объемные или спектральные методы [2], сеточно-характеристические схемы [3]. Практически важными свойствами расчетных схем являются аппроксимация и устойчивость. Усилия различных исследователей направлены на построение схем повышенного порядка аппроксимации. Например, хорошего результата удалось добиться использованием схем MUSCL [4] для линейных задач [5, 6]. В работе [7] продемонстрирована возможность их адаптации на нелинейные задачи с использованием корректной реконструкции потоков [8]. В работе [9] представлены граничные условия третьего порядка точности для нелинейных законов сохранения в одномерной постановке. Для обеспечения монотонности схемы возможно применение лимитеров [10] или процедуры гибридизации решения [11].

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

В настоящей работе рассматривается вопрос задания граничных и контактных условий в рамках сеточно-характеристического подхода, обеспечивающих сохранение повышенного порядка аппроксимации вплоть до границы включительно. Данное исследование является продолжением работ авторов по сеточно-характеристическим методам. В частности, вопрос сохранения порядка при построении многомерных схем рассмотрен в работах [14, 15], а для задач с кусочно-постоянными коэффициентами – в работе [16]. Используемое в работе акустическое приближение часто применяется на практике в задачах сейсмической разведки ввиду меньшей вычислительной сложности задачи. Таким образом, рассматриваемая в работе проблема обладает прикладной значимостью.

2. ОСНОВЫ СЕТОЧНО-ХАРАКТЕРИСТИЧЕСКОГО МЕТОДА

В данном разделе представлено достаточно краткое описание основ сеточно-характеристического метода. Для более подробного описания можно порекомендовать следующие работы научной группы Петрова Игоря Борисовича (МФТИ): [17–20].

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

(1)
${{{\mathbf{q}}}_{t}} + {{A}_{1}}{{{\mathbf{q}}}_{x}} + {{A}_{2}}{{{\mathbf{q}}}_{y}} + {{A}_{3}}{{{\mathbf{q}}}_{z}} = {\mathbf{f}}(x,y,z,t),$
где компонентами вектора ${\mathbf{q}} = {\mathbf{q}}(x,y,z,t)$ являются все неизвестные функции координат и времени, ${\mathbf{f}}$ задает известную правую часть. Матрицы ${{A}_{1}}$, ${{A}_{2}}$, ${{A}_{3}}$ определяются параметрами среды и допускают переход к базису из собственных векторов: ${{A}_{i}} = \Omega _{i}^{{ - 1}}{{\Lambda }_{i}}{{\Omega }_{i}}$, что следует из гиперболичности системы уравнений. Далее будем полагать, что компоненты матриц ${{A}_{i}}$ не зависят от координат и времени.

Для решения многомерной системы уравнений применяется подход координатного расщепления (расщепление по направлениям, операторное расщепление, operator splitting). На каждом шаге по времени последовательно решаются несколько одномерных систем уравнений вида ${{{\mathbf{q}}}_{t}} + {{A}_{i}}{{{\mathbf{q}}}_{{{{x}_{i}}}}} = {\mathbf{0}}$ и ${{{\mathbf{q}}}_{t}} = {\mathbf{f}}$. Простейшим алгоритмом для схемы расщепления является следующий

Алгоритм

Шаг 1. Найти ${{{\mathbf{q}}}^{{n + 1/4}}}$, сделав шаг по времени $dt$ для одномерной системы уравнений ${{{\mathbf{q}}}_{t}} + {{A}_{1}}{{{\mathbf{q}}}_{x}} = {\mathbf{0}}$.

Шаг 2. Найти ${{{\mathbf{q}}}^{{n + 2/4}}}$, сделав шаг по времени $dt$ для одномерной системы уравнений ${{{\mathbf{q}}}_{t}} + {{A}_{2}}{{{\mathbf{q}}}_{y}} = {\mathbf{0}}$ с начальными значениями ${{{\mathbf{q}}}^{{n + 1/4}}}$ после шага 1.

Шаг 3. Найти ${{{\mathbf{q}}}^{{n + 3/4}}}$, сделав шаг по времени $dt$ для одномерной системы уравнений ${{{\mathbf{q}}}_{t}} + {{A}_{3}}{{{\mathbf{q}}}_{z}} = {\mathbf{0}}$ с начальными значениями ${{{\mathbf{q}}}^{{n + 2/4}}}$ после шага 2.

Шаг 4. Найти ${{{\mathbf{q}}}^{{n + 1}}}$, сделав шаг по времени $dt$ для системы уравнений ${{{\mathbf{q}}}_{t}} = {\mathbf{f}}(x,y,z,t)$ с начальными значениями ${{{\mathbf{q}}}^{{n + 3/4}}}$ после шага 3.

Отметим, что порядок аппроксимации итоговой многомерной схемы получится не выше первого. Более подробный анализ различных схем расщепления повышенного порядка представлен в работах [14, 15].

Таким образом, базовым элементом сеточно-характеристического метода является численное решение одномерной системы уравнений вида

(2)
${{{\mathbf{q}}}_{t}} + A{{{\mathbf{q}}}_{x}} = {\mathbf{0}}.$
Дифференциальная задача рассматривается в ограниченной области $x \in [a,b]$, $t \in [0,T]$. Следовательно, для полной постановки задачи должны быть заданы начальные условия ${{\left. {\mathbf{q}} \right|}_{{t = 0}}}$ и достаточное число граничных условий.

Для решения задачи во внутренних точках используется разложение матрицы $A = {{\Omega }^{{ - 1}}}\Lambda \Omega $. Столбцы матрицы ${{\Omega }^{{ - 1}}}$ являются собственными векторами матрицы $A$, строки матрицы $\Omega $ являются левыми собственными векторами (“собственными строками”) $A$, $\Omega {{\Omega }^{{ - 1}}} = I$ – единичная матрица. Матрица $\Lambda $ является диагональной, причем на диагонали стоят собственные значения $A$. Подстановка $A = {{\Omega }^{{ - 1}}}\Lambda \Omega $ и введение новых переменных – вектора инвариантов Римана ${\mathbf{w}} = \Omega {\mathbf{q}}$ – приводят (в предположении постоянства $A$) к системе уравнений ${{{\mathbf{w}}}_{t}} + \Lambda {{{\mathbf{w}}}_{x}} = {\mathbf{0}}$. Она является набором независимых уравнений переноса, так как матрица $\Lambda $ диагональная. Для каждого уравнения переноса выполняется характеристическое свойство ${{w}_{i}}(x,t + \tau ) = {{w}_{i}}(x - {{\lambda }_{i}} \cdot \tau ,t)$. Это позволяет находить значения ${{{\mathbf{w}}}^{{n + 1}}}$ с помощью полиномиальной интерполяции на текущем слое по времени ${{\left. {{{{\mathbf{q}}}^{n}}} \right|}_{{t = {{t}^{n}}}}}$. Таким образом, на каждом шаге по времени ${{t}^{n}}$ для каждого узла расчетной сетки ${{x}_{m}}$ выполняем следующее:

1) вычисляем ${{{\mathbf{w}}}^{n}} = \Omega {{{\mathbf{q}}}^{n}}$;

2) вычисляем ${{{\mathbf{w}}}^{{n + 1}}}$ с помощью интерполяции ${{w}_{i}}({{x}_{m}},{{t}^{{n + 1}}}) = {{w}_{i}}(x - {{\lambda }_{i}} \cdot \tau ,t) \approx P(x - {{\lambda }_{i}} \cdot \tau )$, где $P(x)$ – интерполяционный полином;

3) возвращаемся к исходным неизвестным ${{{\mathbf{q}}}^{{n + 1}}} = {{\Omega }^{{ - 1}}}{{{\mathbf{w}}}^{{n + 1}}}$.

Построенный для узла сетки ${{x}_{m}}$ интерполяционный полином $P(x)$ в общем виде определяется на некотором шаблоне $\{ {{x}_{k}}\} $, $k \in \overline {m - M,m + M} $, $M \in \mathbb{N}$. При этом использование многочленов более высокой степени для полиномиальной интерполяции, построенных на более широких шаблонах (с увеличенным $M$), повышает порядок аппроксимации итоговой расчетной схемы. Условие устойчивости соответствует выполнению условия Куранта: $dt < \frac{h}{{{{{\max }}_{i}}{{\lambda }_{i}}}}$.

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

Однако в граничных узлах ${{x}_{0}}$, ${{x}_{N}}$ и в приграничных узлах ${{x}_{1}}$, ${{x}_{{N - 1}}}$, ${{x}_{2}}$, ${{x}_{{N - 2}}}, \ldots $ описанная выше процедура не определена. Далее будут рассмотрены вопросы реализации схемы в этих точках по заданным граничным условиям.

3. СИСТЕМА УРАВНЕНИЙ АКУСТИКИ

В данной работе рассматривалась задача распространения волн в акустической среде, которую можно описать системой вида

(3)
$\begin{gathered} \rho {{{\mathbf{v}}}_{t}} = - \nabla p + {\mathbf{F}}, \\ {{p}_{t}} = - \rho {{c}^{2}}\nabla \cdot {\mathbf{v}}. \\ \end{gathered} $
Здесь неизвестными функциями координат и времени являются вектор скорости частиц среды ${\mathbf{v}}$ и скалярное давление $p$. Вектор ${\mathbf{F}}$ определяет внешнее воздействие на среду, а плотность $\rho $ и скорость распространения акустических волн $c$ являются параметрами среды. Далее дополнительно используется обозначение $Z = \rho c$ для акустического импеданса.

Система уравнений акустики в одномерном однородном случае принимает вид

(4)
$\begin{gathered} \rho {{{v}}_{t}} = - {{p}_{x}}, \\ {{p}_{t}} = - \rho {{c}^{2}}{{{v}}_{x}}. \\ \end{gathered} $
Переход к инвариантам Римана ${\mathbf{w}} = ({{w}^{ - }}\;{{w}^{ + }}{{)}^{{\text{т}}}}$ определяется соотношениями
(5)
${{w}^{ - }} = ({v} - p{\text{/}}Z){\text{/}}2,\quad {{w}^{ + }} = ({v} + p{\text{/}}Z){\text{/}}2.$
Обратный переход к исходным неизвестным осуществляется по формулам

(6)
${v} = {{w}^{ + }} + {{w}^{ - }},\quad p = Z({{w}^{ + }} - {{w}^{ - }}).$

4. РЕАЛИЗАЦИЯ ГРАНИЧНЫХ УСЛОВИЙ

Как было отмечено выше, задача (2) рассматривается в ограниченной области $(x,t) \in [a,b] \times [0,T]$. Следовательно, для корректности задачи должны быть заданы начальные и граничные условия. Чтобы дифференциальная задача была корректной по Адамару, для системы уравнений акустики граничные условия должны быть заданы на обеих границах $x = a$ и $x = b$.

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

Также могут быть поставлены линейные граничные условия с явно заданной правой частью, записываемые в виде ${{\left. {{{G}_{a}}{\mathbf{q}}} \right|}_{{x = a}}} = {{{\mathbf{g}}}_{a}}(t)$, ${{\left. {{{G}_{b}}{\mathbf{q}}} \right|}_{{x = b}}} = {{{\mathbf{g}}}_{b}}(t)$, где ${{G}_{a}}$, ${{G}_{b}}$ – прямоугольные матрицы. В этом случае можно показать, что число условий, равное числу компонент вектора ${\mathbf{g}}$, должно быть равно числу выходящих из области характеристик. Таким образом, для уравнений акустики необходимо задать одно условие на каждой границе, для уравнений упругости – два условия. Например, широко используется условие “свободной границы”. Для акустической среды оно имеет вид ${{\left. p \right|}_{{x = X}}} = 0$ [21], для упругой среды вид $\sum\nolimits_j {{\left. {({{\sigma }_{{ij}}}{{n}_{j}})} \right|}_{{x = X}}} = {\mathbf{0}}$ [22]. Здесь $X$ означает левую границу $x = a$ или правую границу $x = b$, ${\mathbf{n}}$ – вектор внешней нормали, $\sigma $ – тензор напряжений.

Существует множество подходов к заданию в расчетном алгоритме граничных условий. Например, можно использовать следующие:

1) предзаполнение ghost-узлов;

2) модификацию численной схемы в узлах вблизи границы;

3) специальные корректировки после шага схемы во внутренних точках.

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

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

Третий подход предполагает, что после основного расчетного цикла (возможно, с подсчетом предварительного значения на границе) выполняется специальная корректировка, после которой заданное граничное условие будет выполнено.

Отметим, что эти подходы не являются взаимоисключающими и могут использоваться совместно.

Для задания поглощающих граничных условий могут использоваться подходы выше или изменение уравнений/коэффициентов в приграничной области. Например, широко распространено добавление области с затуханием или perfectly matched layer (PML) [23]. В данной работе рассматриваются только способы без введения приграничной области.

4.1. Понижение порядка аппроксимации

Задание неотражающих граничных условий в сеточно-характеристическом методе возможно при помощи обнуления значений, переносимых на характеристиках, которые выходят из расчетной области. Это может быть реализовано в рамках подхода № 1 со следующим правилом заполнения ghost-узлов. Значения в ghost-узлах заполняются значением функции в граничном узле на текущем шаге по времени ${\mathbf{q}}_{m}^{n}$. Такое заполнение приближает интерполируемую величину $w{\kern 1pt} *$ ее значением в узле ${{w}_{m}}$, благодаря чему “добавка” $w{\kern 1pt} *\; - {{w}_{m}} \approx 0$. После этого численная схема применяется ко всем узлам от 0 до $N$, а дополнительные корректировки не используются. Недостаток данного подхода кроется в том, что при использовании широкого шаблона для интерполяции в формулу расчета $w{\kern 1pt} *$ входят как ghost-узлы, так и значения внутри области, вследствие чего $w{\kern 1pt} *\; - {{w}_{m}} \ne 0$.

Задание линейных граничных условий возможно с помощью следующей процедуры:

1) заполнить ghost-узлы для неотражающей границы;

2) выполнить шаг схемы по времени во всех точках;

3) применить специальные корректировки в граничных узлах.

Формулы для корректировок в граничных узлах выводятся в предположении, что перед их применением в узлах хранятся значения, посчитанные только по внутренним характеристикам, обозначаемые ${{{\mathbf{q}}}^{{{\text{in}}}}}$. В этом случае формула корректировки может быть записана следующим образом [24]:

(7)
${{{\mathbf{q}}}^{{n + 1}}} = {{{\mathbf{q}}}^{{{\text{in}}}}} + {{\Omega }^{{{\text{out}}}}}{{[G \cdot {{\Omega }^{{{\text{out}}}}}]}^{{ - 1}}} \cdot ({\mathbf{g}} - G{{{\mathbf{q}}}^{{{\text{in}}}}}).$
В этой формуле прямоугольная матрица ${{\Omega }^{{{\text{out}}}}}$ состоит из тех столбцов матрицы ${{\Omega }^{{ - 1}}}$, которые соответствуют выходящим из расчетной области характеристикам.

Существование обратной матрицы ${{[G \cdot {{\Omega }^{{{\text{out}}}}}]}^{{ - 1}}}$ необходимо и достаточно для корректности заданных граничных условий. Матрицы в формулах корректировки могут быть предвычислены аналитически. Для их применения возможно использование безматричных операций, что обусловливает высокую эффективность этой процедуры.

Данная формула корректировки не вносит в алгоритм дополнительной численной ошибки. Следовательно, точность метода определяется точностью схемы и точностью задания вектора ${{{\mathbf{q}}}^{{{\text{in}}}}}$. При этом одномерная схема может иметь высокий порядок во внутренних узлах, но при интерполяции значений в ghost-узлах вносится существенная ошибка. Результаты тестовых расчетов, подтверждающие данные утверждения, приведены в соответствующей главе.

Для повышения общего (включая граничные точки) порядка схемы можно выбрать один из двух путей:

1) вывести формулы заполнения ghost-узлов, не требующие дополнительных последующих корректировок;

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

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

4.2. Сохранение порядка аппроксимации

Рассмотрим процедуру задания линейных граничных условий на примере уравнения акустики и границы с заданным давлением $P(t)$. При этом “свободная граница” ($p \equiv 0$) является частным случаем рассматриваемых граничных условий. Возможно провести аналогичные построения для задания нормальной к границе компоненты скорости $V(t)$.

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

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

Без ограничения общности рассмотрим правую границу. Точка ${{x}_{m}}$ лежит на границе, точки ${{x}_{{m - 1}}}$, ${{x}_{{m - 2}}}$, … – внутри расчетной области, точки ${{x}_{{m + 1}}}$, ${{x}_{{m + 2}}}$, … – вне расчетной области, т.е. задают ghost-узлы.

На самой границе согласно (6) ${{p}_{m}} = Z({{w}^{ + }} - {{w}^{ - }}) = P(t) \Rightarrow $

(8)
$\begin{gathered} w_{m}^{ + }(t) = w_{m}^{ - }(t) + P(t){\text{/}}Z, \\ w_{m}^{ - }(t) = w_{m}^{ + }(t) - P(t){\text{/}}Z. \\ \end{gathered} $

Рассмотрим $(m + k)$ узел в предположении, что расчетная область была продолжена до него. С одной стороны, по формуле (5) имеем

$w_{{m + k}}^{ - } = ({{{v}}_{{m + k}}} - {{p}_{{m + k}}}{\text{/}}Z){\text{/}}2.$
С другой стороны, получаем
(9)
$\begin{gathered} w_{{m + k}}^{{ - (n)}} \equiv {{w}^{ - }}({{x}_{m}} + kh,{{t}^{{(n)}}}) = {{w}^{ - }}({{x}_{m}},{{t}^{{(n)}}} + kh{\text{/}}c) = \\ \mathop = \limits^{(8)} {{w}^{ + }}({{x}_{m}},{{t}^{{(n)}}} + kh{\text{/}}c) - P({{t}^{{(n)}}} + kh{\text{/}}c){\text{/}}Z = {{w}^{ + }}({{x}_{m}} - kh,{{t}^{{(n)}}}) - P({{t}^{{(n)}}} + kh{\text{/}}c){\text{/}}Z \equiv \\ \equiv w_{{m - k}}^{{ + (n)}} - P({{t}^{{(n)}}} + kh{\text{/}}c){\text{/}}Z = ({{{v}}_{{m - k}}} + {{p}_{{m - k}}}{\text{/}}Z){\text{/}}2 - P({{t}^{{(n)}}} + kh{\text{/}}c){\text{/}}Z. \\ \end{gathered} $
В формуле выше использовались: перенос значений вдоль характеристик, связь исходных функций и инвариантов Римана, заданное давление на границе $P(t)$ при $t \geqslant {{t}^{{(n)}}}$. Из двух предыдущих цепочек равенств получаем

(10)
$({{{v}}_{{m + k}}} - {{p}_{{m + k}}}{\text{/}}Z){\text{/}}2 = ({{{v}}_{{m - k}}} + {{p}_{{m - k}}}{\text{/}}Z){\text{/}}2 - P({{t}^{{(n)}}} + kh{\text{/}}c){\text{/}}Z.$

Аналогичные формулы можно выписать для $P(t)$, $t \leqslant {{t}^{{(n)}}}$:

$w_{{m + k}}^{ + } = ({{{v}}_{{m + k}}} + {{p}_{{m + k}}}{\text{/}}Z){\text{/}}2,$
(11)
$\begin{gathered} w_{{m + k}}^{{ + (n)}} \equiv {{w}^{ + }}({{x}_{m}} + kh,{{t}^{{(n)}}}) = {{w}^{ + }}({{x}_{m}},{{t}^{{(n)}}} - kh{\text{/}}c) = \\ \mathop = \limits^{(8)} {{w}^{ - }}({{x}_{m}},{{t}^{{(n)}}} - kh{\text{/}}c) + P({{t}^{{(n)}}} - kh{\text{/}}c){\text{/}}Z = {{w}^{ - }}({{x}_{m}} - kh,{{t}^{{(n)}}}) + P({{t}^{{(n)}}} - kh{\text{/}}c){\text{/}}Z \equiv \\ \equiv w_{{m - k}}^{{ - (n)}} + P({{t}^{{(n)}}} - kh{\text{/}}c){\text{/}}Z = ({{{v}}_{{m - k}}} - {{p}_{{m - k}}}{\text{/}}Z){\text{/}}2 + P({{t}^{{(n)}}} - kh{\text{/}}c){\text{/}}Z, \\ \end{gathered} $
(12)
$({{{v}}_{{m + k}}} + {{p}_{{m + k}}}{\text{/}}Z){\text{/}}2 = ({{{v}}_{{m - k}}} - {{p}_{{m - k}}}{\text{/}}Z){\text{/}}2 + P({{t}^{{(n)}}} - kh{\text{/}}c){\text{/}}Z.$

Полученные уравнения (10) и (12) задают следующую систему линейных уравнений на неизвестные в ghost-узлах ${{{v}}_{{m + k}}}$, ${{p}_{{m + k}}}$:

$\begin{gathered} ({{{v}}_{{m + k}}} + {{p}_{{m + k}}}{\text{/}}Z){\text{/}}2 = ({{{v}}_{{m - k}}} - {{p}_{{m - k}}}{\text{/}}Z){\text{/}}2 + P({{t}^{{(n)}}} - kh{\text{/}}c){\text{/}}Z, \\ ({{{v}}_{{m + k}}} - {{p}_{{m + k}}}{\text{/}}Z){\text{/}}2 = ({{{v}}_{{m - k}}} + {{p}_{{m - k}}}{\text{/}}Z){\text{/}}2 - P({{t}^{{(n)}}} + kh{\text{/}}c){\text{/}}Z. \\ \end{gathered} $
Из этой системы находим
(13)
$\begin{gathered} {{p}_{{m + k}}} = - {{p}_{{m - k}}} + P({{t}^{{(n)}}} + kh{\text{/}}c) + P({{t}^{{(n)}}} - kh{\text{/}}c), \\ {{{v}}_{{m + k}}} = {{{v}}_{{m - k}}} - [P({{t}^{{(n)}}} + kh{\text{/}}c) - P({{t}^{{(n)}}} - kh{\text{/}}c)]{\text{/}}Z. \\ \end{gathered} $
Эти формулы верны и для правой границы $k > 0$, и для левой границы $k < 0$ (т.е. ${{x}_{{m + k}}}$ всегда выходит за границу).

Можно сделать важное наблюдение, связанное с гладкостью рассматриваемых функций. Для корректного задания граничных условий достаточно задания значений только на характеристиках определенного знака. Например, значения на отрицательной характеристике для уравнения акустики. Следовательно, существует бесконечное множество способов продолжить функцию на ghost-узлы так, чтобы аналитическое решение задавало требуемое условие на границе. Однако существует единственный способ гладкого продолжения такой функции, удовлетворяющей исходным уравнениям. Формулы (13) задают именно это гладкое продолжение.

Эти формулы можно сравнить с другими известными формулами для филлера, представленными, например, в работе [25]. Для заданного давления ${{p}_{m}} = P(t)$ в работе [25] предлагаются заполнения вида

(14)
${{p}_{{m + k}}} = 2P({{t}^{n}}) - {{p}_{{m - k}}},\quad {{{v}}_{{m + k}}} = {{{v}}_{{m - k}}},$
которые отличаются от выведенных нами формул (13). Можно заметить, что формулы (14) задают разрыв в производной ${{{v}}_{x}}$. Метод конечных объемов успешно преодолевает эту проблему при явном вычислении потоков, но конечно-разностный и сеточно-характеристический методы в этом случае имеют лишь первый порядок аппроксимации. В то же время для “свободной границы” ($P(t) \equiv 0$) приведенные формулы (13) вырождаются в хорошо известные соотношения
${{p}_{{m + k}}} = - {{p}_{{m - k}}},\quad {{{v}}_{{m + k}}} = {{{v}}_{{m - k}}}.$
В этом случае разрыва производной ${{{v}}_{x}}$ в точке ${{x}_{m}}$ не происходит, так как $p$ и ${v}$ удовлетворяют уравнениям акустики и в граничной точке ${{p}_{m}} \equiv 0$ $ \Rightarrow $ ${{p}_{t}}({{x}_{m}},t) = 0$ $ \Rightarrow $ ${{{v}}_{x}}({{x}_{m}},t) = (cZ{{)}^{{ - 1}}}{{p}_{t}}({{x}_{m}},t) \equiv 0$.

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

Несмотря на то, что в данной работе рассматривались только уравнения акустики, обобщение на другие одномерные системы линейных гиперболических уравнений, по всей видимости, возможно. Должны быть рассмотрены инварианты Римана для всех характеристик, выходящих из расчетной области. В акустике это один инвариант: ${{w}^{ - }}$ для правой границы, ${{w}^{ + }}$ для левой. Для каждого из них будет использовано соответствующее собственное значение. Аналоги формул (5), (6) и (8) сохранят линейность. Они могут быть записаны в виде матриц, тогда как для уравнений акустики каждая из формул задается одним линейным уравнением. Полученные системы не вырождены и могут быть разрешены аналитически, приводя к формулам продления вида (13). Данная задача представляется интересной в плане дальнейшего развития выбранного направления исследований.

5. РЕАЛИЗАЦИЯ КОНТАКТНЫХ УСЛОВИЙ

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

1) описывать модели, в различных областях которых параметры среды имеют различные значения;

2) описывать области сложной геометрической формы, используя комбинации структурных сеток;

3) описывать динамическое поведение сред с различной реологией, находящихся в контакте.

Здесь и далее нижними индексами $a$ и $b$ будем обозначать величины, относящиеся к контактирующим средам. Линейные контактные условия записываются следующим образом в виде:

$C\left[ {\begin{array}{*{20}{c}} {{{{\mathbf{q}}}_{a}}} \\ {{{{\mathbf{q}}}_{b}}} \end{array}} \right] = {\mathbf{c}},$
где $C$ – прямоугольная матрица размерами $R \times ({{n}_{a}} + {{n}_{b}})$, ${{n}_{a}}$ – размер вектора ${{{\mathbf{q}}}_{a}}$, ${{n}_{b}}$ – размер вектора ${{{\mathbf{q}}}_{b}}$, $R$ – суммарное число выходящих характеристик для обеих сеток. Записанные в таком виде контактные условия могут быть реализованы с использованием рассмотренных ранее граничных условий.

Для контакта двух акустических сред необходимо задать два условия, так как в каждой среде одна характеристика выходит за границу области. Классическими условиями непроникновения являются равенство давлений и нормальных компонент скорости. В одномерном случае эти условия в точке ${{x}_{m}}$ на границе можно записать в виде

(15)
$\begin{gathered} {{p}_{a}}(t) = {{p}_{b}}(t), \\ {{{v}}_{a}}(t) = {{{v}}_{b}}(t). \\ \end{gathered} $

Для задания контактных условий (15) с помощью корректора граничных условий (см. п. 4.1, (7)) необходимо знать $P({{t}^{{n + 1}}})$, зависимость которого от ${\mathbf{q}}_{a}^{{{\text{in}}}}$, ${\mathbf{q}}_{b}^{{{\text{in}}}}$ может быть выражена аналитически.

Для задания контактных условий (15) с помощью предлагаемого подхода филлеров (см. п. 4.2, (13)) необходимо знать

$P_{a}^{k} = p({{t}^{n}} + k{{h}_{a}}{\text{/}}{{c}_{a}}),\quad P_{b}^{k} = p({{t}^{n}} + k{{h}_{b}}{\text{/}}{{c}_{b}}).$
В формуле выше $k$ определяет ghost-узел, в обозначениях разд. 2 $k = M$ – полуширина шаблона. Для сеточно-характеристических схем 3 -го и 4-го порядков точности мы использовали $M = 2$. Таким образом, $k \in \{ 1,2, - 1, - 2\} $.

Пусть среда $a$ слева, среда $b$ справа. Найдем $P_{a}^{k}$, $k > 0$.

Из контактных условий (15)

${{{v}}_{a}}({{t}^{n}} + k{{h}_{a}}{\text{/}}{{c}_{a}}) = {{{v}}_{b}}({{t}^{n}} + k{{h}_{a}}{\text{/}}{{c}_{a}}).$
Из связи исходных функций и инвариантов Римана (6)
$w_{a}^{ + }({{x}_{m}},{{t}^{n}} + k{{h}_{a}}{\text{/}}{{c}_{a}}) + w_{a}^{ - }({{x}_{m}},{{t}^{n}} + k{{h}_{a}}{\text{/}}{{c}_{a}}) = w_{b}^{ + }({{x}_{m}},{{t}^{n}} + k{{h}_{a}}{\text{/}}{{c}_{a}}) + w_{b}^{ - }({{x}_{m}},{{t}^{n}} + k{{h}_{a}}{\text{/}}{{c}_{a}}).$
Из формул (8) для $P_{a}^{k} = {{p}_{a}}({{t}^{n}} + k{{h}_{a}}{\text{/}}{{c}_{a}}) = {{p}_{b}}({{t}^{n}} + k{{h}_{a}}{\text{/}}{{c}_{a}})$ (знаем $w_{a}^{ + }$, $w_{b}^{ - }$; не знаем $w_{a}^{ - }$, $w_{b}^{ + }$)
$2w_{a}^{ + }({{x}_{m}},{{t}^{n}} + k{{h}_{a}}{\text{/}}{{c}_{a}}) - P_{a}^{k}{\text{/}}{{Z}_{a}} = 2w_{b}^{ - }({{x}_{m}},{{t}^{n}} + k{{h}_{a}}{\text{/}}{{c}_{a}}) + P_{a}^{k}{\text{/}}{{Z}_{b}}.$
Из постоянства значений инвариантов на характеристиках
$2w_{a}^{ + }({{x}_{m}} - k{{h}_{a}},{{t}^{n}}) - P_{a}^{k}{\text{/}}{{Z}_{a}} = 2w_{b}^{ - }\left( {{{x}_{m}} + k\frac{{{{с}_{b}}}}{{{{c}_{a}}}}{{h}_{a}},{{t}^{n}}} \right) + P_{a}^{k}{\text{/}}{{Z}_{b}}.$
Отсюда получаем итоговую формулу при $k > 0$:
(16)
$P_{a}^{k} = \frac{{2{{Z}_{a}}{{Z}_{b}}}}{{{{Z}_{a}} + {{Z}_{b}}}}\left[ {w_{{a,m - k}}^{{ + (n)}} - w_{b}^{{ - (n)}}\left( {{{x}_{m}} + k\frac{{{{c}_{b}}}}{{{{c}_{a}}}}{{h}_{a}}} \right)} \right],\quad k > 0.$
Аналогично, при $k < 0$ получим
(17)
$P_{a}^{k} = \frac{{2{{Z}_{a}}{{Z}_{b}}}}{{{{Z}_{a}} + {{Z}_{b}}}}\left[ {w_{b}^{{ + (n)}}\left( {{{x}_{m}} + \left| k \right|\frac{{{{c}_{b}}}}{{{{c}_{a}}}}{{h}_{a}}} \right) - w_{{a,m - |k|}}^{{ - (n)}}} \right],\quad k < 0.$
Для материала $b$ формулы принимают вид
(18)
$P_{b}^{k} = \frac{{2{{Z}_{a}}{{Z}_{b}}}}{{{{Z}_{a}} + {{Z}_{b}}}}\left[ {w_{a}^{{ + (n)}}\left( {{{x}_{m}} - k\frac{{{{c}_{a}}}}{{{{c}_{b}}}}{{h}_{b}}} \right) - w_{{b,m + k}}^{{ - (n)}}} \right],\quad k > 0,$
(19)
$P_{b}^{k} = \frac{{2{{Z}_{a}}{{Z}_{b}}}}{{{{Z}_{a}} + {{Z}_{b}}}}\left[ {w_{{b,m + |k|}}^{{ + (n)}} - w_{a}^{{ - (n)}}\left( {{{x}_{m}} - \left| k \right|\frac{{{{c}_{a}}}}{{{{c}_{b}}}}{{h}_{b}}} \right)} \right],\quad k < 0.$
В правых частях формул выше используются значения вида $w_{i}^{{ \pm (n)}}\left( {{{x}_{m}} - k\frac{{{{c}_{a}}}}{{{{c}_{b}}}}{{h}_{b}}} \right)$. Данные точки, в общем случае, не попадают в узлы расчетной сетки. Для их определения мы предлагаем использовать полиномиальную интерполяцию.

Таким образом, зная значения функций на текущем шаге по времени в обеих сетках, мы можем посчитать значения $P_{i}^{k}$. Затем они должны быть подставлены в филлеры (13) для заполнения ghost-узлов – своих для каждой сетки. После этого выполнение шага схемы во всех точках в обеих сетках обеспечит автоматическое выполнение контактных условий между двумя средами.

Фиг. 1.

Схематическое изображение шаблона.

6. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ

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

Авторами проведена серия расчетов с уменьшающимися в два раза шагом сетки $h$ и шагом по времени $dt$. Таким образом, число Куранта оставалось постоянным. Для каждого расчета была посчитана ошибка метода как норма разности численного и аналитического решений в финальный момент времени. Использовались интегральная норма ${{\left\| f \right\|}_{{{{L}_{1}}}}} = \sum\nolimits_i \left| {{{f}_{i}}} \right| \cdot h$ и более строгая норма ${{\left\| f \right\|}_{{{{L}_{\infty }}}}} = {{\max }_{i}}\left| {{{f}_{i}}} \right|$. Для каждой последовательной пары расчетов эмпирическая оценка порядка сходимости схемы была рассчитана как логарифм отношения ошибок $\mathcal{P} = {{\log }_{2}}({{E}_{{2h}}}{\text{/}}{{E}_{h}})$, где ${{E}_{h}}$ – ошибка метода при шаге сетки $h$.

6.1. Тестирование граничных условий

Рассматривалась акустическая задача (3) в одномерной постановке в области $x \in [0,1000]$ м, $t \in [0,T]$, $T = 0.5$ с. Параметры среды задавались следующими: $\rho = 1000$ кг/м3, $c = 1500$ м/с. Были заданы нулевые начальные условия ${\mathbf{q}}{{{\text{|}}}_{{t = 0}}} = {\mathbf{0}}$. На левой границе было задано давление ${{\left. {p(x,t)} \right|}_{{x = 0}}} = P(t)$, где $P(t) = {{\sin }^{4}}(\omega t)$ при $t > 0$ и $P(t) = 0$ при $t \leqslant 0$, $\omega = 2\pi f$, $f = 10$ Гц. Заданное таким образом давление на левой границе порождает волны, распространяющиеся от границы вправо. На правой границе области задавалось условие $p = 0$, но за расчетное время $T = 0.5$ с волна не достигала правой границы. Аналитическое решение задачи может быть записано в виде

$p(x,t) = \left( {\begin{array}{*{20}{l}} {0,\quad x - ct > 0,} \\ {\mathop {\sin }\nolimits^4 \left( {\frac{\omega }{c}(x - ct)} \right),\quad x - ct \leqslant 0.} \end{array}} \right.$

Для решения отдельных уравнений переноса использовалась сеточно-характеристическая схема 3 -го порядка аппроксимации. В табл. 1–3 приведены результаты численного исследования порядка сходимости схемы.

Таблица 1.  

Результаты численного исследования сходимости для задачи с граничными условиями. Расчет с использованием граничных корректоров (7)

Шаг сетки $h$ Ошибка по норме ${{L}_{1}}$ Ошибка по норме ${{L}_{\infty }}$ Порядок ${{\mathcal{P}}_{1}}$ Порядок ${{\mathcal{P}}_{\infty }}$
10.0000 132.2322 0.4012 – –
5.00000 69.892 0.2112 0.9199 0.9257
2.50000 34.2352 0.0991 1.0297 1.0911
1.25000 16.8945 0.0468 1.0189 1.0824
0.62500 8.4230 0.0230 1.0042 1.0238
0.31250 4.2091 0.0115 1.0008 1.0059
0.15625 2.1043 0.0057 1.0001 1.0015
Таблица 2.  

Результаты численного исследования сходимости для задачи с граничными условиями. Расчет с использованием формул (14), предложенных LeVeque [25]

Шаг сетки $h$ Ошибка по норме ${{L}_{1}}$ Ошибка по норме ${{L}_{\infty }}$ Порядок ${{\mathcal{P}}_{1}}$ Порядок ${{\mathcal{P}}_{\infty }}$
10.0000 162.0691 0.4381 – –
5.00000 92.5387 0.2573 0.8085 0.7678
2.50000 46.6821 0.1312 0.9872 0.9722
1.25000 23.2507 0.0640 1.0056 1.0358
0.62500 11.6126 0.0317 1.0016 1.0133
0.31250 5.8050 0.0158 1.0003 1.0037
0.15625 2.9024 0.0079 1.0000 1.0009
Таблица 3.  

Результаты численного исследования сходимости для задачи с граничными условиями. Предложенный в работе подход

Шаг сетки $h$ Ошибка по норме ${{L}_{1}}$ Ошибка по норме ${{L}_{\infty }}$ Порядок ${{\mathcal{P}}_{1}}$ Порядок ${{\mathcal{P}}_{\infty }}$
10.0000 93.4185 0.3694 – –
5.00000 32.2373 0.1400 1.5350 1.4001
2.50000 5.9058 0.0285 2.4485 2.2939
1.25000 0.7904 0.0039 2.9014 2.8877
0.62500 0.0997 0.0005 2.9873 2.9842
0.31250 0.0125 0.0001 2.9991 2.9979
0.15625 0.0016 0.0000 3.0002 2.9998

Условия “свободной границы” являются частным случаем условия заданного давления $P(t)$ при $P(t) \equiv 0$. Результаты тестирования порядка схемы на задаче о плоской волне, падающей на “свободную границу”, показали аналогичные результаты. Таким образом, предложенный в работе подход приводит к уменьшению ошибки метода и повышению порядка сходимости схемы.

6.2. Тестирование контактных условий

Рассматривалась задача контакта двух сред, обозначаемых нижними индексами $a$ и $b$. Параметры сред задавались следующими: ${{\rho }_{a}} = {{\rho }_{b}} = 1000$ кг/м3, ${{c}_{a}} = 1500$ м/с, ${{c}_{b}} = 2000$ м/с. Размеры областей, описывающих среды $a$ и $b$, были равны 300 и 400 м соответственно. В среде $a$, находящейся слева, начальные условия задавались в виде плоской волны давления с формой импульса ${{\sin }^{4}}\left( {2\pi fx} \right)$, $f = 0.01$ м–1. Была задана единичная амплитуда скорости; соответственно, амплитуда давления определялась значением импеданса ${{Z}_{a}} = {{c}_{a}}{{\rho }_{a}} = 1.5 \times {{10}^{6}}$. Данная волна движется вправо, достигает границы раздела сред и разделяется на прошедшую волну, продолжающую движение вправо в среде $b$, и отраженную волну, движущуюся обратно в среде $a$. Коэффициенты отражения ${{R}_{{{\text{reflection}}}}} = \frac{{{{Z}_{b}} - {{Z}_{a}}}}{{{{Z}_{b}} + {{Z}_{a}}}}$ и прохождения ${{R}_{{{\text{transmitted}}}}} = \frac{{2{{Z}_{b}}}}{{{{Z}_{a}} + {{Z}_{b}}}}$ определяют амплитуды данных волн. Для данной задачи известно аналитическое решение. Физическое время расчета было равно $T = 0.2$ с. Для решения отдельных уравнений переноса использовалась сеточно-характеристическая схема 3 -го порядка аппроксимации.

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

Таблица 4.  

Результаты численного исследования сходимости для задачи с контактом двух сеток. Использован контакт на основе граничных корректоров (7)

Шаг сетки $h$ Ошибка по норме ${{L}_{1}}$ Ошибка по норме ${{L}_{\infty }}$ Порядок ${{\mathcal{P}}_{1}}$ Порядок ${{\mathcal{P}}_{\infty }}$
10.0000 6.92e+07 8.81e+05 – –
5.00000 3.36e+07 5.35e+05 1.0414 0.7201
2.50000 1.14e+07 2.03e+05 1.5671 1.3991
1.25000 2.18e+06 4.21e+04 2.3802 2.2684
0.62500 5.77e+05 1.16e+04 1.9189 1.8564
0.31250 2.53e+05 4.66e+03 1.1911 1.3191
0.15625 1.24e+05 2.22e+03 1.0247 1.0729
Таблица 5.  

Результаты численного исследования сходимости для задачи с контактом двух сеток. Использован предложенный в работе подход для постановки контактных условий

Шаг сетки $h$ Ошибка по норме ${{L}_{1}}$ Ошибка по норме ${{L}_{\infty }}$ Порядок ${{\mathcal{P}}_{1}}$ Порядок ${{\mathcal{P}}_{\infty }}$
10.0000 6.59e+07 7.63e+05 – –
5.00000 2.63e+07 3.74e+05 1.3273 1.0287
2.50000 7.07e+06 9.93e+04 1.8930 1.9124
1.25000 1.08e+06 1.51e+04 2.7143 2.7221
0.62500 1.39e+05 1.94e+03 2.9500 2.9594
0.31250 1.75e+04 2.43e+02 2.9929 2.9907
0.15625 2.19e+03 3.05e+01 2.9989 2.9978

6.3. Расчет с использованием множества сеток с одинаковыми параметрами

Иногда возникает необходимость описывать расчетную область несколькими расчетными сетками, даже если параметры среды в каждой из подобластей-сеток являются одинаковыми. Мы провели тестовый расчет такой задачи в одномерном случае. Начальные условия задают плоскую волну с формой импульса ${{\sin }^{4}}\left( {2\pi fx} \right)$, $f = 0.01$ м–1, и единичной амплитудой скорости (амплитуда давления равна $c \cdot \rho $, $c = 1500$ м/с, $\rho = 1000$ кг/м3). Размер расчетной области был равен 600 м. Аналитическое решение в момент времени $T = 0.3$ с представлено на фиг. 2. Видно, что численное решение cеточно-характеристической схемой третьего порядка аппроксимации с использованием одной расчетной сетки с шагом $h = 2.5$ м несколько отличается от аналитического решения. Это объясняется наличием у схемы численной вязкости. Во втором расчете область была покрыта десятью сетками с одинаковыми параметрами $c$, $\rho $; при этом в расчете явно выделялись все контакты между сетками. На фиг. 2 приведены результаты, полученные схемой с граничными корректорами (7) и новым предложенным в работе способом. Видно, что первый подход приводит к заметной потере волной амплитуды. При этом решение по новой предложенной схеме визуально неотличимо от решения, полученного с использованием одной единой сетки.

Фиг. 2.

Сравнение аналитического решения с численными, полученными различными способами: без выделения контактов (а), с выделением контактов при помощи корректоров (7) (б), с выделением контактов при помощи заполнения ghost-узлов (16) (в).

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

В работе рассмотрены вопросы постановки граничных и контактных условий при численном решении динамических задачах деформируемых твердых тел. Рассмотрено семейство сеточно-характеристических методов на структурных расчетных сетках, которые успешно применяются для решения систем линейных гиперболических уравнений. Использование сеточно-характеристических схем повышенного порядка аппроксимации позволяет построить расчетный алгоритм с порядком сходимости выше первого для внутренних точек расчетной области, в том числе и для многомерных задач. Однако важной задачей является сохранение повышенного порядка и на границах области. В работе предложена процедура заполнения ghost-узлов, обеспечивающая достижение данного свойства. На примере одномерной акустической задачи подтверждены снижение ошибки метода и повышение порядка сходимости предложенных схем для задач с заданным условием на границе и с наличием явно выделенной контактной границы. Описанный в работе подход представляется возможным расширить на общий класс линейных гиперболических систем уравнений. Важным вопросом является исследование возможности его обобщения на многомерный случай с использованием метода расщепления по пространственным направлениям. Развиваемые в настоящей работе сеточно-характеристические методы могут быть полезны при решении практических прямых и обратных задач сейсмической разведки [26].

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

  1. Zhou H., Liu Y., Wang J. Elastic Wave Modeling With High-Order Temporal and Spatial Accuracies by a Selectively Modified and Linearly Optimized Staggered-Grid Finite-Difference Scheme // IEEE Transactions on Geoscience and Remote Sensing. 2022. V. 60. P. 1–22.

  2. Reinarz A., Charrier D.E., Bader M., Bovard L., Dumbser M., Duru K., Fambri F., Gabriel A.-A., Gallard J.-M., Koppel S., Krenz L., Rannabauer L., Rezzolla L., Samfass P., Tavelli M., Weinzierl T. ExaHyPE: An engine for parallel dynamically adaptive simulations of wave problems // Computer Physics Communications. 2020. V. 254.

  3. Магомедов К.М., Холодов А.С. Сеточно-характеристические численные методы: учебное пособие для вузов. М: Юрайт, 2023.

  4. Yang H.Q., Harris R.E. Development of Vertex-Centered High-Order Schemes and Implementation in FUN3D // AIAA Journal. 2016. V. 54. I. 12. P. 1–19.

  5. Van Leer B., Nishikawa H. Towards the ultimate understanding of MUSCL: Pitfalls in achieving third-order accuracy // J. Comput. Phys. 2021. V. 446. P. 110640.

  6. Nishikawa H. On False Accuracy Verification of UMUSCL Scheme // Communications in Comput. Physics. 2021. V. 30. I. 4. P. 1037–1060.

  7. Padway E., Nishikawa H. Resolving Confusions over Third-Order Accuracy of Unstructured MUSCL // AIAA Journal. 2022. V. 60. I. 3. P. 1415–1439.

  8. Nishikawa H. Economically high-order unstructured-grid methods: Clarification and efficient FSR schemes // Intern. Journal for Numerical Methods in Fluids. 2021. V. 93. I. 11. P. 3187–3214.

  9. Nishikawa H., Van Leer B. Towards High-Order Boundary Procedures for Finite-Volume and Finite-Difference Schemes // AIAA SCITECH 2023 Forum.

  10. Ладонкина М.Е., Неклюдова О.А., Тишкин В.Ф. Исследование влияния лимитера на порядок точности решения разрывным методом Галеркина // Матем. моделирование. 2012. Т. 24. № 12. С. 124–128.

  11. Guseva E.K., Golubev V.I., Petrov I.B. Linear, Quasi-Monotonic and Hybrid Grid-Characteristic Schemes for Hyperbolic Equations // Lobachevskii Journal of Mathematics. 2023. V. 44. I. 1. P. 296–312.

  12. Golubev V.I., Muratov M.V., Guseva E.K., Konov D.S., Petrov I.B. Thermodynamic and Mechanical Problems of Ice Formations: Numerical Simulation Results // Lobachevskii Journal of Mathematics. 2022. V. 43 I. 4.

  13. Guseva E.K., Beklemysheva K.A., Golubev V.I., Epifanov V.P., Petrov I.B. Investigation of Ice Rheology Based on Computer Simulation of Low-Speed Impact // Mathematical Modeling and Supercomputer Technologies. 2022. P. 176-184.

  14. Petrov I.B., Golubev V.I., Shevchenko A.V. Higher-Order Grid-Characteristic Schemes for the Acoustic System // 2021 Ivannikov Memorial Workshop (IVMEM). 2021. P. 61–65.

  15. Golubev V.I., Shevchenko A.V., Petrov I.B. Raising convergence order of grid-characteristic schemes for 2D linear elasticity problems using operator splitting // Computer Research and Modeling. 2022. V. 14. I. 4. P. 899–910.

  16. Golubev V.I., Shevchenko A.V., Khokhlov N.I., Petrov I.B., Malovichko M.S. Compact Grid-Characteristic Scheme for the Acoustic System with the Piece-Wise Constant Coefficients // International Journal of Applied Mechanics. 2022. V. 14. I. 2. P. 2250002.

  17. Petrov I.B., Kholodov A.S. Numerical study of some dynamic problems of the mechanics of a deformable rigid body by the mesh-characteristic method // Computational Mathematics and Mathematical Physics. 1984. V. 24. I. 3. P. 61–73.

  18. Петров И.Б., Тормасов А.Г., Холодов А.С. Об использовании гибридизированных сеточно-характеристических схем для численного решения трехмерных задач динамики деформируемого твердого тела // Ж. вычисл. матем. и матем. физ. 1990. Т. 30. № 8. С. 1237–1244.

  19. 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 // Geophysical Prospecting. 2018. V. 66. I. 8. P. 1485–1502.

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

  21. Исакович М.А. Общая акустика. М.: Наука, 1973.

  22. Sofronov I., Zaitsev N., Dovgilovich L. Multi-block finite-difference method for 3D elastodynamic simulations in anisotropic subhorizontally layered media // Geophysical Prospecting. 2015. V. 63. P. 1142–1160.

  23. Komatitsch D., Tromp J. A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation // Geophysical Journal International. 2003. V. 154. I. 1. P. 146–153.

  24. Челноков Ф.Б. Численное моделирование деформационных динамических процессов в средах со сложной структурой: Дис. … канд. физ. -матем. наук. М.: Моск. физ.-техн. ин-т (гос. ун-т), 2005.

  25. LeVeque R.J. Finite Volume Methods for Hyperbolic Problems. Cambridge Texts in Applied Mathematics. Cambridge University Press, 2002.

  26. Golubev V.I., Nikitin I.S., Vasyukov A.V., Nikitin A.D. Fractured inclusion localization and characterization based on deep convolutional neural networks // Procedia Structural Integrity. 2023. V. 43. P. 29–34.

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

Инструменты

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

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

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