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

Журнал вычислительной математики и математической физики, 2021, T. 61, № 8, стр. 1363-1377

Расчет сейсмостойкости ледового острова сеточно-характеристическим методом на комбинированных расчетных сетках

И. Б. Петров 1, А. В. Фаворская 1, *

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

* E-mail: aleanera@yandex.ru

Поступила в редакцию 24.11.2020
После доработки 24.12.2020
Принята к публикации 14.01.2021

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

Аннотация

Целью работы являлись разработка и применение численных методов из семейства сеточно-характеристических, позволяющих рассчитать повреждения оффшорного объекта (ледового острова), рассчитывая распространение сейсмических волн от гипоцентра землетрясения, расположенного на глубине нескольких километров с учетом специфики cеверных морей (небольшая глубина). Для этого был применен сеточно-характеристический метод на комбинированных расчетных сетках. В одном расчете использовалось свыше 25 отдельных расчетных сеток, отличающихся друг от друга решаемой системой уравнений, типом расчетной схемы, шагами интегрирования по координатам. Часть этих расчетных сеток представляли собой систему из вложенных друг в друга иерархических сеток с кратным шагом по координате, остальные являются конформными друг относительно друга. Также использовалось улучшенное неотражающее условие, предполагающее резкое увеличение шага по координате и применение диссипативной разностной схемы. Особое внимание в работе уделено описанию алгоритмов для расчета на границах и контактных границах этих отдельных расчетных сеток, позволяющих реализовать на практике предложенный подход с использованием сеточно-характеристического метода на комбинированных расчетных сетках. Библ. 40. Фиг. 13. Табл. 3.

Ключевые слова: сеточно-характеристический метод, численное моделирование, ледовый остров, арктические задачи, сейсмостойкость, комбинированные расчетные сетки, иерархические расчетные сетки.

1. ВВЕДЕНИЕ

Ледовые острова (искусственно возведенные) – это, зачастую, единственный способ провести разведочное бурение в мелководных районах cеверных морей, куда из-за суровой ледовой обстановки невозможно доставить обычную платформу. Данный подход успешно реализован в Канаде (см. [1], [2]).

Возведение ледового острова – это дорогостоящий процесс, требующий предварительных расчетов, в том числе для оценки их сейсмостойкости. В современных программных пакетах для оценки сейсмостойкости используют в основном метод конечных элементов (см. [3]–[5]). А для распространения сейсмических волн от гипоцентра землетрясения применяют конечно-разностный метод (см. [6]–[9]), метод спектральных элементов (см. [10]–[12]), разрывный метод Галеркина (см. [13]–[16]) и сеточно-характеристический численный метод (см. [17], [18]).

2. ПОСТАНОВКА ЗАДАЧИ

Рассматриваемая геологическая модель и ее геометрические параметры представлены на фиг. 1 (не в масштабе), соответствие рассматриваемым средам приведено на фиг. 2. Математические формулировки отмеченных на фиг. 1 граничных условий приведены в разд. 4.

Фиг. 1.

Геологическая модель.

Фиг. 2.

Используемые обозначения: 1 – ледовый остров, 2 – осадочные породы, 3 – гранитный слой, 4 – базальтовый слой, 5 – верхняя мантия, 6 – вода, 7 – неотражающие граничные условия, 8 – контактное условие полного слипания, 9 – свободная граница, 10 – нулевое давление, 11 – свободное скольжение, 12 – контактное условие между упругой и акустической средами, 13 – контактное условие между иерархическими вложенными сетками с кратным шагом по координате, 14 – контактное условие между подобластями с различными расчетными схемами с кратным шагом по координате.

Упругие и акустические параметры рассматриваемых сред приведены в табл. 1: плотность, скорости продольных (pressure, P-) волн и поперечных (shear, S-) волн.

Таблица 1.  

Упругие и акустические параметры рассматриваемых сред

Среда Скорость
P-волн, м/c 		Скорость
S-волн, м/c 		Плотность, кг/м3
Ледовый остров 3940 2493 917
Вода 1500 – 1025
Осадочные породы 2250 1000 2000
Гранитный слой 5700 2500 2600
Базальтовый слой 6800 3000 3000
Верхняя мантия 8000 3500 3300

Расчетная область разделена на отдельные расчетные сетки, представленные на фиг. 2. Между ними установлены соответствующие граничные и контактные условия, отмеченные на фиг. 2 различными линиями. Цветом отмечены используемые граничные и контактные условия (фиг. 3).

Фиг. 3.

Отдельные расчетные сетки. Граничные и контактные условия.

Также ледовый остров окружен системой вложенных иерархических сеток (фиг. 4).

Фиг. 4.

Система вложенных иерархических сеток: (а) – без уровня, (б) – уровень –1, (в) – уровень –2, (г) – уровень 0, (д) – уровень 1, (е) – уровень N – 2, (ж) – уровень N – 1.

В табл. 2 приведены параметры всех отдельных расчетных сеток.

Таблица 2.  

Параметры расчетных сеток, пример

Среда Уро-вень Номер сетки Разностная схема Общее количе-ство точек Шаг по оси OX, м Шаг по оси OY, м Размер вдоль оси OX, км Размер вдоль оси OY, км
Вода – 1 КИР 8844 125 0.4 50 0.008
Вода – 2 КИР 8844 125 0.4 50 0.008
Осадочные породы – 3 КИР 129 444 125 6.25 50 2
Осадочные породы – 4 КИР 129 444 125 6.25 50 2
Гранитный слой – 5 КИР 161 604 125 12.5 50 5
Гранитный слой – 6 КИР 161 604 125 12.5 50 5
Базальтовый слой – 7 КИР 290 244 125 25 50 18
Базальтовый слой – 8 КИР 290 244 125 25 50 18
Верхняя мантия – 9 КИР 161 604 125 125 50 50
Верхняя мантия – 10 КИР 161 604 125 125 50 50
Вода – 11 Схема Русанова 50 669 25 0.4 55 0.008
Осадочные породы – 12 Схема Русанова 711 569 25 6.25 55 2
Гранитный слой – 13 Схема Русанова 887 809 25 12.5 55 5
Базальтовый слой – 14 Схема Русанова 1 592 769 25 25 55 18
Верхняя мантия – 15 КИР 885 204 25 125 55 50
Вода –2 I Схема Русанова 2415 12.5 0.4 1.275 0.008
Осадочные породы –2 II Схема Русанова 33 915 12.5 6.25 1.275 2
Гранитный слой –2 III Схема Русанова 2415 12.5 12.5 1.275 0.25
Вода –1 I Схема Русанова 2921 6.25 0.4 0.775 0.008
Осадочные породы –1 II Схема Русанова 5207 6.25 6.25 0.775 0.2375
Вода 0 I Схема Русанова 3933 3.125 0.4 0.525 0.008
Осадочные породы 0 II Схема Русанова 6669 3.125 3.125 0.525 0.1125
Вода 1 I Схема Русанова 5865 1.5625 0.4 0.39375 0.008
Осадочные породы 1 II Схема Русанова 8415 1.5625 1.5625 0.39375 0.046875
Вода 2 I Схема Русанова 529 0.78125 0.4 0.015625 0.008
Вода 2 II Схема Русанова 529 0.78125 0.4 0.015625 0.008
Осадочные породы 2 III Схема Русанова 9821 0.78125 0.78125 0.33125 0.015625
Ледовый остров 2 IV Схема Русанова 10 836 0.78125 0.4 0.3 0.01

