Параллельная реализация алгоритма для расчета генерации длинных поверхностных волн цунами движением оползня Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 14, № 1, 2009

Параллельная реализация алгоритма

для расчета генерации длинных поверхностных

*

волн цунами движением оползня

Е. Н. Березин Кемеровский государственный университет, Россия e-mail: ben@kemsu.ru

С. А. БЕйЗЕЛЬ

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: beisel_s@ngs.ru

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

Ключевые слова: численное моделирование, методы граничных элементов, идеальная несжимаемая жидкость, свободная граница, потенциальные течения, нестационарные течения.

Введение

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

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

В статьях [3, 4] исследовалась возможность использования приближенных математических моделей гидродинамики для моделирования механизма движения оползня, при

* Работа выполнена при поддержке Российского фонда фундаментальных исследований (гранты № 05-05-64460, 06-05-64869, 06-05-72014, 07-05-13583), Программы междисциплинарных интеграционных исследований СО РАН (проекты 2006-28, 2006-113), Программы интеграционных фундаментальных исследований СО РАН (проекты 28, 113), Программы государственной поддержки научных исследований, проводимых ведущими научными школами Российской Федерации НШ-9886.2006.9, проекта INTAS 06-1000013-9236.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2009.

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

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

Увеличение производительности вычислительных систем ведется в двух основных направлениях — совершенствование элементной базы и развитие параллельных вычислительных алгоритмов. Однако всегда существует предел быстродействия, объясняемый конкретным этапом развития производства, кроме того, передовые технологии всегда более дорогостоящие. И тогда на первый план выступают высокопроизводительные вычислительные алгоритмы.

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

1. Постановка задачи

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

Д^ = 0, х(х,у) е О, (1)

Г1 — свободная граница жидкости, Г2 — поверхность дна, на которой располагается оползень, Г3 — твердые границы бассейна. На границе Г3 ставится условие непротекания

Щ = 0, х(х,у) е Г3. (2)

На границе оползня

дП = ип, х(х,у) е Г2. (3) На свободной границе ставятся кинематическое и динамическое условия

^х = х(х,у) е Г1, (4)

^ - 2 1^12 + у = 0, х(х,у) е Г1. (5)

Рис. 1. Схема расчетной области

Для удобства численной реализации уравнения (1)-(5) приведены к безразмерному виду. В качестве характерных размерных величин были выбраны ускорение свободного падения д и глубина бассейна Н. Все геометрические размеры и обозначения показаны на рис. 1.

Величина й характеризует заглубление тела до начала движения. В начальный момент времени £ = 0 точка уреза совпадает с началом декартовой системы координат хОу. Граница Г2 образует с горизонтом угол в, а ее форма задается функцией Л*(х,£) [4, 8]:

Л* (х, ¿) = —х 1ап в + ЛЛ-

1 + 1апЬ

хх хх с (£)

Б

+ ь

1 — 1апЬ

'X 'X с (¿)

Б

-ь

(6)

[1+1апЬ(Ь)]2

Здесь ЛЛ — максимальная высота оползня; Б — параметр, характеризующий крутизну поверхности оползня; величина Ь равна расстоянию между точками перегиба кривой, описывающей границу оползня; хс(£) — координата точки максимальной толщины тела. В точке х* эта граница соединяется с областью постоянной глубины Н. Выбор такого, достаточно сложного, выражения в правой части (6) обусловлен необходимостью обеспечения гладкости профиля донной поверхности.

2. Алгоритм движения по времени

Краевая задача (1)-(5) является нестационарной, но время явно входит только в граничные условия (4), (5), представляющие собой обыкновенные дифференциальные уравнения первого порядка, для интегрирования которых используется явный метод Эйлера первого порядка точности.

Пусть в некоторый момент времени ¿к заданы положение свободной границы Г и распределение потенциала фк на ней. Необходимо решить уравнение Лапласа (1) с условием фк на Г! и условиями (2), (3) на Г2 и Г3. Новое положение свободной границы и распределение потенциала на ней для момента времени ¿к + т можно вычислить, используя условия (4) и (5), дискретный аналог которых расписывается по схеме Эйлера следующим образом:

хк+1 = хк + (уф)к т,

