Содержание к диссертации
Введение
Глава 1. Математические и вычислительные алгоритмы 20
1 Общая постановка задач 20
1.1 Уравнение неразрывности 20
1.2 Уравнения движения 21
1.3 Постановка нестационарной задачи 22
2 Метод граничных элементов 25
2.1 Вычисление интегралов 29
2.2 Решение систем линейных алгебраических уравнений 30
2.3 Алгоритм движения по времени 32
3 Кинематические и динамические характеристики 35
3.1 Вычисление компонент вектора скорости 35
3.2 Вычисление гидродинамических характеристик 36
Глава 2. Взаимодействие поверхностных волн с препятствием 40
1 Тестирование вычислительных алгоритмов 40
1.1 Тестирование МГЭ методом пробных функций 40
1.2 Нестационарное движение уединенной волны по бассейну с ровным дном 41
1.3 Накат солитона на вертикальную стенку 43
1.4 Движение уединенной волны над прямоугольным выступом . 46
2 Взаимодействие поверхностных волн с частично погруженным в жидкость телом 49
2.1 Постановка задачи 49
2.2 Численные результаты 50
3 Численное моделирование взаимодействия солитона с подводным препятствием 63
3.1 Постановка задачи 63
3.2 Численные результаты 64
Глава 3. Численное моделирование генерации поверхностных волн движе нием оползня 76
1 Схема модельной области и механизмы движения оползня 76
1.1 Постановка задачи 76
1.2 Схемы движения оползня 78
2 Численные результаты 81
Глава 4. Информационные технологии в численных расчетах 96
1 Реализация параллельного метода граничных элементов 96
1.1 Эффективность и ускорение . 97
1.2 Схема последовательного алгоритма метода граничных элементов и его распараллеливание 98
1.3 Алгоритм решения СЛАУ методом Гаусса 102
1.4 Тестирование параллельного алгоритма 106
2 Информационная система сопровождения численного эксперимента 111
2.1 Структура информационной системы 111
2.2 Интерфейс обмена данными 112
2.3 Логическая схема базы данных 116
2.4 Хранилище данных 119
2.5 Оболочка информационной системы 121
Литература
- Уравнения движения
- Нестационарное движение уединенной волны по бассейну с ровным дном
- Схемы движения оползня
- Схема последовательного алгоритма метода граничных элементов и его распараллеливание
Введение к работе
Диссертационная работа посвящена исследованию двумерных течений идеальной однородной несжимаемой жидкости с поверхностными гравитационными волнами методом граничных элементов, и вопросам эффективного использования современных вычислительных технологий и методов параллельного программирования.
В современной науке широко используется методология математического моделирования [76]. Сущность этой методологии состоит в замене исходного объекта его математической моделью и в исследовании современными вычислительными средствами математических моделей. Этот метод сочетает в себе многие достоинства как теории, так и эксперимента. Работа с моделью объекта дает возможность, относительно быстро и без существенных затрат исследовать его свойства и поведение. Широкое применение вычислительного эксперимента с математическими моделями объектов позволяет, опираясь на современные вычислительные методы и технические ресурсы, подробно и глубоко изучать объекты в достаточно полном объеме, что недоступно аналитическим подходам.
Широкий круг проблем, стоящих перед современной наукой и техникой, связан с решением уравнений механики жидкости. Большинство течений жидкости имеет природное (океаны, моря, ветер) и техногенное происхождение (самолеты, машины, биоинженерия). Существует потребность в моделировании проблем течения жидкостей, с целью лучшего понимания сложных явлений и повышения качества технологий. Результаты исследований течений жидкости со свободными поверхностями находят многочисленные технические приложения, такие как физическая океанология, гидротехника, кораблестроение и других. Эти задачи традиционно считаются непростыми, поскольку к нелинейности краевой задачи добавляется еще и дополнительная сложность, связанная с определением заранее неизвестной формы свободной
границы. Примерами таких течений являются: нестационарное движение волн над неровным дном, выход волн на мелководье, эволюция свободной поверхности под действием силы тяжести, движение тел в жидкости, взаимодействие поверхностных волн с препятствиями, распространение волн цунами и т.д. Учет подобных явлений усложняет математическую постановку задач и порождает самостоятельные проблемы при их решении. Однако и в этом случае основу задачи составляют классические уравнения механики жидкости.
Фактически единственным эффективным способом решения сложных систем нелинейных дифференциальных уравнений, являются численные методы, реализуемые на быстродействующих вычислительных машинах. То есть, на основе математической модели при помощи непосредственного численного решения соответствующих уравнений количественно определяется поведение течений жидкости в тех или иных условиях.
В отличие от аналитических [73] и инженерных [78,79] методов, где зачастую для каждой задачи разрабатываются свои самостоятельные приемы решения, численные методы отличаются большой универсальностью и применимы для исследования широкого класса явлений.
Среди множества численных методов, применяемых для решения задач потенциальных течений со свободными границами, широкое развитие получил метод граничных элементов (МГЭ). Он составил удачную конкуренцию таким популярным среди исследователей методам, как метод конечных разностей (МКР) [77] или метод конечных элементов (МКЭ) [56,86,87]. Привлекательность МГЭ обусловлена, прежде всего тем, что в МКР, как и МКЭ, требуется разбиение всей области течения, в то время как в МГЭ дискретизации подвергается лишь граница области. Для реализации такой возможности в МГЭ требуется переход от исходной краевой задачи для дифференциальных уравнений, описывающих некоторый процесс, к соотношениям связывающим неизвестные функции на границе области.
Полный обзор технологии метода граничных элементов можно найти в
5.
монографиях П. Беиерджи, Р. Баттерфилда [28] и К. Бреббии, Ж. Теллеса, Л. Вроубела [27].
Обзор состояния проблемы и методы решения различных задач идеальной жидкости со свободными границами изложены в монографиях М.И. Гу-ревича [43], О.М. Киселева, Л.М. Котляра [53]. Теории волновых движений жидкости посвящена монография Л.Н. Сретенского [80]. В этих работах используются, в основном, аналитические методы, которые применимы лишь для ограниченного круга задач.
Если жидкость весомая, а на свободной границе присутствуют сильно нелинейные деформации, то применение аналитических методов практически невозможно. В этом случае используются различные модификации численно -аналитических методов (эти направления нашли развитие в работах В.П. Жит-никова [47^9], Д.В.; Маклакова [61-65]; Р.А. Рузиева, Г.С. Хакимзянова [74]).
Решению задач в точной постановке, выяснению особенностей и разработке методов исследования посвящена монография A.M. Лаврентьева, Б.В. Шабата [57].
Многие из задач течения идеальной однородной несжимаемой жидкости описываются уравнением Лапласа с нелинейными условиями на свободной поверхности. Значительное место в этих задачах занимает волновая тематика. В литературе приводится немалое количество примеров уединенных волн (солитонов). Отметим некоторые из них, полученные аналитически по различным линейным и приближенным нелинейным теориям с помощью численного анализа точных и приближенных нелинейных уравнений: Д.В. Маклаков [63,64], Б.Е. Протопопов [72], Препринт по ред. Ю.И. Шокина [92], М. Tanaka [131], Ан.Г. Марчук, Л.Б. Чубаров, Ю.И. Шокин [60], Е.А. Karabut [115]; а также путем моделирования волн различными подвижками боковых стенок, дна, или с помощью создания локального возвышения уровня жидкости; A.M. Франк [98], найденные в эксперименте: СВ. Манойлин [59], В.И. Букреев, Н.П. Туранов [30]. Обзор методов численного решения стационарных и неста-
ционарных задач со свободными границами приведен в работе И.В. Стуро-вой [84].
Подход к численному решению задачи о течении идеальной несжимаемой жидкости с поверхностными волнами в каналах со сложным очертанием берегов методом конечных-разностей рассматривается в работах [22,101]. В статье [102] предложен итерационный алгоритм расчета на криволинейных адаптивных сетках стационарных течений жидкости с поверхностными гравитационными волнами в речных руслах с островами в рамках модели мелкой воды с учетом неровности дна. Обзор работ, посвященных разработке конечно-разностных методов на адаптивных сетках [58] для расчета течений идеальной жидкости со свободной границей в рамках нелинейных моделей мелкой воды и модели потенциальных течений дан в работе Г.С. Хакимзянова [99].
Исследование задач взаимодействия поверхностных волн с препятствиями является одним из наиболее интересных приложений нелинейной гидродинамики в связи с важностью вопросов по определению воздействия этих волн на гидротехнические сооружения (плавучие доки, волноломы, платформы и т.д.) [114,116] и акватории портов.
Вопросам взаимодействия поверхностных волн с препятствиями посвящено множество теоретических и экспериментальных работ: взаимодействию с вертикальной стенкой [72,108,111,117,122,130], с подводным препятствием (ступенька, выступ, полуцилиндр): теоретические и экспериментальные исследования [128], численное моделирование выполнено в [31-33,74,92,97,104, 107,109,118,126]. В работе [128] - использовалась неявная конечно разностная схема для нелинейно-дисперсионной модели мелкой воды, в [97] - дискретная модель несжимаемой жидкости, в [74,92,104,126] - модель потенциальных течений, в [31,32] - одномерная модель мелкой воды.
В работе [33] в рамках одномерных моделей мелкой воды второго приближения разработан метод расчета транскритических течений над неровным дном, позволяющий учесть опрокидывание волн и возникновение приповерх-
ностного турбулентного слоя. Для нелинейно дисперсионных моделей мелкой воды в работе [109] рассматривается решение задачи о движении уединенной волны в канале с криволинейным дном. Обсуждаются результаты сравнительного анализа основных свойств волны и применимость моделей Перегрина, Железияка-Пелиновского, Ким-Рейд-Витакера, Федотовой-Пашкова и классических уравнений линейной и нелинейной теории мелкой воды.
Для модели плоских потенциальных течений идеальной несжимаемой жидкости со свободной границей в работах [100,103] методом конечных разностей предложен подход к решению задачи о взаимодействии солитона с подводным или частично погруженным препятствием.
В полной нелинейной постановке изучены различные режимы обтекания подводных препятствий, например, обтекание полуцилиндра [3,4,19], взаимодействие уединенной волны с вертикальной и наклонной стенкой [18,21,81], взаимодействие уединенной волны с частично или полностью погруженным в жидкость препятствием [8,9].
В диссертационной работе проводятся исследования о генерации поверхностных волн (цунами) движением оползня. Многообразие причин, порождающих волны цунами, обилие факторов, определяющих характер их трансформации при распространении по океану и в береговой зоне, выделяют феномен цунами и обуславливают необходимость комплексного подхода к его изучению.
Основные направления исследований, сформулированы в работах Е.Ф. Са-варенского [75], развитые затем в трудах С. С. Войта [35], А.С. Алексеева [2], Ю. И. Шоішна [60,92]. Обзор посвященный проблеме цунами дан в работе Л.А. Островского и Е.Н. Пелиновского [69], Н.А. Щетникова [94].
В работах [113,132] показано, что для изучения закономерностей волнообразования необходимо исследовать зависимость характеристик процесса от основных параметров задачи: длины и ширины оползня, глубины его заглубления и закона его движения. Такой подход при соответствующей параметри-
зации представляет адекватную схематическую картину реальных оползневых процессов в широком диапазоне изменения определяющих характеристик.
В статьях [46,90] исследовалась возможность использования приближенных математических моделей гидродинамики для моделирования механизма движения оползня, при этом анализировалась необходимость учета вертикальной структуры течения. Было показано, что для детального и качественного описания явлений в обширных водоемах следует использовать модели волновой гидродинамики, хорошо описывающие дисперсию и отражающие неоднородность процесса в вертикальном направлении. Для полной модели, в работах [46,132,133] получены результаты наиболее близкие к лабораторным экспериментам.
В работе [93] выполнен комплекс многопараметрических расчетов с помощью иерархии моделей волновой гидродинамики, учитывающих изменение во времени донной поверхности: линейные и нелинейные модели мелкой воды, слабо нелинейные дисперсионные модели, полученные в [45] и совпадающие в случае ровного дна с моделями Мея - Меоте и Перегрина [123,127], упрощенные варианты модели Грина - Нагди [112] и Нвогу [125], одно -двухслойные нелинейно-дисперсионные модели (НЛД) Лью - Линетта [120], а также полная модель течения течения идеальной жидкости [100]. Более полное описание моделей приведено в работе [91].
Появление быстродействующих ЭВМ и создание вычислительных методов внесло новую струю в исследования задач со свободными границами, которые требует для своего решения значительных вычислительных ресурсов. Повышенные требования к производительности и памяти обусловлены сложными нелинейными моделями среды, описываемыми большим числом уравнений, пространственным характером задачи и нестационарностью протекающих процессов [18,37,129]. Многим из них требуется высокое быстродействие, а также необходима обработка и хранение большого объема информации, что предъявляет повышенные требования к вычислительным системам.
В 70-80-х годах в России велись исследования, направленные на создание параллельных вычислительных систем. Примерами таких систем являются PHOENIX [1], ПС-2000 [71], а так же многопроцессорные вычислительные комплексы Эльбрус [55]. Принципы, заложенные в основу структурной организации упомянутых машин, находят свое применение и в настоящие время. Одновременно с разработкой параллельных вычислительных систем учеными велись работы по распараллеливанию алгоритмов сложных задач, например, [34,38] посвященные общим вопросам распараллеливания алгоритмов, [41,68] - распараллеливанию численных методов линейной алгебры. В работах [51,70,95] - рассматриваются проблемы численного решения краевых задач для обыкновенных дифференциальных уравнений и уравнений в частных производных.
Одно из наиболее интересных направлений развития современной вычислительной техники на сегодняшний день представляют многопроцессорные вычислительные системы [20,38,105].
Вычислительные системы (кластеры) представляют собой мультикомпью-теры, состоящие из множества отдельных компьютеров или рабочих станций общего назначения (узлов), связанных между собой единой коммуникационной системой. Каждый узел имеет свою локальную оперативную память (общей физической оперативной памяти для узлов не существует). Для вычислительных кластеров используются стандартные для рабочих станций операционные системы, чаще всего, свободно распространяемые - Limrx/FreeBSD, вместе со специальными средствами поддержки параллельного программирования и распределения нагрузки. Программирование, выполняется на основе модели передачи сообщений (MPI) [85]. Привлекательность использования вычислительных кластеров для решения задач со свободными границами рассматриваются в работах [5,6,20,83].
Численный решение задач гидродинамики не ограничивается только применением современных вычислительных технологий и методов параллельного
программирования. Большую важность при численном моделировании принимает корректность и простота ввода начальных данных задачи, а также обработка результатов численного расчета [14,15].
Поэтому актуальной задачей является разработка полнофункциональной информационной системы [12], с возможностью автоматизированной подготовки данных для численного эксперимента, хранения и систематизации, как полученных результатов, так и самих постановок задач, входных данных для каждого из расчетов и необходимых алгоритмов.
Основу численного эксперимента составляет триада модель - алгоритм -программа [76].
На первом этапе вычислительного эксперимента строится модель исследуемого объекта, отражающая в математической форме важнейшие его свойства.
Второй этап связан с выбором вычислительного алгоритма для реализации модели на компьютере.
На третьем этапе создается программное обеспечение для реализации модели и алгоритма на компьютере.
Опираясь на триаду модель - алгоритм - программа, исследователь получает в руки универсальный, гибкий и недорогой инструмент, который вначале отлаживается, тестируется и калибруется на решении содержательного набора пробных задач.
В настоящие время существует ряд пакетов прикладных программ (ППП) предназначенных для решения научных, инженерных и прикладных задач. Описание некоторых пакетов можно найти в работах [36,40]. Среди систем автоматизирующих полный цикл решения некоторой задачи можно выделить:
ANSYS () - предназначенный для исследования: за
дач статического и динамического анализа конструкций с учетом геомет
рической и физической нелинейности; задач ползучести и пластичности;
задач линейной и нелинейной устойчивости конструкций; стационарных
и нестационарных задач теплофизики с учетом фазового перехода; задач гидро-газодинамики; электромагнитных полей (в т.ч. высокочастотный анализ); задач акустики; связанных задач (например, взаимодействие жидкости с конструкцией), BEASY (). Пакет состоит из четырех основных модулей: Mechanical Design - предназначен для решения задач механики; Fatigue and Crack Growth - анализ усталости материалов и процесса образования трещин; Acoustic Design - решение задач акустики; Corrosion and Cathodic Protection - задачи коррозии и защиты от коррозии. О предмете и содержании диссертации.
Диссертация состоит из введения, четырех глав, приложения и списка литературы.
Диссертация посвящена численному исследованию двумерных нестационарных задач потенциальных течений идеальной однородной несжимаемой жидкости со свободными границами с использованием новых информационных технологий. Рассматриваются традиционные для этой области исследования:
взаимодействие уединенной волны с частично или полностью погруженным в жидкость препятствием прямоугольного сечения,
генерация поверхностных волн движением оползня, в зависимости от типа движения и геометрических параметров,
применение современных информационных технологий для разработки инструментария автоматизации численного эксперимента и средств хранения информации.
Решение задач со свободными границами в полной нелинейной постановке выполнялось численным методом граничных элементов.
Цель диссертации заключается в решении двумерных нелинейных задач гидродинамики идеальной однородной жидкости со свободными границами,
а также построение эффективного параллельного численного алгоритма метода граничных элементов, и автоматизации обработки результатов численных расчетов.
Содержание работы
Первая глава посвящена некоторым аспектам использования МГЭ для решения нелинейных задач движения жидкости со свободными границами.
В первом параграфе приводятся постановка плоской нестационарной задачи идеальной однородной несжимаемой жидкости со свободными границами.
Второй параграф посвящен краткому описанию метода граничных элементов основанном на третьей формуле Грина. Рассматриваются вопросы сведения исходной дифференциальной задачи к интегральным соотношениям, заданным на границе области, вычисления полученных интегралов и построения матриц.
В третьем параграфе приводятся формулы высокого порядка точности для дифференцирования функций, заданных на границе области. Рассматривается алгоритм решения нестационарной задачи, и алгоритм автоматического выбора шага по времени. Для движения по времени и выбор шага по времени используется метод Эйлера. Данный алгоритм движения по времени впервые был предложен A.M. Лаврентьевым и Б.В. Шабатом [57] и хорошо зарекомендовал себя при решении задач со свободными границами [3,7-9,54,106]. Тестовые расчеты, представленные в первом параграфе второй главы, демонстрируют достаточную точность используемого метода расчета для моделирования нестационарных волновых задач.
Приводятся алгоритмы вычисления гидродинамических и энергетических характеристик - давления, массы кинетической и потенциальной энергии.
Во второй главе в полной нелинейной постановке приводится решение двумерных нелинейных задач гидродинамики идеальной однородной несжимаемой жидкости со свободными границами, на примере которых демонстрируется эффективность алгоритмов и методов.
Вопросам взаимодействия солитонов с препятствиями посвящено большое количество теоретических и экспериментальных работ в связи с важностью вопросов по определению воздействия этих волн на гидротехнические сооружения и акватории портов. При проведении экспериментов возникают существенные трудности, связанные с тем, что волнопродуктор создает волну с дисперсионным хвостом из волн малой амплитуды, для удаления которых приходится использовать специальные волногасители, которые вместе с гашением дисперсионного хвоста могут изменить амплитуду и форму основной (рабочей волны). Кроме того, регистрация проистекающих процессов требует специальной аппаратуры и является задачей сложной и дорогостоящей. Поэтому разработка и тестирование численных алгоритмов для решения задач взаимодействия солитонов с препятствиями является актуальной задачей.
Задачи о взаимодействии уединенной волны с частично или полностью погруженным в жидкость телом выполнены в рамках договоров ИВТ СО РАН и совместной лабораторией Кемеровского государственного университета ИВТ СО РАН.
В первом параграфе на ряде тестовых задач демонстрируется эффективность применения методов граничных элементов для решения задач движения идеальной несжимаемой жидкости.
Рассматриваются следующие задачи:
тестирование методом пробных функций;
распространение уединенной волны по бассейну с ровным дном;
взаимодействие уединенной волны с вертикальной твердой преградой;
взаимодействие уединенной волны с подводной ступенькой. Точность аппроксимации интегральных уравнений демонстрируется при помощи теста предложенного Петровым А.Г. и Смоляниным В.Г. в работе [66]. Здесь найденные значения искомой функции сравниваются с известным аналитическим решением. Показывается сходимость решения к точному при увеличении количества точек разбиения.
Хорошим критерием проверки качества численного метода является задача о движении уединенной волны по ровному дну, поскольку известно, что во все время движения волна должна сохранять свои основные характеристики: амплитуду, скорость, массу. Сохранение характеристик уединенной волны зависит не только от точности численного метода, но и от качества начального приближения. Б случае задания начальной формы волны из известных приближений (Бусинеска-Лейтона, Овсянникова и т.п.) во время ее движения образуется заметный диспергирующий волновой цуг волн малой амплитуды [3,18,72]. Эти волны,оказывают влияние на характеристики основной волны (амплитуду, скорость волны, массу и т.д.). Для уменьшения таких эффектов, начальное приближение уединенной волны выбирается из численного решения стационарной задачи [18,19,82].
На примере задачи о взаимодействии уединенной волны с вертикальной твердой преградой и с подводной ступенькой демонстрируется эффективность применения прямого метода Эйлера с автоматическим выбором шага по времени для решения нестационарной задачи.
Во втором параграфе представлены результаты расчетов движение жидкости при набегании уединенной волны на тело прямоугольного сечения, частично погруженного в жидкость. В работе Г.С. Хакимзянова, Ю.И. Шокина, В.Б. Барахнина, Н.Ю. Шокиной [100] данная задача исследовалась с помощью конечно - разностных методов расчета на адаптивных сетках.
В данной работе моделирование осуществлялось методом граничных элементов и модели потенциальных течений идеальной несжимаемой жидкости. Исследовано влияние ширины и заглубления тела, а также амплитуды набегающей волны на высоту заплеска, волновую картину перед телом и за ним, динамическую нагрузку на тело.
В третьем параграфе представлены результаты исследования трансформации волн, распространяющихся над препятствием, используется метод граничных элементов и модели потенциальных течений идеальной несжимаемой
жидкости. Постановка задачи о распространении солитона над выступом расположенным на дне была взята из работы М.Г. Хажоян, Г.С. Хакимзянова [103], где показано, что для изучения закономерностей волнообразования необходимо исследовать зависимость характеристик процесса от основных параметров задачи: протяженности, высоты и амплитуды набегающей волны.
Проводиться анализ численных результатов, при движении солитона над подводным препятствием прямоугольного сечения, расположенного на горизонтальном дне. Изучается изменение динамической нагрузки при взаимодействии жидкости с твердыми границами. Исследуется влияние основных параметров на характеристики течения.
В третьей главе рассматриваются результаты исследования генерации поверхностных волн движением оползня. Задача о генерации поверхностных волн движением оползня выполнена в рамках интеграционного проекта фундаментальных исследований Объединенного ученого совета по механике и энергетике СО РАН (2006-2008) по теме "Численное моделирование нестационарного взаимодействия сложных упругих конструкций с жидкостью или газом".
Первый параграф содержит постановку задачи о генерации поверхностных волн цунами движением оползня. Рассматриваются разные законы движения оползня [91] называемые слайдами и слампами [1.33].
Во втором параграфе рассматривается решение задачи о генерации поверхностных волн движением оползня. Приводятся результаты сравнения значений мареограмм изменения уровня свободной поверхности полученные с помощью алгоритмов, основанных на иерархии моделей волновой гидродинамики, включающей уравнения теории мелкой воды в приближениях, учитывающих нелинейные и дисперсионные эффекты [91], и полные уравнения гидродинамики идеальной жидкости [100] решенные с помощью конечно -разностных методов расчета на адаптивных сетках..
Четвертая глава посвящена обсуждению ряда вопросов построения эф-
фективных алгоритмов ускорения и автоматизации решения задач со свободными границами. Получено свидетельство об официальной регистрации программы "Пакет приоадных программ для задач идеальной несжимаемой жидкости со свободными границами (AKORD)" для ЭВМ №2002611284 от 31 июля 2002г. (РОСПАТЕНТ).
В первом параграфе рассматриваются обсуждаются вопросы геометрической декомпозиции расчетной области и балансировки нагрузки процессоров, способы распределения данных по процессорам, а также особенности параллельной реализации метода граничных элементов, применяемого для решения задач со свободными границами. Приводятся результаты тестирования параллельного алгоритма метода граничных элементов на кластере Кемеровского государственного университета и Института вычислительных технологий СО РАН.
Во втором параграфе приводится описание проблемно - ориентированной оболочки предназначенной для информационной поддержки всех основных этапов численного эксперимента, начиная с постановки, и заканчивая графическим анализом полученных результатов. Приводятся форматы файлов обмена данными между приложениями системы. Рассматриваются вопросы создания удаленной базы данных расчетов необходимой для хранения представляющей интерес информации.
Основные результаты, выносимые на защиту:
1, Решение нестационарных задач волновой гидродинамики при взаимодействии нелинейных поверхностных волн с препятствиями, частично или полностью погруженными в жидкость. Установлено влияние амплитуды волны А, высоты (или зазора) h и протяженности а препятствия на следующие характеристики: максимальный заплеск, амплитуда прошедшей волны Ор, амплитуда отраженной волны а0 и нагрузки на твердых границах Ps. Для различных значений амплитуды А набегающей водны, протяженности препятствия а и высоты (или зазора) h получены табли-
цы изменения параметров максимального заплеска, амплитуды прошедшей а0 и отраженной волны ар, а также изменения динамической нагрузки Ps. Выявлены эффекты сильной деформации свободной границы "опрокидывание волн" в зависимости от геометрических характеристик препятствия.
Решение нестационарной задачи о генерации волн цунами движением оползня для разных законов движения и его геометрических параметров. Получены волновые режимы в зависимости от начального залегания d, толщины h, протяженности Ь и законов движения оползня. Установлено, что изменение геометрических параметров и схемы движения оползня, приводит к разным картинам формирования цуга волн различной амплитуды, движущихся к берегу и от него. Выявлены эффекты "опрокидывания волн" в зависимости от геометрических характеристик оползня.
Параллельная реализация метода граничных элементов, применяемого для решения задач механики жидкости.
Проблемно - ориентированная оболочка для информационной поддержки вычислительного эксперимента для задач идеальной несжимаемой жидкости со свободными границами.
Результаты работы по мере их получения докладывались и обсуждались на V Сибирской школе-семинаре "Математические проблемы механики сплошных сред" (Новосибирск, 2001г.), Международной конференции молодых ученых по математическому моделированию и информационным технологиям (Новосибирск, 2002г.), IV Международной конференции молодых ученых по математическому моделированию и информационным технологиям (Красноярск, 2003г.), V Международной конференции молодых ученых по математическому моделированию и информационным технологиям (Новосибирск, 2004г.), I региональной научно-практической конференции "Информационные недра Кузбасса" (Кемерово 2001г.), II региональной научно-практической конференции "Информационные недра Кузбасса" (Кемерово 2003г.), III региональной
научно-практической конференции "Информационные недра Кузбасса" (Кемерово 2004г.), Международной конференции "High Speed Hydrodynamics" (Чебоксары 2004г.), VI Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям (с участием иностранных ученых) (Кемерово 2005г.), IX Международной летней научной школе "Гидродинамика больших скоростей и численное моделирование" (Кемерово 2006г.), на научном семинаре "Численные методы решения задач механики сплошной среды" Кемеровского университета (Кемерово, 2000-2006гг.).
По теме диссертации опубликовано 11 работ, в том числе (в скобках в числителе указан общий объем этого типа публикации, в знаменателе - объем, принадлежащий автору) 2 статьи в изданиях, рекомендованных ВАК для представления кандидатской диссертации (1,6/1,2 печ. л.), 9 - в трудах международных и всероссийских конференций (3,3/1,6 печ. л.). Основные результаты диссертации опубликованы в работах [8-11,14,15,23-26,106]. При выполнении работ [14,15], опубликованных совместно с научным руководителем и другими соавторами, автор диссертации принимал участие в разработке и реализации проблемно-ориентированной оболочки, интерфейса обмена данными и базы данных расчетов. В соавторстве с научным руководителем [8-11,106] автору диссертации принадлежит участие в постановке задачи, исследование численных алгоритмов, обсуждении полученных результатов, подготовке и представлении статей и докладов на конференциях. Автором выполнена программная реализация численных алгоритмов, проведены расчеты тестовых задач и значительный цикл вычислительных экспериментов.
Автор выражает глубокую признательность научному руководителю - доктору физико-математических наук, профессору К.Е. Афанасьеву за постоянное внимание, творческие идеи и огромную помощь при выполнении работы.
Уравнения движения
Рассмотрим идеальную однородную несжимаемую жидкость, находящуюся в однородном поле сил тяжести. Предположим, что течение является потенциальным, т.е. для вектора скорости v — и(х, t) существует функция (х, і) такая, что v = gradip. Здесь х = х(ж,т/) - радиус-вектор точки области течения D. Тогда из (1.6) получаем, что для потенциального течения однородной несжимаемой жидкости внутри и на границе области (D — Г U D) справедливо уравнение Лапласа Др(х,г) = 0, xeD. (1.7)
Если движение происходит под действием сил тяжести, направленных противоположно оси у, то давление Р, плотность р и потенциал скорости р удовлетворяют уравнению движения, записанному в форме интеграла Коши-Лагранжа g + iVH2 + 5y + = c(t), (1.8) где д - ускорение свободного падения, c(t) - функция от времени, равная значению левой части уравнения (1.8) в некоторой точке пространства. Если на бесконечности жидкость покоится и давление на уровне у = 0 равно нулю, то ф) = 0. Если применить соотношение (1.4) для функции ip, то уравнение (1.8) перепишется в виде d -P 1 ,„ ,9 Р (1.9)
Постановка нестационарной задачи Для решения дифференциального уравнения (1.9) относительно ip и давления Р необходимо задавать искомые функции в некоторый момент времени (начальные условия) и удовлетворять определенным условиям на границе области (граничные условия). Граничные условия задаются либо в виде условия Дирихле (потенциал скорости принимает заданные значения на границе области), либо в виде условия Неймана (задается нормальная производная от потенциала скорости).
Приведем пример постановки нестационарной задачи движения жидкости в бассейне конечной глубины. Рассмотрим движение идеальной однородной несжимаемой жидкости под действием силы тяжести в области D, ограниченной свободной поверхностью Т\ и твердой границей Г 2 (рис. 1.1.1). В этом случае потенциал поля скоростей, зависящий от времени р(х, t), удовлетворяет уравнению Лапласа вида: Др(х, ) = 0, xe_D. (1.10) На твердых границах выполняются условия непротекания: п 0, хєГ2. (1.11) on дір дп На свободной поверхности должны выполняться кинематическое cJx. - = V , хєГь (1.12) и динамическое условия -\Ш2 + 9У + - = ф), хГ:. (1.13) at 2 р
Здесь с() - функция от времени, равная значению левой части в некоторой точке пространства. Если жидкость на бесконечности покоится и давление на уровне у = 0 равно нулю, то с() — 0. Из решения нелинейной стационарной задачи в качестве начальных условий задается начальное положение свободной границы Гі в момент времени t = 0 и распределение потенциала (р на ней. Постановка нелинейной стационарной задачи и ее решение рассмотрены в работах К.Е. Афанасьева, СВ. Стуколова [3,18,81].
Для удобства численной реализации задача приводится к безразмерному виду. В качестве характерных размерных величин выбираются ускорение свободного падения g и глубина бассейна Я. Безразмерные переменные х, у, t и С, записываются следующим образом: x = x/Ht у = у/Н, t = tM, C = = = F, (1.14) где С - скорость волны и F - число Фруда. При этом краевая задача (1.10)-(1.12) остается без изменений, а уравнение (1.13) принимает вид: -\\М + У = Ь хеГь (1.15)
Тильда над безразмерными переменными в формуле (1.15) и далее опускается. Требуется определить положение свободной поверхности и распределение поля скоростей на ней в последующие моменты времени. Компоненты вектора скорости определяются по формулам:
Нестационарное движение уединенной волны по бассейну с ровным дном
Данный тест представляет собой взаимодействие уединенной волны амплитуды А = 0.2; 0.3; 0.4; 0.5; 0.6 с вертикальной стенкой имеющей абсциссу х — 5. Для расчета была взята область D = {-15 х 5; — 1 у г/о}. На границе области было задано 552 узла, из них 251 на свободной поверхности. Остальные параметры были взяты из предшествующего теста. Численные результаты сравнивались с результатами полученными СВ. Стукоповым с помощью комплексного метода граничных элементов (КМГЭ) [18].
На рисунке 2,1.2а приведены профили свободной поверхности при накате солитона амплитудьгА — 0.4 на вертикальную стенку для нескольких моментов безразмерного времени, При і = 0.0, утах = 0.4 - первоначальная форма солитона; t — 8.91, утах = 0.946 и угэ = 0.945 - форма свободной поверхности в момент максимального заплеска; t = 19.12, утах — 0.393 и г/гэ = 0.392 - амплитуда восстановленного солитона, вершина которого находится в точке с абсциссой х = —7.5.
На рисунке 2.1.26 приведены профили свободной поверхности при накате солитона амплитуды А = 0.6 на вертикальную стенку для нескольких моментов безразмерного времени. При і — 0.0, утах — 0.6 - первоначальная форма солитона; t = 8.72, утах - 1.748 и 1гэ = 1.747 - форма свободной поверхности в момент максимального заплеска; t = 18.98, утах = 0.547 и угэ = 0.546 - амплитуда восстановленного солитона, вершина которого находится в точке с абсциссой х = -8.
Анализ численных результатов позволяет утверждать, что отражение волны от вертикальной стенки приводит к изменению амплитуды волны и формированию хвоста из вторичных волн малой амплитуды, что находится в полном соответствии с результатами численных исследований [72],
На рисунке 2,1.3(2,6 представлены зависимости кинетической, потенциальной и полной энергий в зависимости от времени. Видно, что максимум по тенциальной энергии и максимальный заплеск на стенку смещены по времени друг относительно друга (рис. 2.1.36), где показано изменение энергий при накате относительно высокой волны А = 0.6. Эта временная несогласованность может служить объяснением фазового сдвига набегающей и отраженной волн. В таблице 2.1.3 приведено сравнение с результатами расчетов из работы Б.Е. Протопопова [72]; J.W. Kim, K.J. Bai, R.C. Ertekin, W.C. Webster [117]. Теоретические результаты получены по формуле [130] 1 3 тахф, t) = 2А + -А2 + -А3. (2.1) 2 8
Для расчета была выбрана область D = {—35 х 35; — 1 у г/о}. В начальный момент времени t — 0.0 вершина волны находилась в точке х = 0. На границе области было задано 540 узлов, из них 301 на свободной поверхности. Из приведенных результатов видно совпадение результатов во втором и третьем знаке. Когда амплитуда волны А 0.4 значения для метода граничных элементов возрастают, но разница с теоретическими результатами находится в пределах 1%. Когда амплитуда волны А 0.4, то разница равна 10%. Это видно для амлитуды А = 0.6. Сравнение результатов с численными результатами других авторов показывает, достаточно высокую точность метода граничных элементов.
Кроме кинематических особенностей возникающих течений интерес представляет задача определения динамического воздействия волн на вертикальную стенку. На рисунке 2.1.4 для волн различных амплитуд показаны хроно граммы волнового давления в точке стенки, совпадающей с начальным урезом жидкости. Кривые на рисунке 2.1.4 обозначенные цифрами 1,2,3,4,5 соответствуют амплитуде волны А = 0.2, А = 0.3, А — 0.4, А — 0.5, А — 0.6 соответственно.
Схемы движения оползня
При моделировании рассматривают два вида оползней, называемые слайдами и слампами [133]. Оползень вида (слайд) можно определить как тонкий слой грунта. В работах [132,133] описан следующий закон движения слайда s(t) — SQIU cosh м. (3.7) где характерные время и расстояние вычисляются следующим образом t0 = Ui/aQ, s0 = u?/ao, t0 = щ/а0.
При данных допущениях уравнение (3.8) обеспечивает хорошую аппроксимацию начального разгона тела, поле чего движение становиться равномерным и продолжается до окончания расчета или остановки оползня. aot2 (3.8) 8(t) =
При естественных физических условиях начальное ускорение и скорость можно определить в виде a,Q = 0.3gsin щ = 1.16 ybgsm (/3). Второй вид оползня (сламп) продвигаются на угол ф по дуге окружности радиуса R. Полагая значение ф малым, a R относительно большим, эту дугу вследствие ее малой кривизны можно приблизить отрезком прямой. Тогда, закон движения можно аппроксимировать формулой (3.9), где SQ = О.ЬЯф, to — 1.8Ay/R/g, время движения 7го s(t) = s0 — COS t (3.9)
Для изучения волновых процессов в работе рассматриваются различные законы движения оползня (рис. 3.1.2) обсуждаемые в работе [91] : 1. слайд 1 (сді) - разгон, равномерное движение, остановка, покой; 2. слайд 2 (сд2) - разгон, равномерное движение, торможение, покой; 3. слайд 3 (сдЗ) - разгон, равномерное движение; 4. сламп 1 (cnl) - разгон, торможение, покой; 5. сламп 2 (сп2) - разгон, торможение, покой. 2. равномерное движение по закону xc(t) = xco+u(ti) + - до момента времени tCTOn; 3. остановка по закону xc(t) — XCQ -\- u(tQTOa - ti) 4- р- Тип движения "слайд 2" состоит из этапов: 1. разгона xc(t) = XCQ + Щ- до момента времени х; 2. равномерного движения xt:(t) = XCQ -\- u(t — t\) + - - до t%, 3. торможения xc{t) xCQ + u(t - t{) + Y— 2 до (стоп; 4. окончательной остановки xc(t) — жсо + u(iCTOri - ) + — 2&г Тип движения "слайд 3" состоит из этапов: 1. разгон xc{t) жсо 4- Щ- до ii; 2. равномерное движение жс() = х($+и( — i)+ ;р- до окончания расчета конед Для типа движения "слайд" приняты следующие значения параметров: а — 0.3sin( ), и = 1.16л/зш( Тип движения "сламп 1" включает следующие шаги: . 1. разгон xc(t) = Х( + SQ{1 — cos j-) до момента времени f.CTOn = тг ts\\ 2. покой xc(t) = XCQ 4- 2SQ.
Для исследования волновых режимов, порожденных различными законами движения оползня (рис. 3.1.2), была проведена серия вычислительных экспериментов в полной нелинейной постановке. При проведении расчетов выполнялось сравнение с результатами полученными по нелинейно-дисперсионным (НЛД) моделям: Мея-Меоте, Перегрина, Грина-Нагди, одно-двухслойной модели Лью-Линета и модели теории мелкой воды. Численные результаты для иерархии моделей волновой гидродинамики были предоставлены сотрудниками Института вычислительных технологий СО РАН [46,90,91,93].
Расчеты проводились в области с координатами XQ — 1,0 и хп = 41.0. Расчет проводился до времени і 50. Для изучения волновой картины были установлены семь мареографов с координатами: хмо = XCQ — 1, хм\ = &c0 = 2.38, %МІ ХМІ-1 + 2, і — 2..6. Для моделирования движения оползня используются следующие значения параметров: Ah — 0.05, Ь = 1.0, Я = 2.3, хс = 2.38, /3-6.
На рисунках (3.2,1)-(3.2.2) приведены графики мареограмм, рассчитанные для первого и седьмого мареографа. Сравнение выполнялось для закона движения оползня типа "слайд 1" и "сламп 1" соответственно. Результаты приводятся для первого (рис. 3.2.1,3.2.2а), второго (рис. 3.2.1,3.2.26) и седьмого (рис. 3.2.1,3.2.2е) мареографа. Первый мареограф фиксирует волну, распространяющуюся к берегу, а седьмой в открытую часть расчетной области D. Из рисунков видно, для закона движения "слайд 1" первый мареограф фиксирует волну в виде впадины, что характеризует начало движения оползня. Затем следует пологая волна повышения, связанная с этапом равномерного движения оползня. Далее, формируются волны порожденные остановкой оползня и отражением волн от левой береговой стенки (рис. 3.2.1д). Эти эффекты подтверждаются мареограммами второго мареографа (рис. 3.2.16), в котором "волна остановки" и "волна отражения" полностью разделены.
Схема последовательного алгоритма метода граничных элементов и его распараллеливание
В основе разработанного параллельного алгоритма для многопроцессорных систем с распределенной памятью лежит изложенный в первой главе последовательный алгоритм метода граничных элементов.
Для нестационарных задач со свободными границами характерен следующий алгоритм расчета, блок-схема которого приведена на рисунке 4.1.1.
В начальный момент времени известна геометрия границы расчетной области и граничные условия. Далее на каждом шаге по времени производиться вычисление интегралов, подстановка граничных условий, формирование системы линейных алгебраических уравнений (СЛАУ) и решение СЛАУ методом Гаусса с выбором ведущего элемента. Вычисляется шаг по времени, находится новое положение свободной границы и находится новое распределение значений граничных условий. В ходе расчета результаты записываются в файл для дальнейшей обработки или возобновления расчета.
Ввод данных
Для выполнения расчета необходимы следующие граничные и начальные данные: 1- #і) Уй і = 1)п " узловые координаты точек на границе области; 2. Code i — 1,п - тип граничных условий в г-ом узле, если Codei = О, то в узле задана функция потенциала р, если Codei = 1 - функция нормальной производной др/дп; 3. Ui,Bi,i — 1,п - заданное значение функции (р и дір/дп соответственно в г-ом узле.
Распараллеливание этого блока не требуется, так как ввод всех необходимых начальных и граничных данных происходит во время расчета лишь один раз. Для избежания межпроцессорных пересылок считывание данных целесообразно осуществить на каждом процессоре вычислительного комплекса.
Вычисление интегралов и решение СЛАУ
На этом шаге алгоритма формируется СЛАУ матрица А размерностью N х N я вектор правой части F (1.24). Следует привести некоторые соображения: каждая строка результирующей матрицы А строится независимо от других, следовательно, степень параллелизма данного блока алгоритма максимальна и равна п. Далее, на этом же этапе следует позаботиться о цикличном слоистом распределении колонок матрицы А по процессорам для ее дальнейшего решения.
Так как во втором блоке алгоритма распределение матрицы А по процессорам было колоночное, то для LU - факторизации целесообразно использовать колоночно-ориентированный алгоритм, приведенный далее.
Построение результирующей функции
Если алгоритм решения получившейся в результате LU - факторизации системы уравнений распараллелен, то результат решения распределен по процессорам, если нет - то он собран на одном процессоре. Обозначим вектор решения через СІ,І = 1, п. Суть построения результирующей функции заключается в следующем: если в г-и узле задан потенциал, то есть Code,- = 0, то значение вектора решения С\ равно результирующей функции В{ — С,-; если в г-м узле задана функция нормальной производной, то есть Codei — 1, то значение вектора решения d равно результирующей функции Щ = Q;