Примечание. Схема КИР – схема Куранта–Изаксона–Риса.

В каждой из отдельных расчетных сеток в зависимости от рассматриваемой среды решают либо акустическое волновое уравнение (в водном слое) (см. [19]):

(2.1)
${{\rho }}\frac{\partial }{{\partial t}}{\mathbf{v}}\left( {{\mathbf{r}},t} \right) = - \nabla p\left( {{\mathbf{r}},t} \right),$
(2.2)
$\frac{\partial }{{\partial t}}p\left( {{\mathbf{r}},t} \right) = - {{\rho }}{{c}^{2}}\left( {\nabla \cdot {\mathbf{v}}\left( {{\mathbf{r}},t} \right)} \right),$
либо упругое (в остальных слоях):

(2.3)
${{\rho }}\frac{\partial }{{\partial t}}{\mathbf{v}}\left( {{\mathbf{r}},t} \right) = {{\left( {\nabla \cdot {\mathbf{\sigma }}\left( {{\mathbf{r}},t} \right)} \right)}^{{\text{т}}}}$,
(2.4)
$\frac{\partial }{{\partial t}}{\mathbf{\sigma }}\left( {{\mathbf{r}},t} \right) = \left( {{{\rho }}c_{{\text{P}}}^{{\text{2}}} - 2{{\rho }}c_{{\text{S}}}^{{\text{2}}}} \right)\left( {\nabla \cdot {\mathbf{v}}\left( {{\mathbf{r}},t} \right)} \right)I + {{\rho }}c_{{\text{S}}}^{{\text{2}}}\left( {\nabla \otimes {\mathbf{v}}\left( {{\mathbf{r}},t} \right) + {{{\left( {\nabla \otimes {\mathbf{v}}\left( {{\mathbf{r}},t} \right)} \right)}}^{{\text{т}}}}} \right)$.

Здесь и далее ${\mathbf{v}}\left( {{\mathbf{r}},t} \right)$ – векторное поле скорости (производная смещения по времени), $p\left( {{\mathbf{r}},t} \right)$ – скалярное поле давления, ${\mathbf{\sigma }}\left( {{\mathbf{r}},t} \right)$ – тензорное поле симметричного тензора напряжений второго ранга Коши, с, ${{c}_{{\text{P}}}}$, ${{c}_{{\text{S}}}}$ – скорость звука в акустической среде, скорости продольной (P-) и поперечной (S-) волн в упругой среде соответственно, ${{\rho }}$ – плотность, $I$ – единичный тензор второго ранга, $ \otimes $ – знак тензорного произведения векторов ${{\left( {{\mathbf{a}} \otimes {\mathbf{b}}} \right)}_{{ij}}} = {{a}_{i}}{{b}_{j}}$. Для моделирования разрушения ледового острова использовались критерий по главному напряжению и модель трещин (см. [20]). Критическое главное напряжение бралось равным 1 МПа (см. [21]–[25]).

В соответствии с условиями устойчивости (см. [26]) шаг по времени брался равным 0.1 мс, общее время расчета – 24 с, что соответствует 240 000 шагам по времени. Заметим, что при определении шага по времени в соответствии с условиями устойчивости для каждой из расчетных сеток следует определить выражение для случая упругого волнового уравнения

$\frac{{\min \left( {{{h}_{{\text{X}}}},{{h}_{{\text{Y}}}}} \right)}}{{{{c}_{{\text{P}}}}}}$
и следующее выражение для случая акустического волнового уравнения:

$\frac{{\min \left( {{{h}_{{\text{X}}}},{{h}_{{\text{Y}}}}} \right)}}{c}.$

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

$k_{A}^{i} = \frac{{h_{A}^{i}}}{{h_{{\text{Y}}}^{{{\text{ICE}}}}}}\frac{{c_{{\text{P}}}^{{{\text{ICE}}}}}}{{c_{{\text{P}}}^{i}}}k_{{\text{Y}}}^{{{\text{ICE}}}}.$

Здесь $k_{{\text{Y}}}^{{{\text{ICE}}}}$ – число Куранта в направлении OY для сетки, покрывающей ледовый остров; $c_{{\text{P}}}^{{{\text{ICE}}}}$ – скорость продольных волн во льде; $h_{{\text{Y}}}^{{{\text{ICE}}}}$ – шаг вдоль оси OY в расчетной сетке, покрывающей ледовый остров; $A \in \left\{ {{\text{X}},{\text{Y}}} \right\}$; i – индекс рассматриваемой расчетной сетки; $h_{{\text{X}}}^{i}$, $h_{{\text{Y}}}^{i}$ – шаги по координате рассматриваемой расчетной сетки; $c_{{\text{P}}}^{i}$ – скорость продольных волн в рассматриваемой расчетной сетке.

В качестве источника используется модель очага землетрясения, задающаяся как начальное условие, примененная в [27], [28], c максимальной скоростью 20 м/c. Возможно использование более сложных моделей очагов землетрясений (см. [6]).

3. СЕТОЧНО-ХАРАКТЕРИСТИЧЕСКИЙ МЕТОД

3.1. Акустическое волновое уравнение

Осуществляется переход к инвариантам Римана, используя следующие выражения (на примере оси OX):

${{{{\omega }}}_{{1,2}}} = {{{v}}_{{\text{X}}}} \pm \frac{p}{{c{{\rho }}}},$
${{{{\omega }}}_{3}} = {{{v}}_{{\text{Y}}}}.$

Здесь и далее индексы X, Y соответствуют компонентам вектора скорости.

После перехода к инвариантам Римана решается система независимых уравнений переноса с коэффициентами $\left\{ {c, - c,0} \right\}$ соответственно.

Затем осуществляется обратный переход, используя выражения

${\mathbf{v}} = \left[ {\begin{array}{*{20}{c}} {\frac{1}{2}\left( {{{{{\omega }}}_{1}} + {{{{\omega }}}_{2}}} \right)} \\ {{{{{\omega }}}_{3}}} \end{array}} \right],$
$p = \frac{1}{2}c{{\rho }}\left( {{{{{\omega }}}_{1}} - {{{{\omega }}}_{2}}} \right).$

3.2. Упругое волновое уравнение

Осуществляется переход к инвариантам Римана, используя следующие выражения (на примере оси OX):