Фк+1 = Фк + (0, 5 |(Уф)к|2 - ук)т, к = 0,1, 2..., (7)

где хк,фк — значения функций на к-м шаге по времени. В результате, на каждом шаге по времени ¿к решается смешанная краевая задача для уравнения Лапласа (1) с граничными условиями (2), (3) и (5).

Для описания траекторий частиц (точек) свободной границы применяется метод Лагранжа. Приведенная выше схема Эйлера условно устойчива и требует выполнения некоторых дополнительных условий на шаг по времени, для чего авторами используется хорошо зарекомендовавший себя алгоритм автоматического выбора шага по времени, описанный в работе [6]. Контроль за точностью численных расчетов осуществляется проверкой выполнения законов сохранения полной энергии Е = Ек + Ер и массы М, где Ек — кинетическая, Ер — потенциальная энергия.

3. Реализация параллельного метода граничных элементов

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

Ускорением параллельного алгоритма называют отношение времени выполнения алгоритма на одном процессоре ко времени выполнения алгоритма на системе из р процессоров:

Sp = Т1 /Тр. (8)

Эффективностью параллельного алгоритма Ер называют отношение его ускорения к числу процессоров, на котором это ускорение получено:

Е = ^ ■ №

Последовательный алгоритм метода граничных элементов и его распараллеливание

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

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

С (£,х) = -2-Ь (г). (10)

В этом случае справедливо следующее интегральное уравнение:

в (0 Ф (£) + / Ф(х) ^СП^ аг(*) = / Щх) С (е,х) dг(x), (11)

г г

где Ф(х) — гармоническая функция, п(х) — внешняя по отношению к области О единичная нормаль границы Г(х). Параметр ) определяется следующим образом: ) = 2п

для внутренней точки, ) = п для точки на гладкой границе, ) = а для угловой точки границы (а — угол при вершине). Таким образом, на каждом шаге по времени решается интегральное уравнение (11) по границе Г, к которому сводится уравнение Лапласа внутри области D.

Для решения нестационарных задач со свободными границами используется параллельный алгоритм [9].

Ввод данных. Для расчета необходимы следующие граничные и начальные данные:

1) xi, yi,i = 1, п, — узловые координаты точек на границе области;

2) Codei, i = 1, n, — тип граничных условий в i-м узле: если Codei = 0, то в узле задана функция потенциала р, если Codei = 1, то задана функция нормальной производной др/дп;

3) Ui,Bi,i = 1, п, — заданные значения функций р и др/дп соответственно в i-м узле.

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

Вычисление интегралов и решение СЛАУ. Весьма трудоемкими процессами являются вычисление интегралов и решение системы линейных алгебраических уравнений

AX = F, (12)

где A — полностью заполненная матрица размерности п х n, X — вектор неизвестных значений функций р и др/дп на границе Г, F — вектор правой части.

На этом шаге алгоритма система (12) формируется из решения интегрального уравнения (11) по границе Г. Каждая колонка результирующей матрицы A строится независимо от других, следовательно, степень параллелизма этого блока алгоритма максимальна и равна п. Здесь следует позаботиться о циклично-слоистом распределении колонок матрицы A по процессорам для ее дальнейшего решения с помощью метода Гаусса с выбором ведущего элемента. Так как на этом шаге распределение матрицы A по процессорам колоночное, то для LU-факторизации целесообразно использовать колоночно-ориентированный алгоритм [9].

Построение результирующей функции. Если алгоритм решения получившейся в результате LU-факторизации системы уравнений распараллелен, то результат решения распределен по процессорам, если нет — то он собирается на одном процессоре. Обозначим вектор решения через C = (Ci),i = 1, п. Суть построения результирующей функции заключается в следующем: если в i-м узле задан потенциал р, то значение компоненты вектора решения Ci равно результирующей функции Bi; если в i-м узле задана функция нормальной производной др/дп, то значение компоненты вектора решения Ci равно результирующей функции Ui.

Вычисление скоростей. На этом шаге алгоритма по известным узловым координатам xi,yi,i = 1,п, и найденным функциям Ui,Bi вычисляются узловые компоненты вектора скорости vxi ,vyi. Степень параллелизма этого блока равна п. При распараллеливании следует учесть, что массивы xi,yi,i = 1, п, уже имеются на каждом процессоре.

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

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

