МЕТОД РАСЧЕТА КИНЕМАТИКИ ВОЛНОВОГО ФРОНТА ЦУНАМИ В СЕТОЧНОЙ ОБЛАСТИ Текст научной статьи по специальности «Физика»

Научная статья

УДК 550.344.42

DOI 10.25205/1818-7900-2022-20-1-57-66

Метод расчета кинематики волнового фронта цунами

в сеточной области

Андрей Гурьевич Марчук

Институт вычислительной математики и математической геофизики Сибирского отделения Российской академии наук Новосибирск, Россия

Новосибирский государственный университет Новосибирск, Россия

mag@omzg.sscc.ru

Аннотация

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

волна цунами, кинематика волнового фронта, цифровая сеточная батиметрия, времена вступления цунами Благодарности

Исследование выполнено в рамках Государственной бюджетной программы с ИВММГ СО РАН (№ 03152019-0005) Для цитирования

Марчук Ан. Г. Метод расчета кинематики волнового фронта цунами в сеточной области // Вестник НГУ. Серия: Информационные технологии. 2022. Т. 20, № 1. С. 57-66. БО! 10.25205/1818-7900-2022-20-1-57-66

Method for Tsunami Wave Kinematics Modeling in a Gridded Area

Andrey G. Marchuk

Institute of Computational Mathematics and Mathematical Geophysics of the Siberian Branch of the Russian Academy of Sciences Novosibirsk, Russian Federation

Novosibirsk State University Novosibirsk, Russian Federation

mag@omzg.sscc.ru

Abstract

The new method for the tsunami travel times calculation to the nodes of rectangular grid was developed and tested. The algorithm is based on constructing of the wave front segment inside the grid cell using tsunami arrival times to the cell corners. Testing the method proposed shows the better quality of the wave-front line approximation as compared to the existing method developed by the author earlier.

© Марчук Ан. Г., 2022

Keywords

tsunami wave, kinematics of the wave front, digital gridded bathymetry, tsunami travel times Acknowledgements

This work was carried out under state contract with ICMMG SB RAS (no. 0315-2019-0005) For citation

Marchuk An. G. Method for Tsunami Wave Kinematics Modeling in a Gridded Area. Vestnik NSU. Series: Information Technologies, 2022, vol. 20, no. 1, pp. 57-66. (in Russ.) DOI 10.25205/1818-7900-2022-20-1-57-66

Введение

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

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

1. Основы кинематики цунами

Волна цунами относится к классу длинных волн, для описания которых используется система дифференциальных уравнений мелкой воды [7]. В случае отсутствия внешних сил, кроме гравитации, эти уравнения имеют вид

Н + (иН )х + (уН) у = 0, и + иих + уиу + %НХ = фх, (1)

V + иух + + 8Ну = , Здесь H = D + п представляет полную толщину водного слоя, п - вертикальное смещение водной поверхности относительно невозмущенного уровня, u и v - компоненты горизонтальной скорости водного потока, D - глубина, g - ускорение силы тяжести. Из уравнений (1) следует, что скорость волнового фронта длинной волны не зависит от ее параметров (длины и высоты) и выражается формулой [7]

c =

JgD. (2)

Из той же системы уравнений (1) следует, что точки волнового фронта продвигаются с указанной в формуле Лагранжа (2) скоростью в направлении нормали к линии этого волнового фронта. Здесь следует заметить, что из-за амплитудной нелинейности скорость гребня волны несколько отличается от скорости распространения фронта и составляет

С = + , (3)

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

рыми равно Ь, а глубина линейно меняется от величины Нх до величины Н2. Введем вспомогательную величину - угол наклона дна tg(a) = |(Н2 — Нх)|/ Ь . Тогда время пробега вдоль отрезка длиной Ь запишется в виде

L dl i LL h Y1/2 /, H л

l н--1

Г = = fil н

Wg • (H +1 • tga) yjg • tga i V tgaj

2 Г, H^

л/g • tga

1/2

l + Hl V tgaJ

d

V tga j

2 VH2 — 4H1

л/g • tga yftgâ

(4)

2 H2 — Hx 2L

^■ tga ^Н + л/Н

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

2. Описание метода