${{{{\omega }}}_{{1,2}}} = {{{v}}_{{\text{X}}}} \mp \frac{1}{{{{c}_{{\text{P}}}}{{\rho }}}}{{{{\sigma }}}_{{{\text{XX}}}}},$
${{{{\omega }}}_{{3,4}}} = {{{v}}_{{\text{Y}}}} \mp \frac{1}{{{{c}_{{\text{S}}}}{{\rho }}}}{{{{\sigma }}}_{{{\text{XY}}}}},$
${{{{\omega }}}_{5}} = {{{{\sigma }}}_{{{\text{YY}}}}} - \frac{{{{c}_{3}}}}{{{{c}_{{\text{P}}}}}}{{{{\sigma }}}_{{{\text{XX}}}}}.$

Здесь и далее индексы XX, YY и XY соответствуют компонентам симметричного тензора напряжений Коши, а ${{c}_{3}}$ задается выражением

${{c}_{3}} = \left( {1 - 2{{{\left( {\frac{{{{c}_{{\text{S}}}}}}{{{{c}_{{\text{P}}}}}}} \right)}}^{2}}} \right){{c}_{{\text{P}}}}.$

После перехода к инвариантам Римана решается система независимых уравнений переноса с коэффициентами $\left\{ {{{c}_{{\text{P}}}}, - {{c}_{{\text{P}}}},{{c}_{{\text{S}}}}, - {{c}_{{\text{S}}}},0} \right\}$ соответственно.

Затем осуществляется обратный переход, используя выражения

${\mathbf{v}} = \frac{1}{2}\left[ {\begin{array}{*{20}{c}} {{{{{\omega }}}_{1}} + {{{{\omega }}}_{2}}} \\ {{{{{\omega }}}_{3}} + {{{{\omega }}}_{4}}} \end{array}} \right],$
${\mathbf{\sigma }} = \frac{1}{2}\left[ {\begin{array}{*{20}{c}} {{{\rho }}{{c}_{{\text{P}}}}\left( {{{{{\omega }}}_{2}} - {{{{\omega }}}_{1}}} \right)}&{{{\rho }}{{с}_{{\text{S}}}}\left( {{{{{\omega }}}_{4}} - {{{{\omega }}}_{3}}} \right)} \\ {{{\rho }}{{с}_{{\text{S}}}}\left( {{{{{\omega }}}_{4}} - {{{{\omega }}}_{3}}} \right)}&{{{\rho }}{{с}_{3}}\left( {{{{{\omega }}}_{2}} - {{{{\omega }}}_{1}}} \right) + {{{{\omega }}}_{5}}} \end{array}} \right].$

4. ГРАНИЧНЫЕ И КОНТАКТНЫЕ УСЛОВИЯ

Так как шаблоны используемых разностных схем выходят за пределы области интегрирования, то вдоль границы располагаются слои дополнительных точек в соответствии с фиг. 5. На фиг. 5б, г и далее эти дополнительные точки отмечены крестиками.

Фиг. 5.

Шаблоны разностных схем (случай положительного коэффициента) и расположение дополнительных точек вдоль левой границы: (а) – шаблон разностной схемы Русанова; (б) – расположение дополнительных точек, схема Русанова; (в) – шаблон разностной схемы КИР; (г) – расположение дополнительных точек, схема КИР.

Вдоль каждой из координатных осей вычисления проводятся в три этапа.

1. Заполнение значений в дополнительных точках.

2. Расчет во внутренних узлах в соответствии с формулами из разд. 3. Значения в граничных точках, рассчитанные на данном этапе, далее в тексте отмечают индексом IN.

3. Действие граничных и контактных корректоров. Значения в граничных точках, рассчитанные на данном этапе, далее в тексте отмечают индексом $n + 1$.

Формулы в данном разделе приведены таким образом, что будут справедливы как для 2D, так и для 3D случая, а также как для регулярных структурированных сеток (см. [27]–[30]), так и для структурированных криволинейных (см. [17], [30], [31]).

Приведенные в данном разделе формулы для коррекции значений в граничных узлах выводятся из соображений восстановления инвариантов Римана (соответствующих выходящим за пределы рассматриваемой сетки характеристик) таким образом, чтобы выполнялaсь математическая формулировка граничного или контактного условия. Подробнее с алгоритмом вывода граничных и контактных условий можно ознакомиться в [32], [33].

4.1. Неотражающие граничные условия

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

(4.1)
${\mathbf{v}}\left( {{{{\mathbf{r}}}_{{m - 2}}},{{t}_{n}}} \right) = {\mathbf{v}}\left( {{{{\mathbf{r}}}_{{m - 1}}},{{t}_{n}}} \right) = {\mathbf{v}}\left( {{{{\mathbf{r}}}_{m}},{{t}_{n}}} \right),$
(4.2)
$p\left( {{{{\mathbf{r}}}_{{m - 2}}},{{t}_{n}}} \right) = p\left( {{{{\mathbf{r}}}_{{m - 1}}},{{t}_{n}}} \right) = p\left( {{{{\mathbf{r}}}_{m}},{{t}_{n}}} \right),$
(4.3)
${\mathbf{\sigma }}\left( {{{{\mathbf{r}}}_{{m - 2}}},{{t}_{n}}} \right) = {\mathbf{\sigma }}\left( {{{{\mathbf{r}}}_{{m - 1}}},{{t}_{n}}} \right) = {\mathbf{\sigma }}\left( {{{{\mathbf{r}}}_{m}},{{t}_{n}}} \right).$

Здесь нижние индексы соответствуют индексам по координате и индексам по времени в соответствии с фиг. 5. Выражения (4.1), (4.2) используют в случае акустического волнового уравнения (2.1), (2.2), а выражения (4.1), (4.3) – в случае упругого волнового уравнения (2.3), (2.4).

4.2. Контактное условие полного слипания, случай упругого волнового уравнения

Данное контактное условие записывается следующим образом:

(4.4)
${{{\mathbf{v}}}^{{\text{L}}}} = {{{\mathbf{v}}}^{{\text{R}}}},$
(4.5)
${{{\mathbf{\sigma }}}^{{\text{L}}}} \cdot {\mathbf{m}} = {{{\mathbf{\sigma }}}^{{\text{R}}}} \cdot {\mathbf{m}}.$

Здесь и далее m – внешняя нормаль к левой подобласти интегрирования (когда речь идет о контактных условиях), верхние индексы L и R соответствуют левой и правой контактирующим расчетным сеткам.

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