Для определения шага необходимы следующие величины: £тах — длина максимального граничного элемента свободной границы, £тщ — длина минимального элемента, ^тах — модуль максимальной скорости частиц. Каждый процессор вычисляет временной шаг, исходя из условия [6]:

дт = 2«,

V

где К — заданное максимальное перемещение частиц за один временной шаг, V - модуль максимальной скорости частиц. Величина 7 характеризует некоторую относительную меру дискретизации границы расчетной области.

Перестроение свободной границы. Новое положение свободной границы и распределение потенциала на ней находятся по схеме Эйлера (7). Для этого необходимы следующие данные: , vУi, г = 1 , п. Так как в результате предыдущих шагов ал-

горитма все массивы, кроме vxi, целиком хранятся на каждом процессоре, то распараллеливание этого блока обусловлено хранением vxi, vУi на процессорах. Степень параллелизма этого блока также максимальна и равна п.

Когда распараллеливание этого блока выполнено, для продолжения расчета необходимо разослать на все процессоры следующие данные: Х^у^, Бг,г = 1, п, — новое положение свободной поверхности, распределение потенциала и нормальной производной на ней. Далее необходимо выполнить проверку условия останова вычислений.

4. Тестирование вычислительного алгоритма

Для проверки эффективности и степени ускорения параллельного алгоритма МГЭ при решении двумерных задач идеальной однородной несжимаемой жидкости была проведена серия экспериментов на задаче о трансформации уединенной волны при взаимодействии с вертикальной стенкой. Исследование производительности и эффективности программного кода выполнялось в зависимости от размерности задачи и числа процессоров на кластерах кафедры ЮНЕСКО Кемеровского государственного университета (КемГУ) и Института вычислительных технологий СО РАН (ИВТ СО РАН). Параметры оборудования вычислительного кластера КемГУ

— 8 персональных компьютеров;

— сеть Fast Ethernet. Характеристики каждого компьютера

— процессор Intel Pentium-III с тактовой частотой 667 МГц;

— 64 Мбайт ОЗУ;

— жесткий диск емкостью 20 Гбайт;

— 512 Кбайт КЭШ.

Коммуникационная сеть Fast Ethernet

— пиковая пропускная способность — 100 Мбайт/с;

— коммутатор, полный дуплекс. Программное обеспечение

— ОС Linux на основе дистрибутива Red Hat Linux 8;

— компилятор Intel C/C++, FORTRAN 8.1;

— отладчик Gnu debugger 6.4;

— MPICH 1.2.5.

Параметры оборудования вычислительного кластера ИВТ СО РАН

— 4 вычислительных узла: двухпроцессорный узел на материнской плате Intel SE7501CW2 с процессорами Intel Xeon 3 ГГц, 2 Гбайта ОЗУ и жестким диском емкостью 40 Гбайт;

— i управляющий узел: двухпроцессорный узел на материнской плате Intel SE7501WV2 с процессорами Intel Xeon 3 ГГц, 2 Гбайта ОЗУ и жестким диском емкостью i46 Гбайт на шине SCSI;

— соединение узлов GigE-сетью на основе коммутатора 3Com 3C17401;

— источник бесперебойного питания APC Smart UPS на 3000 В •А.

Программное обеспечение

— управление очередью задач посредством ПО SUN Grid Engine 5.3;

— ОС Linux на основе дистрибутива Red Hat Linux 9;

— компилятор Intel C/C++ версии 8.1;

— отладчик Gnu debugger 6.4;

— MPI MPICH версии 1.2.7;

— библиотека FFTW (Fastest Fourier Transform in the West) версии 2.1.5.

Нестационарное движение уединенной волны по бассейну с накатом на вертикальную стенку

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

Проверка численного алгоритма выполнялась на задаче о взаимодействии уединенной волны амплитуды A = O.l, 0.2, 0.3, 0.4, 0.5, 0.6 с вертикальной стенкой, расположенной в точке x = 35. В качестве характерных безразмерных параметров выбирались глубина канала H = l и ускорение силы тяжести g = l. Шаг по времени подбирался автоматически. Для расчета была взята область D = {-35 < x < 35; — l < y < y0}. В начальный момент времени t = 0.0 вершина волны находилась в точке x = 0. На границе области было задано 540 узлов, из них 301 — на свободной поверхности.