Фронт волны цунами продвигается в направлении внешней нормали к линии этого фронта, которую можно определить как границу раздела между точками акватории, куда возмущение к этому моменту уже пришло, и всеми остальными точками области [8]. В публикации [9] предложен численный метод ортогонального продвижения волнового фронта для расчета кинематики цунами. Но для его реализации требуются значения глубины в любой точке области, которые вычисляются билинейным интерполированием сеточной батиметрии.

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

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

L

0

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

Fig. 1. The location of the rectangular grid nodes representing the initial circular wave front. The length of the spatial steps along the abscissa and the ordinate axes are Ax and Ay, respectively

(i-lj)

(i-lj-1)

(ij-l)

(i+lj-1) (MJ-1)

(id)

(U-l)

(i+lj-1)

б

в

Рис. 2. Шаблоны для расчета времени прихода в точку (i, j): а - шаблон № 1; б - шаблон № 2; в - шаблон № 3 Fig. 2. The stencils for the computation of the wave arrival time to the grid node (i, j) а - stencil no. 1; b - stencil no. 2; c - stencil no. 3

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

иным образом. В шаблоне № 1 (рис. 2, а) узел сетки с индексами (1 + 1,/ - 1) может отсутствовать. Далее опишем алгоритм расчета времени прихода цунами в узел расчетной сетки при конфигурации, приведенной на рис. 2, а. Рисунок 5 показывает схему этого расчета, основанную на простых геометрических построениях.

Рис. 3. Схема вычисления времени прихода в точку A на основе известных времен прихода волны в точки B, C и D. Во всех углах ячейки сетки глубины заданы

Fig. 3. The scheme for calculating the arrival time at the point A based on the known arrival times of the wave at the points B, C and D. At the angles of the grid cell depths are given

Рассмотрим ячейку расчетной сетки, где требуется определить время вступления волны ТА в точку А при известных временах прихода цунами Тв, ТС и Т0 в остальные три угла ячейки В, С и D (рис. 2, в). Значения глубины НА, Нв, НС и Н0 известны во всех углах ячейки. Пусть для определенности Тв > Т0. В этом случае точка пересечения отрезка волнового фронта, соответствующего моменту времени Т0, и отрезка ВС (вместе с этим и угол а) находится с использованием формулы (4). Принимая во внимание то, что скорость движения точки пересечения фронта с отрезком ВС больше скорости волнового фронта в 1/ео8(а) раз, величина ео8(а) вычисляется по формуле

С05(й)= (5)

2Ау

Далее, легко найти расстояние I от точки А до отрезка волнового фронта DE:

I = Ау ■ сов(а) (6)

Теперь для вычисления времени вступления цунами в точку А достаточно определить глубину в точке F, определяемую ее координатами. Эти координаты ХР и УР привязаны к координатам Х0 и У0 точки D по формулам

Хр. = — /2 = — I ■ зт(а) , (7)

у, = + ^ = Уд +1 ■ соэ(а). (8)

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

Н = НС (ХП - ХЕ )+ И и ' ХЕ УБ - ^ , НВ (-ХП - ХЕ )+ НА ' ХЕ ^

Р Лх Ау Лх Ау В результате искомое время вступления волны в точку А выразится в виде

21

ТЛ ТЦ + ' (10)

где l и глубина HF вычисляются по формулам (6) и (9) на основе известных значений глубины во всех углах ячейки и значений времени прихода цунами в три из них (в точки В, С и D). Аналогично с помощью геометрических построений вычисляются времена вступления волны в узлы сетки для других конфигураций узлов с известными и неизвестными значениями времени (рис. 2, б, в).

3. Тестирование метода

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

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

Fig. 4. Sequential positions of the wave front when calculating the kinematics of the oblique front in an area with a constant depth (the zone of initial disturbance where t(i, j) = 0 is colored by green)

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

Н = 0.001-(у +100)2, 1=1,...,1000, j=1,..., 1000. (11)

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

Рис. 5. Положения фронта цунами над параболическим рельефом дна через каждые 200 с Fig. 5. Tsunami wave front positions every 200 sec above the parabolic bottom topography

Из рис. 5 видно, что волновой фронт остается прямолинейным на протяжении всего времени распространения цунами. Это идеально соответствует точному решению для такого модельного распределения глубин в расчетной области [10].