(4.6)
$\begin{gathered} {\mathbf{V}} = \frac{1}{{{{{{\rho }}}^{{\text{L}}}}c_{{\text{S}}}^{{\text{L}}} + {{{{\rho }}}^{{\text{R}}}}c_{{\text{S}}}^{{\text{R}}}}}\left\{ {\mathop {{{{{\rho }}}^{{\text{L}}}}}\limits_{}^{} \left( {c_{{\text{P}}}^{{\text{L}}} - c_{{\text{S}}}^{{\text{L}}}} \right)\left( {{\mathbf{m}} \cdot {\mathbf{v}}_{{{\text{IN}}}}^{{\text{L}}}} \right){\mathbf{m}} + {{{{\rho }}}^{{\text{L}}}}c_{{\text{S}}}^{{\text{L}}}{\mathbf{v}}_{{{\text{IN}}}}^{{\text{L}}} + } \right. \\ + \;{{{{\rho }}}^{{\text{R}}}}\left( {c_{{\text{P}}}^{{\text{R}}} - c_{{\text{S}}}^{{\text{R}}}} \right)\left( {{\mathbf{m}} \cdot {\mathbf{v}}_{{{\text{IN}}}}^{{\text{R}}}} \right){\mathbf{m}} + {{{{\rho }}}^{{\text{R}}}}c_{{\text{S}}}^{{\text{R}}}{\mathbf{v}}_{{{\text{IN}}}}^{{\text{R}}} + \left( {{\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{R}}} - {\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{L}}}} \right) \cdot {\mathbf{m}} - \\ \left. { - \;\frac{{{{{{\rho }}}^{{\text{L}}}}\left( {c_{{\text{P}}}^{{\text{L}}} - c_{{\text{S}}}^{{\text{L}}}} \right) + {{{{\rho }}}^{{\text{R}}}}\left( {c_{P}^{R} - c_{{\text{S}}}^{{\text{R}}}} \right)}}{{{{{{\rho }}}^{{\text{L}}}}c_{{\text{P}}}^{{\text{L}}} + {{{{\rho }}}^{{\text{R}}}}c_{{\text{P}}}^{{\text{R}}}}} \cdot \left( {{\mathbf{m}} \cdot \left( {\left( {{{{{\rho }}}^{{\text{L}}}}c_{{\text{P}}}^{{\text{L}}}{\mathbf{v}}_{{{\text{IN}}}}^{{\text{L}}} + {{{{\rho }}}^{{\text{R}}}}c_{{\text{P}}}^{{\text{R}}}{\mathbf{v}}_{{{\text{IN}}}}^{{\text{R}}}} \right) + \left( {{\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{R}}} - {\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{L}}}} \right) \cdot {\mathbf{m}}} \right)} \right){\mathbf{m}}} \right\}, \\ \end{gathered} $
(4.7)
${\mathbf{v}}_{{n + 1}}^{{\text{L}}} = {\mathbf{v}}_{{n + 1}}^{{\text{R}}} = {\mathbf{V}},$
(4.8)
$\begin{gathered} {\mathbf{\sigma }}_{{n + 1}}^{{\text{L}}} = {\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{L}}} + {{{{\rho }}}^{{\text{L}}}}\left\{ {\left( {\left( {{\mathbf{V}} - {\mathbf{v}}_{{{\text{IN}}}}^{{\text{L}}}} \right) \cdot {\mathbf{m}}} \right)\left( {\left( {c_{{\text{P}}}^{{\text{L}}} - 2c_{{\text{S}}}^{{\text{L}}} - c_{3}^{{\text{L}}}} \right)\left( {{\mathbf{m}} \otimes {\mathbf{m}}} \right) + c_{3}^{{\text{L}}}{\mathbf{I}}} \right)} \right. + \\ \left. { + \;c_{{\text{S}}}^{{\text{L}}}\left( {{\mathbf{m}} \otimes \left( {{\mathbf{V}} - {\mathbf{v}}_{{{\text{IN}}}}^{{\text{L}}}} \right) + \left( {{\mathbf{V}} - {\mathbf{v}}_{{{\text{IN}}}}^{{\text{L}}}} \right) \otimes {\mathbf{m}}} \right)} \right\}, \\ \end{gathered} $
(4.9)
$\begin{gathered} {\mathbf{\sigma }}_{{n + 1}}^{{\text{R}}} = {\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{R}}} + {{{{\rho }}}^{{\text{R}}}}\left\{ {\left( {\left( {{\mathbf{v}}_{{{\text{IN}}}}^{{\text{R}}} - {\mathbf{V}}} \right) \cdot {\mathbf{m}}} \right)\left( {\left( {c_{{\text{P}}}^{{\text{R}}} - 2c_{{\text{S}}}^{{\text{R}}} - c_{3}^{{\text{R}}}} \right)\left( {{\mathbf{m}} \otimes {\mathbf{m}}} \right) + c_{3}^{{\text{R}}}{\mathbf{I}}} \right)} \right. + \\ \left. { + \;c_{{\text{S}}}^{{\text{R}}}\left( {{\mathbf{m}} \otimes \left( {{\mathbf{v}}_{{{\text{IN}}}}^{{\text{R}}} - {\mathbf{V}}} \right) + \left( {{\mathbf{v}}_{{{\text{IN}}}}^{{\text{R}}} - {\mathbf{V}}} \right) \otimes {\mathbf{m}}} \right)} \right\}. \\ \end{gathered} $

4.3. Условие свободной границы

Данное граничное условие записывается следующим образом:

${\mathbf{\sigma }} \cdot {\mathbf{m}} = 0.$

Здесь и далее m – внешняя нормаль к границе (когда речь идет о граничных условиях).

Для реализации данного граничного условия первый этап осуществляется аналогично описанному в п. 4.2. На третьем этапе корректируются значения в граничных точках по формулам

${{{\mathbf{v}}}_{{n + 1}}} = {{{\mathbf{v}}}_{{{\text{IN}}}}} - \frac{1}{{{{\rho }}{{{\text{c}}}_{{\text{S}}}}}}\left( {{{{\mathbf{\sigma }}}_{{{\text{IN}}}}} \cdot {\mathbf{m}}} \right) + \frac{1}{{{\rho }}}\left( {\frac{1}{{{{c}_{{\text{S}}}}}} - \frac{1}{{{{c}_{{\text{P}}}}}}} \right)\left( {{\mathbf{m}} \cdot {{{\mathbf{\sigma }}}_{{{\text{IN}}}}} \cdot {\mathbf{m}}} \right){\mathbf{m}},$
${{{\mathbf{\sigma }}}_{{n + 1}}} = {{{\mathbf{\sigma }}}_{{{\text{IN}}}}} - 2{{{\mathbf{\sigma }}}_{{{\text{IN}}}}} \cdot \left( {{\mathbf{m}} \otimes {\mathbf{m}}} \right) - \frac{{{\mathbf{m}} \cdot {{{\mathbf{\sigma }}}_{{{\text{IN}}}}} \cdot {\mathbf{m}}}}{{c_{{\text{P}}}^{2}}}\left( {c_{3}^{2}{\mathbf{I}} - 2\left( {c_{{\text{P}}}^{2} - c_{{\text{S}}}^{2}} \right)\left( {{\mathbf{m}} \otimes {\mathbf{m}}} \right)} \right).$