Размещение правой стенки на удалении 35 единиц от начального положения обеспечивало свободный пробег волны протяженностью от одной ее длины (при A = O.l) до трех длин (при A = 0.6). Длина волны l определялась величиной отрезка по оси x, на котором выполняется условие y(t) У 0.0lA(t) [11].

Серия численных экспериментов выполнялась для разного количества узлов N на границе расчетной области (от 400 до 1600 с шагом 300) и разного количества процессоров, задействованных в расчете (2, 4, 8). Программа запускалась 3-5 раз, для каждого запуска считалось среднее время выполнения одной итерации метода, затем выбиралось среднее, наименьшее и наибольшее время из всех запусков.

В табл. 1 приведены значения временных затрат на реализацию одного шага по времени на двух, четырех и восьми процессорах соответственно. При увеличении размерности задачи N У T00 время выполнения алгоритма на двух и четырех процессорах

значительно увеличивается, что предполагает необходимость использования большего числа процессоров для решения задачи.

В табл. 2 и 3 приведены результаты ускорения и эффективности параллельного алгоритма МГЭ при его реализации на двух, четырех и восьми процессорах.

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

Для количества узлов N = 1600 на границе Г области Д и восьми процессоров в табл. 2 показано, что ускорение параллельного алгоритма равнялось 3.5 на обоих

Таблица 1. Временные затраты на реализацию одного временного шага

N КемГУ ИВТ СО РАН

2 4 8 2 4 8

400 1.0294 0.9427 1.0631 0.2152 0.2808 0.4350

700 4.0587 3.1010 2.9414 0.8256 0.6128 0.6919

1000 10.1185 7.1709 6.2155 2.2360 1.4797 1.3627

1300 21.5993 13.6901 11.1241 4.6265 2.9007 2.6842

1600 38.0553 23.5927 16.3035 8.3150 5.0353 4.6026

Таблица 2. Ускорение параллельного алгоритма

N КемГУ ИВТ СО РАН

2 4 8 2 4 8

400 1.0819 1.1814 1.0475 1.0929 0.8167 0.5405

700 1.3440 1.7591 1.8546 1.5141 2.0399 1.8066

1000 1.4688 2.0725 2.3911 1.7486 2.6424 2.8692

1300 1.4617 2.3062 2.8382 1.9258 3.0716 3.3194

1600 1.5321 2.4713 3.5763 1.9386 3.2013 3.5024

Таблица 3. Эффективность параллельного алгоритма

N КемГУ ИВТ СО РАН

2 4 8 2 4 8

400 0.5409 0.2953 0.1309 0.5464 0.2041 0.0675

700 0.6720 0.4397 0.2318 0.7570 0.5099 0.2258

1000 0.7344 0.5181 0.2988 0.8743 0.6606 0.3586

1300 0.7308 0.5765 0.3547 0.9629 0.7679 0.4149

1600 0.7660 0.6178 0.4470 0.9693 0.8003 0.4378

Таблица 4. Значения максимального заплеска

А 0.1 0.2 0.3 0.4 0.5 0.6

Теория 0.206 0.426 0.665 0.928 1.22 1.54

МГЭ 0.206 0.426 0.663 0.925 1.25 1.69

[11] 0.205 0.424 0.660 0.922 1.22 1.58

[13] 0.206 0.424 0.660 0.924 1.24 1.69

кластерах. На двух процессорах ускорение получено равным 1.5 на кластере КемГУ и 1.9 на кластере ИВТ СО РАН.

По результатам расчетов составлена табл. 3 эффективности Ер алгоритма. Для двух и четырех процессоров была получена достаточно хорошая эффективность распараллеливания алгоритма, равная 70-60 % (кластер КемГУ) и 80-90 % (кластер ИВТ СО РАН). Таким образом, можно говорить об эффективности применения параллельных алгоритмов для решения задач со свободными границами.

В табл. 4 приведены теоретические и численные результаты максимального заплеска уединенной волны на вертикальную стенку. Теоретические результаты получены по формуле [12]:

Численные результаты, полученные с помощью МГЭ, сравниваются с результатами из работ [11, 13].