Еще одним тестом для проверки корректности метода является расчет кинематики изначально кругового волнового фронта в области с постоянной глубиной. Теоретически в этом случае линия волнового фронта в любой момент времени будет иметь форму окружности. Для тестирования предложенного метода в центре расчетной области 1000 х 1000 узлов располагается очаг цунами, ограниченный окружностью с радиусом 5000 м. Длина шага сетки составляет 100 м в обоих направлениях, а глубина во всей области постоянна и равна 1000 м. Описанным алгоритмом с учетом трех возможных конфигураций узлов сетки с известными и неизвестными временами вступления туда волны (см. рис. 2) производятся многократные переборы всех расчетных узлов и рассчитываются времена прихода цунами в узлы сетки. Для корректной работы алгоритма направления перебора узлов меняются на противополож-

ные через раз. На рис. 6 визуализированы изолинии поля времен вступления волны в узлы сетки.

Рис. 6. Последовательные положения фронта волны от круглого источника в области с постоянной глубиной (1000 м), найденные описанным методом (a) и последовательные положения волнового фронта от точечного источника, построенные сеточным методом (b), реализованным автором в 1988 г. [8]

Fig. 6. The isolines of the calculated wave arrival times from a round source to the nodes of the calculation grid when using the method proposed (a) and as a result of the calculation (b) using the algorithm developed by author in 1988 [8]

Из рис. 6, a видно, что форма вычисленного волнового фронта в любой момент времени близка к окружности. Это еще раз подтверждает корректность предложенного метода. Для сравнения на рис. 6, Ь результаты тестирования той же задачи с использованием сеточного метода, основанного на принципе Гюйгенса [8], представлены в том же виде. Можно видеть, что с помощью алгоритма, описанного в этой статье, полученные изохроны ближе к точному решению, чем положения волнового фронта, полученные сеточным методом [8]. Количественная разница (в секундах) между полученными временами вступления волны в узлы сетки и точным решением визуализирована на рис. 7 в виде изолиний и цветовой легенды.

Рис. 7. Распределение (в виде изолиний и цветовой палитры) отклонения от точного решения вычисленных времен прихода в узлы регулярной сетки цунами, генерированного источником круглой формы Fig. 7. Distribution (in the form of isolines and a color palette) of the deviation from the exact solution of the calculated travel times at the nodes of the regular grid when tsunami generated by the circular source

Из рис. 7 видно, что в некоторых направлениях разница между вычисленным и точным решением кинематической задачи близка к нулю (< 1 с) вплоть до границ расчетной области. Но имеются 4 азимута, где эта разница постепенно нарастает при удалении от источника, но не превышает 2 % от времени движения фронта волны в этих направлениях. Это можно объяснить неточным вычислением направлений сегментов волнового фронта при некоторых конфигурациях узлов с вычисленными и неизвестными временами вступления туда цунами.

Заключение

Разработан новый метод нахождения времени вступления волны в узлы прямоугольной расчетной сетки, где глубина задана. Он основан на нахождении направления внешней нормали к сегменту волнового фронта и продвижении волнового фронта в этом направлении. Тестирование метода на нескольких тестовых заданиях показало его достоверность и эффективность. Относительные ошибки в вычисленных временах вступления в узлы прямоугольной сетки не превышают 1-2 %.

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

1. Shokin Yu. I., Chubarov L. B., Novikov V. A., Sudakov A. N. Calculations of tsunami travel time charts in the Pacific Ocean - Models, Algorithms, Techniques, Results. Science of Tsunami Hazards, 1987, vol. 5, no. 2, pp. 85-113.

2. Kenji Satake. Effects of bathymetry on tsunami propagation: Application of ray tracing to tsunamis. Pure and Applied Geophysics, 1988, vol. 126, iss. 1, pp. 27-36.

3. Кремлев А. Н. Преломление плоской волны на выпуклом и вогнутом углах в приближении геометрической акустики // Сиб. журн. вычисл. математики. 2017. Т. 20, № 3. С.251-271.

4. Titov V. V., Gonzalez F. Implementation and Testing of the Method of Splitting Tsunami (MOST). Washington DC, 1997. (Technical Memorandum / National Oceanic and Atmospheric Administration; ERL PMEL-112)