Здесь и далее для двух симметричных тензоров ${\mathbf{A}}$ и ${\mathbf{B}}$ второго ранга выполняется условие

${{\left( {{\mathbf{A}} \cdot {\mathbf{B}}} \right)}_{{ij}}} = \sum\limits_k {{{A}_{{ik}}}{{B}_{{kj}}}} = {{\left( {{\mathbf{B}} \cdot {\mathbf{A}}} \right)}_{{ij}}}.$

4.4. Граничное условие нулевого давления

Данное граничное условие записывается следующим образом:

$p = 0.$

Для реализации данного граничного условия первый этап осуществляется аналогично описанному в п. 4.2. На третьем этапе корректируются значения в граничных точках по формулам

${{{\mathbf{v}}}_{{n + 1}}} = {{{\mathbf{v}}}_{{{\text{IN}}}}} + \frac{{{{p}_{{{\text{IN}}}}}}}{{c{{\rho }}}}{\mathbf{m}},$
${{p}_{{n + 1}}} = 0.$

4.5. Контактное условие свободного скольжения

Данное контактное условие записывается следующим образом:

${{{\mathbf{v}}}^{{\text{L}}}} \cdot {\mathbf{m}} = {{{\mathbf{v}}}^{{\text{R}}}} \cdot {\mathbf{m}},$
${{{\mathbf{\sigma }}}^{{\text{L}}}} \cdot {\mathbf{m}} = {{{\mathbf{\sigma }}}^{{\text{R}}}} \cdot {\mathbf{m}},$
${{{\mathbf{\sigma }}}^{{\text{L}}}} \cdot {\mathbf{m}} = \left( {{\mathbf{m}} \cdot {{{\mathbf{\sigma }}}^{{\text{L}}}} \cdot {\mathbf{m}}} \right){\mathbf{m}},$
${{{\mathbf{\sigma }}}^{{\text{R}}}} \cdot {\mathbf{m}} = \left( {{\mathbf{m}} \cdot {{{\mathbf{\sigma }}}^{{\text{R}}}} \cdot {\mathbf{m}}} \right){\mathbf{m}}.$

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

${\mathbf{f}} = \frac{{{{{{\rho }}}^{{\text{L}}}}c_{{\text{P}}}^{{\text{L}}}{{{{\rho }}}^{{\text{R}}}}c_{{\text{P}}}^{{\text{R}}}\left( {{\mathbf{v}}_{{{\text{IN}}}}^{{\text{R}}} - {\mathbf{v}}_{{{\text{IN}}}}^{{\text{L}}}} \right) \cdot {\mathbf{m}} + {{{{\rho }}}^{{\text{L}}}}c_{{\text{P}}}^{{\text{L}}}\left( {{\mathbf{m}} \cdot {\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{R}}} \cdot {\mathbf{m}}} \right) + {{{{\rho }}}^{{\text{R}}}}c_{{\text{P}}}^{{\text{R}}}\left( {{\mathbf{m}} \cdot {\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{L}}} \cdot {\mathbf{m}}} \right)}}{{{{{{\rho }}}^{{\text{L}}}}c_{{\text{P}}}^{{\text{L}}} + {{{{\rho }}}^{{\text{R}}}}c_{{\text{P}}}^{{\text{R}}}}}{\mathbf{m}},$
${{{\mathbf{z}}}^{{\text{L}}}} = {\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{L}}} \cdot {\mathbf{m}} - {\mathbf{f}},$
${{{\mathbf{z}}}^{{\text{R}}}} = {\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{R}}} \cdot {\mathbf{m}} - {\mathbf{f}},$
${\mathbf{v}}_{{n + 1}}^{{\text{L}}} = {\mathbf{v}}_{{{\text{IN}}}}^{{\text{L}}} - \frac{1}{{{{{{\rho }}}^{{\text{L}}}}c_{{\text{S}}}^{{\text{L}}}}}{{{\mathbf{z}}}^{{\text{L}}}} + \frac{1}{{{{{{\rho }}}^{{\text{L}}}}}}\left( {\frac{1}{{c_{{\text{S}}}^{{\text{L}}}}} - \frac{1}{{c_{{\text{P}}}^{{\text{L}}}}}} \right)\left( {{{{\mathbf{z}}}^{{\text{L}}}} \cdot {\mathbf{m}}} \right){\mathbf{m}},$
${\mathbf{v}}_{{n + 1}}^{{\text{R}}} = {\mathbf{v}}_{{{\text{IN}}}}^{{\text{R}}} + \frac{1}{{{{{{\rho }}}^{{\text{R}}}}c_{{\text{S}}}^{{\text{R}}}}}{{{\mathbf{z}}}^{{\text{R}}}} - \frac{1}{{{{{{\rho }}}^{{\text{R}}}}}}\left( {\frac{1}{{c_{{\text{S}}}^{{\text{R}}}}} - \frac{1}{{c_{{\text{P}}}^{{\text{R}}}}}} \right)\left( {{{{\mathbf{z}}}^{{\text{R}}}} \cdot {\mathbf{m}}} \right){\mathbf{m}},$
${\mathbf{\sigma }}_{{n + 1}}^{{\text{L}}} = {\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{L}}} - \left( {{{{\mathbf{z}}}^{{\text{L}}}} \otimes {\mathbf{m}} + {\mathbf{m}} \otimes {{{\mathbf{z}}}^{{\text{L}}}}} \right) + \frac{{{{{\mathbf{z}}}^{{\text{L}}}} \cdot {\mathbf{m}}}}{{{{{\left( {с_{{\text{P}}}^{{\text{L}}}} \right)}}^{2}}}}\left\{ {2\left( {{{{\left( {с_{{\text{P}}}^{{\text{L}}}} \right)}}^{2}} - {{{\left( {с_{{\text{S}}}^{{\text{L}}}} \right)}}^{2}}} \right)\left( {{\mathbf{m}} \otimes {\mathbf{m}}} \right) - {{{\left( {с_{{\text{3}}}^{{\text{L}}}} \right)}}^{2}}{\mathbf{I}}} \right\},$
${\mathbf{\sigma }}_{{n + 1}}^{{\text{R}}} = {\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{R}}} - \left( {{{{\mathbf{z}}}^{{\text{R}}}} \otimes {\mathbf{m}} + {\mathbf{m}} \otimes {{{\mathbf{z}}}^{{\text{R}}}}} \right) + \frac{{{{{\mathbf{z}}}^{{\text{R}}}} \cdot {\mathbf{m}}}}{{{{{\left( {с_{{\text{P}}}^{{\text{R}}}} \right)}}^{2}}}}\left\{ {2\left( {{{{\left( {с_{{\text{P}}}^{{\text{R}}}} \right)}}^{2}} - {{{\left( {с_{{\text{S}}}^{{\text{R}}}} \right)}}^{2}}} \right)\left( {{\mathbf{m}} \otimes {\mathbf{m}}} \right) - {{{\left( {с_{{\text{3}}}^{{\text{R}}}} \right)}}^{2}}{\mathbf{I}}} \right\}.$