Из табл. 4 видно совпадение результатов вплоть до второго знака. Если амплитуда волны A < 0.4, то относительная погрешность по сравнению с теоретическими результатами находится в пределах 1 %. Если же амплитуда волны A > 0.4, то относительная погрешность возрастает и составляет 10%, что видно для случая A = 0.6. Сравнение результатов c численными результатами, полученными другими авторами, показывает достаточно высокую точность МГЭ.

5. Численные результаты

Для исследования волновых режимов, порожденных различными законами движения оползня [14], проведена серия вычислительных экспериментов в полной нелинейной постановке на кластере КемГУ. При этом рассматривались следующие типы движения для двух видов оползней — слайда и слампа [2, 5]:

1) слайд 1 — разгон, равномерное движение, остановка, покой;

2) слайд 2 — разгон, равномерное движение, торможение, покой;

3) слайд 3 — разгон, равномерное движение;

4) сламп 1 — разгон, торможение, покой;

5) сламп 2 — разгон, торможение, покой.

Расчеты проводились в области с координатами x\ = 1.0 и xr = 41.0 до момента времени t = 50.0. Для записи характеристик волновых процессов были установлены семь виртуальных мареографов в точках с координатами: хмо = xco — 1, Хмi = xco = 2.38, xMi = xmj-i + 2, i = 2, 6. Для моделирования движения оползня использовались следующие значения параметров: Ah = 0.05, S = cos в/2, b = 1.0, H = 2.3, xc0 = 2.38,

Результаты сравнивались с результатами, полученными по нелинейно-дисперсионным (НЛД) моделям Мея—Меоте [15], Перегрина [16], Грина—Нагди [17], одно- и двухслойной моделям Лью—Линетта [8] и модели теории мелкой воды [3].

На рис. 2 показаны результаты сравнения записей седьмого мареографа, рассчитанных по полной модели и по НЛД-моделям с линейными дисперсионными членами. Как было показано в работе [4], учет линейной дисперсии приводит к существенному усложнению волнового режима при одновременном сохранении базовых характеристик. На рис. 2 можно видеть хорошее совпадение результатов расчетов по полной модели

(13)

в = 6°.

а б

Рис. 2. Результаты, рассчитанные для оползня, движущегося по закону "слайд 1", в седьмом мареографе по НЛД-моделям (жирные линии) с линейной дисперсией: а — модель Перегрина, б — модель Грина—Нагди; для сравнения приведены результаты, рассчитанные с использованием полной модели (тонкие линии)

аб

Рис. 3. Результаты, рассчитанные для оползня, движущегося по закону "слайд 1", в седьмом мареографе по НЛД-моделям (жирные линии) с нелинейной дисперсией: а — однослойная модель Лью—Линетта, б — двухслойная модель Лью—Линетта; для сравнения приведены результаты, рассчитанные с использованием полной модели (тонкие линии)

и моделям Перегрина и Грина—Нагди на начальных стадиях процесса. Расхождение результатов имеет место при T > 30.

На рис. 3 показаны результаты, полученные по НЛД-моделям с нелинейной дисперсией и по полной модели. Видно, что мареограммы для одно- и двухслойной моделей Лью—Линетта достаточно хорошо совпадают с результатами, полученными по полной модели. Различие результатов наблюдается по величине амплитуд. Отметим, что особенно хорошее совпадение кривых получено для полной и двухслойной моделей (рис. 3, б).

Также изучалась зависимость волновой картины от различных значений угла наклона в. На рис. 4 показаны мареограммы для полной модели и двухслойной модели Лью—Линетта. Угол в равнялся 6, 9, 12 и 15 градусам. Закон движения "слайд 3" задавался следующим образом: s(t) = a0t2/2 до момента времени ti, s(t) = u(t — t1) + a0t\/2 для t > t1, где s(t) — путь, который оползень прошел по склону: xc(t) = xc0 + s(t) cos в;

а0 = 0.3$ 8т(в), и = 1.16^/6д 8т(в). Когда глубина (при невозмущенной жидкости) над центром тела хс достигает значения 1.44, оползень останавливается. Значения основных параметров выбирались следующие: ДЛ = 0.05, 6 = 1.0, Н = 2.3, д = 1.0. Начальное заглубление оползня й для разных углов в было одинаковым и равнялось 0.44.