5. Imamura F. Tsunami numerical simulation with the staggered leap-frog scheme (numerical code of TUNAMI-N1 and N2). In: Disaster Control Research Center. Tohoku University, 1995, 33 p.

6. Лаврентьев М. М., Лысаков К. Ф., Марчук Ан. Г., Облаухов К. К. Ускорение расчетов распространения волны цунами с использованием FPGA // Успехи кибернетики. 2020. Т. 1, № 4. С. 42-54.

7. Стокер Дж. Дж. Волны на воде. М.: ИЛ, 1959. 617 с.

8. Марчук Ан. Г. Численные методы расчета кинематики волн цунами // Математические проблемы геофизики: Численные исследования геофизических задач. Новосибирск, 1988. С. 69-90.

9. Marchuk An. G., Vasiliev G. S. The fast method for a rough tsunami amplitude estimation //

Bulletin of the Novosibirsk Computing Center. Series: Math. Model. in Geoph., 2014, no. 17, pp.21-34.

10. Марчук Ан. Г. Оценка высоты цунами, распространяющейся над параболическим дном, в лучевом приближении // Сиб. журн. вычисл. математики. 2017. Т. 20, № 1. С. 23-35.

References

1. Shokin Yu. I., Chubarov L. B., Novikov V. A., Sudakov A. N. Calculations of tsunami travel time charts in the Pacific Ocean - Models, Algorithms, Techniques, Results. Science of Tsunami Hazards, 1987, vol. 5, no. 2, pp. 85-113.

2. Kenji Satake. Effects of bathymetry on tsunami propagation: Application of ray tracing to tsunamis. Pure and Applied Geophysics, 1988, vol. 126, iss. 1, pp. 27-36.

3. Kremlev A. N. Prelomlenie ploskoy volny na vypuklom i vognutom uglakh v priblizhenii geometricheskoy akustiki. Sib. zhurn. vychisl. matematiki, 2017, vol. 20, no. 3, pp. 251-271. (in Russ.)

4. Titov V. V., Gonzalez F. Implementation and Testing of the Method of Splitting Tsunami (MOST). Washington DC, 1997. (Technical Memorandum / National Oceanic and Atmospheric Administration; ERL PMEL-112)

5. Imamura F. Tsunami numerical simulation with the staggered leap-frog scheme (numerical code of TUNAMI-N1 and N2). In: Disaster Control Research Center. Tohoku University, 1995,33 p.

6. Lavrentev M. M., Lysakov K. F., Marchuk An. G., Oblaukhov K. K. Uskorenie raschetov rasprostraneniya volny tsunami s ispol'zovaniem FPGA. Uspekhi kibernetiki, 2020, vol. 1, no. 4, pp. 42-54. (in Russ.)

7. Stoker J. J. Volny na vode. Moscow, 1959, 617 p. (in Russ.)

8. Marchuk An. G. Chislennye metody rascheta kinematiki voln tsunami. In: Matematicheskie problemy geofiziki: Chislennye issledovaniya geofizicheskikh zadach. Novosibirsk, 1988, pp. 69-90. (in Russ.)

9. Marchuk An. G., Vasiliev G. S. The fast method for a rough tsunami amplitude estimation //

Bulletin of the Novosibirsk Computing Center. Series: Math. Model. in Geoph., 2014, no. 17, pp.21-34.

10. Marchuk An. G. Otsenka vysoty tsunami, rasprostranyayushcheysya nad parabolicheskim dnom, v luchevom priblizhenii. Sib. zhurn. vychisl. matematiki, 2017, vol. 20, no. 1, pp. 23-35. (in Russ.)

Информация об авторе

Андрей Гурьевич Марчук, доктор физико-математических наук SPIN 4054-6578

WoS Researcher ID S-9502-2017 Scopus Author ID 57023485400

Information about the Author

Andrey G. Marchuk, Doctor of Sciences (Physics and Mathematics) SPIN 4054-6578

WoS Researcher ID S-9502-2017 Scopus Author ID 57023485400

Статья поступила в редакцию 11.12.2021; одобрена после рецензирования 01.02.2022; принята к публикации 01.02.2022 The article was submitted 11.12.2021; approved after reviewing 01.02.2022; accepted for publication 01.02.2022