4.6. Контактное условие между упругой и акустической средами

Данное контактное условие записывается следующим образом:

(4.10)
${\mathbf{\sigma }} \cdot {\mathbf{m}} + p{\mathbf{m}} = 0,$
(4.11)
${{{\mathbf{v}}}^{{\text{A}}}} \cdot {\mathbf{m}} = {{{\mathbf{v}}}^{{\text{E}}}} \cdot {\mathbf{m}}.$

В (4.10), (4.11) m – внешняя нормаль к упругой среде, здесь и в (4.11) верхние индексы A и E соответствуют акустической и упругой среде соответственно.

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

$P = \frac{{{{c}_{{\text{P}}}}{{{{\rho }}}^{{\text{E}}}}c{{{{\rho }}}^{{\text{A}}}}\left( {{\mathbf{v}}_{{{\text{IN}}}}^{{\text{E}}} - {\mathbf{v}}_{{{\text{IN}}}}^{{\text{A}}}} \right) \cdot {\mathbf{m}} + {{c}_{{\text{P}}}}{{{{\rho }}}^{{\text{E}}}}{{p}_{{{\text{IN}}}}} - c{{{{\rho }}}^{{\text{A}}}}\left( {{\mathbf{m}} \cdot {{{\mathbf{\sigma }}}_{{{\text{IN}}}}} \cdot {\mathbf{m}}} \right)}}{{{{c}_{{\text{P}}}}{{{{\rho }}}^{{\text{E}}}} + c{{{{\rho }}}^{{\text{A}}}}}},$
${\mathbf{v}}_{{n + 1}}^{{\text{A}}} = {\mathbf{v}}_{{{\text{IN}}}}^{{\text{A}}} + \frac{{P - {{p}_{{{\text{IN}}}}}}}{{c{{{{\rho }}}^{{\text{A}}}}}}{\mathbf{m}},$
${{p}_{{n + 1}}} = P,$
${\mathbf{z}} = {{{\mathbf{\sigma }}}_{{{\text{IN}}}}} \cdot {\mathbf{m}} + P{\mathbf{m}},$
${\mathbf{v}}_{{n + 1}}^{{\text{E}}} = {\mathbf{v}}_{{{\text{IN}}}}^{{\text{E}}} - \frac{1}{{{{{{\rho }}}^{{\text{E}}}}{{c}_{{\text{S}}}}}}{\mathbf{z}} + \frac{1}{{{{{{\rho }}}^{{\text{E}}}}}}\left( {\frac{1}{{{{c}_{{\text{S}}}}}} - \frac{1}{{{{c}_{{\text{P}}}}}}} \right)\left( {{\mathbf{z}} \cdot {\mathbf{m}}} \right){\mathbf{m}},$
${{{\mathbf{\sigma }}}_{{n + 1}}} = {{{\mathbf{\sigma }}}_{{{\text{IN}}}}} - \left( {{\mathbf{z}} \otimes {\mathbf{m}} + {\mathbf{m}} \otimes {\mathbf{z}}} \right) + \left( {{\mathbf{z}} \cdot {\mathbf{m}}} \right)\left( {\frac{{{{c}_{3}}}}{{{{c}_{{\text{P}}}}}}{\mathbf{I}} - \left( {1 + \frac{{{{c}_{3}}}}{{{{c}_{{\text{P}}}}}}} \right)\left( {{\mathbf{m}} \otimes {\mathbf{m}}} \right)} \right).$

4.7. Контактное условие между иерархическими вложенными сетками с кратным шагом по координате

Для реализации данного контактного условия значения в дополнительных узлах расчетной сетки с меньшим шагом по координате заполняются значениями из узлов расчетной сетки с большим шагом по координате, либо путем интерполяции по узлам расчетной сетки с большим шагом по координате (см. [27]), в соответствии с фиг. 6а, а действие граничного корректора отсутствует.

Фиг. 6.

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

4.8. Контактное условие между подобластями с различными расчетными схемами с кратным шагом по координате

Для реализации данного контактного условия значения на первом этапе в дополнительных узлах расчетной сетки с меньшим шагом по координате заполняются значениями из узлов расчетной сетки с большим шагом по координате, либо путем интерполяции по узлам расчетной сетки с большим шагом по координате. А значения в дополнительных узлах расчетной сетки с большим шагом по координате заполняются значениями из узлов расчетной сетки с меньшим шагом по координате (фиг. 6б). Обоснование использования резко меняющегося шага по координате дано в [32].

На втором этапе вычисления в расчетной сетке с большим шагом по координате выполняются с помощью диссипативной разностной схемы, например, схемы КИР (см. [34]). А вычисления в расчетной сетке с меньшим шагом по координате выполняются с помощью разностной схемы Русанова (см. [35]).

Так как интерполяция внесет ошибку в вычисления, а шаги в контактирующих сетках могут отличаться в несколько раз, то дополнительно требуется выполнение условия полного слипания на границе. При этом в случае упругого волнового уравнения (2.3), (2.4) будут выполнены выражения (4.4), (4.5). А в случае акустического волнового уравнения (2.1), (2.2) данное контактное условие полного слипания примет следующий вид:

${{{\mathbf{v}}}^{{\text{L}}}} \cdot {\mathbf{m}} = {{{\mathbf{v}}}^{{\text{R}}}} \cdot {\mathbf{m}},$
${{p}^{{\text{L}}}} = {{p}^{{\text{R}}}}.$

Можно показать, что выполнение данных условий полного слипания эквивалентно математически непрерывному счету с точностью до порядка аппроксимации используемых разностных схем (см. [17]).

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

$P = \frac{{p_{{{\text{IN}}}}^{{\text{R}}} + p_{{{\text{IN}}}}^{{\text{L}}} + c{{\rho }}\left( {\left( {{\mathbf{v}}_{{{\text{IN}}}}^{{\text{L}}} - {\mathbf{v}}_{{{\text{IN}}}}^{{\text{R}}}} \right) \cdot {\mathbf{m}}} \right)}}{2},$
${\mathbf{v}}_{{n + 1}}^{{\text{L}}} = {\mathbf{v}}_{{{\text{IN}}}}^{{\text{L}}} + \frac{{p_{{{\text{IN}}}}^{{\text{L}}} - P}}{{c{{\rho }}}}{\mathbf{m}},$
${\mathbf{v}}_{{n + 1}}^{{\text{R}}} = {\mathbf{v}}_{{{\text{IN}}}}^{{\text{R}}} + \frac{{P - p_{{{\text{IN}}}}^{{\text{R}}}}}{{c{{\rho }}}}{\mathbf{m}},$
$p_{{n + 1}}^{{\text{L}}} = p_{{n + 1}}^{{\text{R}}} = P.$

А в случае упругого волнового уравнения (2.3), (2.4) сначала вычисляется величина по формуле