в = 6° : Хсо = 5.2, ¿1 и 11.96, ¿^ор и 30.0;

в = 9° : Хсо = 3.47, ¿1 и 9.77, и 18.0;

в = 12° : Хсо = 2.59, ¿1 и 8.48, ¿^ор и 12.8

в = 15° : Хсо = 2.05, ¿1 и 7.6, ¿^ор и 10.0.

Анализ численных результатов, полученных для полной и двухслойной моделей, показал, что мареограммы достаточно хорошо совпадают в прибрежной зоне при в = 6°, но в мористой зоне расходятся по амплитудам мареограмм. Увеличение угла наклона в приводит к расхождению результатов как в мористой, так и в прибрежной зонах (рис. 4). В этом случае для полной модели характерно затухание волновых процессов, а для модели Лью—Линетта — увеличение амплитуды.

На рис. 5 приведено сопоставление численных результатов решения полных уравнений гидродинамики идеальной жидкости с помощью МГЭ и метода конечных разностей (МКР) на адаптивных сетках [4]. Мареограммы показаны для второго (рис. 5, а) и седьмого (рис. 5, б) мареографов при следующих параметрах: ДЛ = 0.05, Хсо = 2.38, в = 6°, 6 = 1.0. Расчеты проводились в полной нелинейной постановке для законов движения "слайд 1", "сламп 1", "сламп 2".

Из рис. 5 видно, что в прибрежной зоне методы дают близкие результаты. Однако для МГЭ волна понижения больше, чем для МКР. В случае, когда тело движется по закону "слайд 1" и "сламп 2", в мористой зоне методы показывают достаточно хорошее совпадение до момента остановки оползня Т = 30.0 и различаются в поведении "хвоста" мареограмм. Отметим, что амплитуда мареограмм для метода граничных элементов больше, чем для МКР. Для закона движения тела "сламп 1" как в прибрежной, так и в мористой зоне получено хорошее совпадение результатов. Такое поведение мареограмм можно объяснить тем, что закону движения "сламп 1" характерно кратковременное движение тела.

Рис. 4. Результаты, рассчитанные для оползня, движущегося по закону "слайд 3", в = 15° в первом (а) и седьмом (б) мареографах: жирная линия — двухслойная модель Лью—Линетта, тонкая линия — полная модель

"слайд 1"

0.01

-0.01

10 20 30 40 50

10 20 30 40 50

"сламп 1"

"сламп 2"

а б

Рис. 5. Мареограммы, рассчитанные по полной модели МГЭ (сплошная линия) и МКР (штриховая линия) во втором (а) и седьмом (б) мареографах

Для изучения влияния основных параметров задачи (толщины Ак, длины Ъ, заглубления й оползня и угла в) на характеристики волновых режимов была проведена серия расчетов в полной нелинейной постановке методом граничных элементов. На рис. 6, а показаны мареограммы для первого мареографа при следующих параметрах: толщина оползня равна Ак = 0.01 (штриховая линия) и Ак = 0.05 (сплошная линия), начальное положение центра масс хс0 = 2.38, протяженность оползня Ъ =1. На рис. 6, б приведены мареограммы для первого (сплошные линии) и четвертого (штриховые линии) мареографов при следующих параметрах: протяженность оползня равна Ъ = 1 (сплошные и штриховые линии без маркеров) и Ъ = 3 (сплошные и штриховые линии с маркерами), хс0 = 4.38, Ак = 0.05.

Рис. 6. Зависимость от толщины (а) и от протяженности (б) оползня

Рис. 7. Режимы обрушения волн для закона движения "слайд 1": а — АЬ = 0.1, хс = 2.38, в = 6°, Ь = 1; б — АЬ = 0.2, хс = 3.38, в = 10°, Ь = 3

а

Видно, что увеличение толщины и протяженности оползня приводит к росту амплитуды мареограмм, при этом возможно возникновение сильно нелинейных деформаций свободной границы, таких как "опрокидывание" волн (рис. 7). При решении задачи МГЭ происходит перехлест узлов, что приводит к прерыванию расчета. Сначала формируется волна понижения, следующая за оползнем, затем образуется отраженная волна повышения, которая движется в мористом направлении. В вершине этой волны формируется сгусток, переходящий в скользящий бурун (рис. 7).