(4.12)
${\mathbf{V}} = \frac{1}{2}\left( {{\mathbf{v}}_{{{\text{IN}}}}^{{\text{L}}} + {\mathbf{v}}_{{{\text{IN}}}}^{{\text{R}}} + \frac{{\left( {{\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{R}}} - {\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{L}}}} \right) \cdot {\mathbf{m}}}}{{{{\rho }}{{c}_{{\text{S}}}}}} - \frac{{{{c}_{{\text{P}}}} - {{c}_{{\text{S}}}}}}{{{{\rho }}{{c}_{{\text{P}}}}{{c}_{{\text{S}}}}}}\left( {{\mathbf{m}} \cdot \left( {{\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{R}}} - {\mathbf{\sigma }}_{{{\text{IN}}}}^{{\text{L}}}} \right) \cdot {\mathbf{m}}} \right){\mathbf{m}}} \right),$

затем ее подставляют в выражения (4.7)–(4.9) из п. 4.2 с учетом того, что

(4.13)
$c_{{\text{P}}}^{{\text{L}}} = c_{{\text{P}}}^{{\text{R}}} = {{c}_{{\text{P}}}},$
(4.14)
$c_{{\text{S}}}^{{\text{L}}} = c_{{\text{S}}}^{{\text{R}}} = {{c}_{{\text{S}}}},$
(4.15)
$c_{3}^{{\text{L}}} = c_{3}^{{\text{R}}} = {{c}_{3}},$
(4.16)
${{{{\rho }}}^{{\text{L}}}} = {{{{\rho }}}^{{\text{R}}}} = {{\rho }}{\text{.}}$

Заметим, что выражение (4.6) переходит в выражение (4.12), если выполняется (4.13)–(4.16).

5. РЕЗУЛЬТАТЫ ТЕСТОВЫХ РАСЧЕТОВ

В данном разделе приведены результаты двух разных расчетов, параметры которых представлены в табл. 3. Отметим, что в табл. 2 из разд. 2 приведены параметры всех отдельных расчетных сеток для расчета № 2.

Таблица 3.  

Параметры расчетов

№ Расстояние от эпицентра до ледового острова, км Ширина центральных расчетных сеток 11–15, км Общее количество точек
во всех расчетных сетках
1 10 35 4 524 970
2 35 55 5 724 970

Ниже на фиг. 7–12 представлены волновые картины в различные моменты времени в различных масштабах по координатам: крупном, среднем и мелком. Стрелками отмечено положение ледового острова. На фиг. 13 представлены финальные разрушения.

Фиг. 7.

Распространение сейсмических волн в подобластях с различными расчетными схемами с кратным шагом по координате, волновые картины модуля скорости в момент времени 6 с: (а) – расчет № 1 с шириной центральных расчетных сеток 35 км, мелкий масштаб; (б) – расчет № 1 с шириной центральных расчетных сеток 35 км, средний масштаб; (в) – расчет № 2 с шириной центральных расчетных сеток 55 км, мелкий масштаб; (г) – расчет № 2 с шириной центральных расчетных сеток 55 км, средний масштаб.

Фиг. 8.

Распространение сейсмических волн в подобластях с различными расчетными схемами с кратным шагом по координате, волновые картины модуля скорости в момент времени 18 с: (а) – расчет № 1 с шириной центральных расчетных сеток 35 км, мелкий масштаб; (б) – расчет № 1 с шириной центральных расчетных сеток 35 км, средний масштаб; (в) – расчет № 2 с шириной центральных расчетных сеток 55 км, мелкий масштаб; (г) – расчет № 2 с шириной центральных расчетных сеток 55 км, средний масштаб.

Фиг. 9.

Воздействие P-волны на ледовый остров, расчет № 1 с расстоянием от эпицентра до ледового острова 10 км, волновая картина модуля скорости в момент времени 6.2 с: (а) – средний масштаб, (б) – крупный масштаб.

Фиг. 10.

Воздействие S-волны на ледовый остров, расчет № 1 с расстоянием от эпицентра до ледового острова 10 км, волновая картина модуля скорости в момент времени 5.4 с: (а) – средний масштаб, (б) – крупный масштаб.

Фиг. 11.

Воздействие P-волны на ледовый остров, расчет № 2 с расстоянием от эпицентра до ледового острова 35 км, волновая картина модуля скорости в момент времени 6.4 с: (а) – средний масштаб, (б) – крупный масштаб.

Фиг. 12.

Воздействие S-волны на ледовый остров, расчет № 2 с расстоянием от эпицентра до ледового острова 35 км, волновая картина модуля скорости в момент времени 7.6 с: (а) – средний масштаб, (б) – крупный масштаб.

Фиг. 13.

Финальные разрушения: (а) – расчет № 1 с расстоянием от эпицентра до ледового острова 10 км, (б) – расчет № 2 с расстоянием от эпицентра до ледового острова 35 км.

На фиг. 10, 12 видно, что за счет условия скольжения между ледовым островом и осадочными породами поперечные волны (по сравнению с продольными, фиг. 9, 11) отражаются аналогично с границей раздела между осадочными породами и водой. Для определения типов волн использовался метод Wave Logica (см. [36]–[40]).

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

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