Заключение

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

1) создание параллельного алгоритма метода граничных элементов для моделирования генерации поверхностных волн движением затопленного оползня;

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

Авторы выражают признательность д-ру физ.-мат. наук, профессору К.Е. Афанасьеву, канд. физ.-мат. наук З.И. Федотовой и д-ру физ.-мат. наук, профессору Л.Б. Чу-

барову за внимание к работе.

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

[1] Grilli S.T., Watts P. Modeling of waves generated by moving submerged body. Applications to underwater landslide // Eng. Analysis With Boundary Elements. 1999. Vol. 23. P. 645-656.

[2] Watts P., Imamura F., Grilli S. Comparing model simulations of three benchmark tsunami generation cases // Sci. of Tsunami Hazards. 2000. Vol. 18, N 2. P. 107-123.

[3] Чубаров Л.Б., Федотова З.И., Елецкий С.В. Численное моделирование генерации волн движением оползня // Тр. Междунар. конф. по вычисл. математике. Новосибирск: ИВМиМГ СО РАН, 2004. Ч. 2. С. 753-758.

[4] Шокин Ю.И., Федотова З.И., ХАкимзянов Г.С. и др. Моделирование генерации цунами движением оползня с учетом вертикальной структуры течения // Современные методы математического моделирования природных и антропогенных катастроф: Тр. 8-й Всерос. конф. (Кемерово, 26-28 октября 2005 г.). Кемерово, 2005. С. 20-40.

[5] Watts P., Grilli S.T., Kirby J.T. et al. Landslide tsunami case using a Boussinesq model and fully nonlinear tsunami generation model // Natural Hazards and Earth System Sci. 2003. Vol. 3, N 5. P. 391-402.

[6] Афанасьев К.Е., Гудов А.М. Информационные технологии в численных расчетах: учеб. пособие. Кемерово: Изд-во КемГУ, 2001.

[7] Бребия К., Теллес Ж., Вроубел Л. Методы граничных элементов. М.: Мир, 1987.

[8] Lynett P.A., Liu P.L-F. Numerical study of submarine landslide generated waves and runup // Proc. Royal Society of London A. 2002. Vol. 458. P. 2885-2910.

[9] Афанасьев К.Е., Стуколов С.В. Многопроцессорные вычислительные системы и параллельное программирование: учеб. пособие. Кемерово: Изд-во КемГУ, 2003.

[10] Афанасьев К.Е., Стуколов С.В. О наличии трех решений при обтекании препятствий сверхкритическим установившимся потоком тяжелой жидкости // ПМТФ. 1999. № 1. С. 27-35.

[11] Протопопов Б.Е. Численный анализ трансформации уединенной волны при отражении от вертикальной преграды // Изв. АН. СССР. Сер.: МЖГ. 1990. № 5. С. 115-123.

[12] Su C.H., Mirie R.M. On head-on collision between solitary waves //J. Fluid Mech. 1980. Vol. 98. Pt. 3. P. 509-525.

[13] Kim J.W., Bai K.J., Ertekin R.C., Webster W.C. A strongly-nonlinear model for water waves in water of variable depth — the irrotational Green-Naghdi model // J. of Offshore Mech. and Arctic Engin. 2003. Vol. 125. P. 25-32.

[14] Шокин Ю.И., Чубаров Л.Б. О подходах к численному моделированию оползневого механизма генерации волн цунами // Вычисл. технологии. 2006. Т. 11. Специальный выпуск, посвященный 85-летию со дня рождения академика Н.Н. Яненко. Ч. 2. С. 100-111.

[15] Mei C.C., Le Mehaute B. Note on the equations of long waves on uneven bottom //J. Geophys. Res. 1966. Vol. 72, N 2. P. 815-827.

[16] Peregrine D.H. Long waves on a beach // J. Fluid Mech. 1967. Vol. 27, Pt. 4. P. 815-827.

[17] Green A.E., Naghdi D.M. A derivation of equations for wave propagation in water at variable depth // J. Fluid Mech. 1976. Vol. 78, Pt. 2. P. 237-246.

Поступила в редакцию 15 января 2008 г.