Также проведенные расчеты демонстрируют, что использование системы вложенных иерархических сеток с кратным шагом по координате позволило рассчитать распространение сейсмических волн от гипоцентра, расположенного на глубине 14 км до ледового острова высотой 10 м и диаметром 300 м, и определить локализацию разрушений.

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

  1. Crawford A.J., Mueller D.R., Humphreys E.R., Carrieres T., Tran H. Surface ablation model evaluation on a drifting ice island in the Canadian Arctic // Cold Regions Sci. and Technology. 2015. V. 110. P. 170–182.

  2. Crawford A., Crocker G., Mueller D., Desjardins L., Saper R., Carrieres T. The canadian ice island drift, deterioration and detection (CI2D3) database // J. Glaciology. 2018. V. 64. № 245. P. 517–521.

  3. Xunqiang Y., Jianbo L., Chenglin W., Gao L. ANSYS implementation of damping solvent stepwise extraction method for nonlinear seismic analysis of large 3-D structures // Soil Dynam. and Earthquake Engineer. 2013. V. 44. P. 139–152.

  4. Nikolic Z., Zivaljic N., Smoljanovic H., Balic I. Numerical modelling of reinforced concrete structures under seismic loading based on the finite element method with discrete inter element cracks // Earthquake Engineer and Structur. Dynam. 2017. V. 46. № 1. P. 159–178.

  5. Yaghin M.L., Hesari M.A. Dynamic analysis of the arch concrete dam under earthquake force with ABAQUS // J. Appl. Sci. 2008. V. 8. № 15. P. 2648–2658.

  6. Moczo P., Robertsson J.O., Eisner L. The finite-difference time-domain method for modeling of seismic wave propagation // Adv. Geophys. 2007. V. 48. P. 421–516.

  7. Saenger E.H., Shapiro S.A. Effective velocities in fractured media: a numerical study using the rotated staggered finite-difference grid // Geophys. Prospecting. 2002. V. 50. № 2. P. 183–194.

  8. Slawinski R.A., Krebes E.S. Finite-difference modeling of SH-wave propagation in nonwelded contact media // Geophysics. 2002. V. 67. № 5. P. 1656–1663.

  9. Petrov I.B., Favorskaya A.V., Favorskaya M.N., Simakov S.S., Jain L.C. Development and applications of computational methods // Smart Innovation, Systems and Technologies. 2019. V. 133. P. 3–7.

  10. Komatitsch D., Vilotte J.P., Vai R., Castillo-Covarrubias J.M., Sanchez-Sesma F.J. The spectral element method for elastic wave equations-application to 2-D and 3-D seismic problems // Internat. J. for Numer. Meth. Engineer. 1999. V. 45. № 9. P. 1139–1164.

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

  12. Priolo E., Carcione J.M., Seriani G. Numerical simulation of interface waves by high-order spectral modeling techniques // J. of the Acoustical Soc. Am. 1994. V. 95. № 2. P. 681–693.

  13. Kaser M., Dumbser M. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes – I. The two-dimensional isotropic case with external source terms // Geophys. J. Intern. 2006. V. 166. № 2. P. 855–877.

  14. Dumbser M., Kaser M. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes – II. The three-dimensional isotropic case // Geophys. J. Intern. 2006. V. 167. № 1. P. 319–336.

  15. Wilcox L.C., Stadler G., Burstedde C., Ghattas O. A high-order discontinuous Galerkin method for wave propagation through coupled elastic–acoustic media // J. Comput. Phys. 2010. V. 229. № 24. P. 9373–9396.

  16. De Basabe J., Mrinal S., Wheeler M. The interior penalty discontinuous Galerkin method for elastic wave propagation: grid dispersion // Geophys. J. Intern. 2008. V. 175. № 1. P. 83–93.

  17. Favorskaya A.V., Khokhlov N.I., Petrov I.B. Grid-characteristic method on joint structured regular and curved grids for modeling coupled elastic and acoustic wave phenomena in objects of complex shape // Lobachevskii J. of Math. 2020. V. 41. № 4. P. 512–525.

  18. 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.

  19. Favorskaya A.V., Breus A.V., Galitskii B.V. Application of the grid-characteristic method to the seismic isolation model // Smart Innovation, Systems and Technologies. 2018. V. 133. P. 167–181.

  20. Favorskaya A., Golubev V., Grigorievyh D. Explanation the difference in destructed areas simulated using various failure criteria by the wave dynamics analysis // Procedia Computer Sci. 2018. V. 126. P. 1091–1099.

  21. Karulina M., Marchenko A., Karulin E., Sodhi D., Sakharov A., Chistyakov P. Full-scale flexural strength of sea ice and freshwater ice in Spitsbergen Fjords and North-West Barents Sea // Appl. Ocean Res. 2019. V. 90. Paper № 101853.

  22. Schwarz J., Frederking R., Gavrilo V., Petrov I.G., Hirayama K.I., Mellor M., Tryde P., Vaudry K.D. Standardized testing methods for measuring mechanical properties of sea ice // Cold Regions Sci. and Technology 1981. V. 4. № 3. P. 245–253.

  23. Marchenko A., Karulin E., Chistyakov P., Sodhi S., Karulina M., Sakharov A. Three dimensional fracture effects in tests with cantilever and fixed ends beams // Proc. of the 22th IAHR Symp. on Ice. 2014. Paper № 1178.

  24. Jones S.J., Gagnon R.E., Derradji A., Bugden A. Compressive strength of iceberg ice // Canadian J. of Phys. 2003. V. 81. № 1–2. P. 191–200.

  25. Jones S.J. A review of the strength of iceberg and other freshwater ice and the effect of temperature // Cold Regions Sci. and Technology. 2007. V. 47. № 3. P. 256–262.

  26. Favorskaya A.V., Petrov I.B. A study of high-order grid-characteristic methods on unstructured grids // Numer. Anal. and Appl. 2016. V. 9. № 2. P. 171–178.

  27. Petrov I.B., Favorskaya A.V., Khokhlov N.I. Grid-characteristic method on embedded hierarchical grids and its application in the study of seismic waves // Comput. Math. and Math. Phys. 2017. V. 57. № 11. P. 1771–1777.

  28. Breus A., Favorskaya A., Golubev V., Kozhemyachenko A., Petrov I. Investigation of seismic stability of high-rising buildings using grid-characteristic method // Procedia Computer Sci. 2019. V. 154. P. 305–310.

  29. Favorskaya A., Golubev V. Study the elastic waves propagation in multistory buildings, taking into account dynamic destruction // Smart Innovation, Systems and Technologies. 2020. V. 193. P. 189–199.

  30. Favoskaya A.V., Petrov I.B. Calculation the earthquake stability of various structures using the grid-characteristic method // Radioelektronika, Nanosistemy, Informacionnye Tehnologii. 2019. V. 11. № 3. P. 345–350.

  31. Favorskaya A. Computation the bridges earthquake resistance by the grid-characteristic method // Smart Innovation, Systems and Technologies. 2020. V. 193. P. 179–187.

  32. 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. Prospecting. 2018. V. 66. № 8. P. 1485–1502.

  33. Favorskaya A.V., Petrov I.B. Grid-characteristic method // Smart Innovation, Systems and Technologies. 2018. V. 90. P. 117–160.

  34. Courant R., Isaacson E., Rees M. On the solution of nonlinear hyperbolic differential equations by finite differences // Comm. on Pure and Appl. Math. 1952. V. 5. № 3. P. 243–255.

  35. Kholodov A.S., Kholodov Ya.A. Monotonicity criteria for difference schemes designed for hyperbolic equations // Comput. Math. and Math. Phys. 2006. V. 46. № 9. P. 1560–1588.

  36. Favorskaya A., Petrov I. A novel method for investigation of acoustic and elastic wave phenomena using numerical experiments // Theoret. and Appl. Mech. Lett. 2020. V. 10. № 5. P. 307–314.

  37. Favorskaya A.A. Novel method for wave phenomena investigation // Procedia Computer Sci. 2019. V. 159. P. 1208–1215.

  38. Favorskaya A., Kabanova A., Petrov I. Applying Wave Logica method for geophysical prospecting // Geomodel. 2018. P. 1–5.

  39. Favorskaya A.V. Elastic wave scattering on a gas-filled fracture perpendicular to plane P-wave front // Smart Innovation, Systems and Technologies. 2020. V. 173. P. 213–224.

  40. Favorskaya A.V., Petrov I.B. The use of full-wave numerical simulation for the investigation of fractured zones // Math. Models and Computer Simulations. 2019. V. 11. № 4. P. 518–530.

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

Инструменты

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

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

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