Содержание к диссертации
Введение
1. Методы пространственной аппроксимации задач волнопой гидродинамики 18
1.1. Модель малкой воды 18
1.2. Краевые условия 20
1.3. Методы численного решения уравнений мелкой воды 22
1.4. Пространственная аппроксимация 23
2. Модифицированный метод конечных объемов 27
2.1. Основные принципы построения МКС) 27
2.2. ММКО для задач волновой гидродинамики 30
2.3. Интегральные формулы МКО 31
3. ММКО-аппроксимация волнового уравнения 33
3.1. Постановка задачи 33
3.2. Дискретизация по времени 33
3.3. Интегральная форма законов сохранения 31
3.4. ММКО на линейных базисных функциях Зо
3.4.1. Способ построения двойственной сетки 35
3.4.2. Аппроксимация волнового уравнения 3fi
3.4.3. Учет краевых условий 38
3.4.4. Сборка глобальных матриц
3.5. ММКО на квадратичных базисных функциях
3.5.1. Особенности построения двойственной сетки .
3.5.2. Аппроксимация волнового уравнения
3.5.3. Учет краевых условий
3.5.4. Сборка глобальных матриц
4. ММКО-аппроксимация уравнений мелкой поды
4.1. Постановка задачи
4.2. Дискретизация по времени
4.3. Способ построения двойственной сетки
4.4. Интегральная форма уравнений мелкой воды
4.5. Аппроксимация уравнения неразрывности
4.6. Аппроксимация уравнений движения
4.7. Учет краевых условий
4.8. Сборка глобальных матриц
5. Описание программного комплекса
5.1. Основные модули пакетов
5.2. Организация входных и выходных данных пакетов .
6. Численные эксперименты
6.1. Моделирование плоской волны в канале
6.2. Проверка адекватности "открытых" краекых условий .
6.3. Взаимодействие уединённой волны с коническим островом
6.4. Моделирование цунами 1940 года в Японском море .
Заключение
Список литературы
- Методы численного решения уравнений мелкой воды
- Интегральная форма законов сохранения
- Дискретизация по времени
- Организация входных и выходных данных пакетов
Введение к работе
В последнее время численное моделирование играет всё более значительную роль в исследовании реальных явлений. Совместное проведение вычислительных и физических экспериментов при анализе какого-либо явления позволяет как уменьшить количество реальных измерений, так и произвести верификацию и усовершенствование математических моделей. Кроме того, существуют такие глобальные задачи, которые в силу очевидных причин невозможно моделировать экспериментальным образом. Одной из таких важных задач является задача распространения волн цунами в водных бассейнах, имеющих сложную структуру береговой линии и распределения глубин.
Поэтому используемые для моделирования такой задачи вычислительные методы должны предоставлять возможность наиболее точного описания геометрии расчётной области. Это возможно при использовании неортогональных и неструктурированных сеток. Использование неструктурированных сеток позволяет описывать с любой степенью точности многосвязную расчётную область с произвольной конфигурацией границы, а также даёт возможность реализовать локальные сгущения и адаптировать сетки в зависимости от поведения решения, либо распределения глубин. Создание алгоритмов численного моделирования для данного класса задач, использующих неструктурированные сетки, представляется АКТУАЛЬНОЙ ТЕМОЙ ИССЛЕДОВАНИЯ.
Диссертационная работа посвящена разработке технологий метода конечных объёмов (МКО) на неструктурированных сетках для задач волновой гидродинамики, описываемых в рамках теории мелкой воды, и созданию комплексов программ для проведения численных
экспериментов в рамках разработанных технологий. В данной рабо-
J*GC. НАЦИОНАЛЬНАЯ)
„ БИБЛИОТЕКА І
те рассматривается класс модифицированных методов конечных объёмов (ММКО) на финитных базисных функциях первого и второго порядков. Расчёт неизвестных производится в узлах конечноэлемент-ной сетки, т.е. в узлах триангуляции для кусочно-линейных базисных функций и в узлах триангуляции и центрах рёбер для кусочно-квадратичных базисных функций.
Итак, ЦЕЛЬЮ НАСТОЯЩЕЙ РАБОТЫ является разработка и применение технологии ММКО построения дискретных аналогов задач волновой гидродинамики. Для достижения заданной цели были сформулированы следующие задачи исследования:
Разработка новой технологии ММКО-аппроксимации для уравнений в частных производных первого и второго порядка с использованием кусочно-квадратичных базисных функций;
Разработка технологии смешанной ММКО-аппроксимации для решения системы уравнений нелинейной теории мелкой воды;
Создание на основе разработанных технологий комплексов программ, позволяющих адекватно моделировать распространение длинных волн в водном бассейне с геометрически сложными границами.
МЕТОДЫ ИССЛЕДОВАНИЯ. Методы вычислительной математики. Сравнительный анализ результатов моделирования и имеющегося аналитического решения. Расчёты на последовательности сгущающихся конечноэлементных разбиений, с последующим анализом сходимости к экспериментальным данным.
НАУЧНАЯ НОВИЗНА работы состоит в следующем:
Предложена технология учета кусочно-полиномиального представления решения и коэффициентов при пространственной аппроксимации линейной модели мелкой воды, системы нелинейных уравнений мелкой воды методом конечных объёмов. Технология основана на использовании разложения по базису конеч-ноэлементных пространств в терминах барицентрических сим-плициальных координат, с дальнейшим точным интегрированием их одночленов. В рамках предложенной технологии впервые использован ММКО на квадратичных базисных функциях, для чего получены соответствующие интегральные формулы.
Предложен способ построения смешанного модифицированного метода конечных объёмов на совмещённых симплициальных сетках, удовлетворяющего условиям Ладыженской-Бабушки-Брец-ци (ЛББ).
С использованием предложенных технологий аппроксимаций задач волновой гидродинамики в рамках линейной и нелинейной теории мелкой воды созданы комплексы программ для моделирования процессов распространения волн и проведён ряд вычислительных экспериментов, подтверждающих эффективность построенных схем.
СТРУКТУРА И ОБЪЁМ ДИССЕРТАЦИИ. Диссертационная работа состоит из введения, шести глав, заключения, списка литературы, приложений и содержит 777 страниц, включая 41 рисунок и 2 таблицы. Список литературы содержит 65 наименовании.
Методы численного решения уравнений мелкой воды
Адекватность численного решения модели теории мелкой воды зависит от многих вычислительных аспектов. Система уравнений мелкой воды - это система нелинейно связанных дифференциальных уравнений п частных производных, определенных в расчётной области со сложной геометрией границ, с заданными на них граничными условиями, и, возможно, существенно сложным распределением глубин. Кроме того, в эту систему может быть включен учет внешних сил. таких как силы Кориолиса, донного трения, ветрового напряжения [49, 51, 35]. В дополнение к этим физическим факторам, усложняющим алгоритмы моделирования, существуют другие трудности, имеющие математическую природу. Наибапее важная -это связь между полной глубиной и нолем скоростей, которая приводит к возникновению паразитных осцилляции при реализации численных алгоритмов решения уравнений мелкой воды, пс учитывающих эти нюансы специальным образом.
Почти все разработанные методы для численного решения системы уравнений мелкой воды можно условно разделить па дна широких класса, В первом случае производится моделирование непосредственно системы уравнений мелкой поды (1.1)-(1.2). В литературе встречается множестко подходов для численного интегрирования этой модели с использованием техники конечных разностей, конечных члемеитпв. Для устранения нефизичиих пространств em п.тх оеннлляпнй решения, возникающих вследствие нелинейной связи полной глубины и поля скоростей, первоначально использовались такие подходы, как введение п модель мелкой воды исскуственной вязкости [45], использования численного сглаживания решения [20, 64] и введения в модель чрезмерного донного трения 21, Однако более эффективным способом минимизации паразитного шума решения стало использование разнесенных шаблонов для конечных разностей или выбор интерполяционных подпространств разных порядков в случае конечных элементов [44, 12].
Во втором классе методов решения системы уравнений мелкой воды, основоположниками которого были Линч и Грэй, выполняется переформулировка исходной задачи таким образом, что вместо уравнения неразрывности первого порядка (1.1) используется гиперболическое уравнение второго порядка, называемое волновым уравнением мелкой воды (SWWE) [49], решаемое совместно с стандартным уравнением движения (1.2). Численное решение задачи в такой постановке свободно от возникновения кефизичных пространствештых осцилляции, без использования специальных методов.
Традиционно, наиболее популярными методами пространственной аппроксимации для решения системы уравнений мелкой воды являются методы конечных разностей и конечных элементов. Обзор методов, базирующихся на конечно-разностной пространственной аппроксимации, достаточно хорошо представлен в работе [40]. Однако, несмотря на широкую распространённость технологии конечных разностей, использование её для модели мелкой воды в ряде случаев не способно дать необходимую точность численного решения. Как правило, вычисления необходимо реализовать в областях, возможно многосвязных, с существенно сложной геометрией границ акватории, на которых обычно заданы краевые условия полного отражения (1.6) п "открытые" граничные условия (1.7). Поэтому, для получения наиболее точного решения в прибрежных районах необходимо как можно более точное описание береговой линии, что практически невозможно реализовать па ортогональных сетках, используемых в МКР.
Наиболее распространённым способом описания областей со сложной геометрией является использование неструктурированных симилициальных сеток, позволяющих эффективно сгущать либо разрежать узлы сетки при ре тисни и задач с большими градиентами решений [50, 48]. Однако, несмотря на очевидные достоинства неструктурированных сеток, их применение сдерживалось трудоёмкостью расчётов, которая резко возрастала с переходом к ним от структурированных сеток. Это объясняется тем, что положение и количество ненулевых элементов в матрицах дискретных аналогов зависит от смежности узлов сетки, поэтому матрицы приходится хранить с использованием специальных форматов [5, 57]. Все это делает гораздо более трудоёмкими основные операции умножения матрицы па вектор, что существенно замедляет итерационные схемы решения уравнений вычислительной гидродинамики. Разработка и усовершенствование процедур автоматического построения неструктурированных сеток [2, 47, 15], а также появление мощных вычислительных систем, дали толчок к развитию численных методов, базирующихся иа таких сетках. Для применения пространственной дискретизации в этом случае используется метод конечных элементен, и в последнее время - метод конечных объёмов. Например, в работе Линча и Грэя [49] МКЭ был использован для пространственной аппроксимации волнового уравнения мелкой воды и уравнения движения, для временной дискретизации была использована конечнорлл-постная схема. Разновидности схем шіда предиктор-корректор в рамках коисчпо-элементой пространственной аппроксимации использованы некоторыми авторами [25, 411. В работах Перейры, Зинкевича, Моргана [51] и Амброси, Квартапелли [16] для численного реигения системы уравнений мелкой ноды был использован метод Тейлора -Галеркина. В рамках модели мелкой воды Гиральдо использовал метод Лагранжа- Галеркина [35].
В работах Ладыженской, Бабушки и Темами [3, 17, 12] показано, что для устойчивых- смешанных МКЭ-аппроксиманий уравнений Навье Стокса в естественных переменных, выбранные пространства для представления давления и скоростей должны удовлетворять условию ЛББ. Дня выполнения этого условия, Тейлор и Худ [62] доказали необходимость представления поля скоростей функциями формы на единицу большего порядка, чем функции формы для давления. Следовательно, при реализации конечнозлементно-нодобшлх методов для системы уравнений мелкой воды выполнение этих условий для аппроксимации функции свободной поверхности, аналога давления из уравнений Навье Стокса, и поля скоростей будет приводить к устранению нефизичных пространственных осцилляции. Применение алгоритмов типа SIMPLE [53] позволяет достичь аналогичного результатів с использованием функций формы одного порядка.
Интегральная форма законов сохранения
Предположим, что в двумерной области S2 задана некоторая триангуляция, такая, что два ее симплекса (конечных элемента), могут иметь общими лишь узел или ребро. При решении уравнения (3.6) будем считать, что значение волновой поверхности рассчитывается точках, которые совпадают с вершинами конечных элементов (рис. 3, а). Тогда двойственная, по отношению к конечно-элементной, сетка, для аппроксимации уравнения (3.1) на каждом симплексе, может быть образована центроидами треугольников и центрами их сторон [58]. Такт: образом, произвольная вершина триангуляции г на конечном элементе окружена отрезками медиан (рис. 3, а), образующими часть конечного объёма П;, соответствующего этому узлу,
Аппроксимация волнового уравнения
Введем локальную нумерацию узлов на симплексе Iі (рис. 3, а). Вазовые конечные объёмы І1 введем так, что ІІІ = Ullj. Пусть 5; граница конечного объёма з fi;, ассоциированная с г -ым узлом, тогда 5;;- часть границы S,, образованная отрезком медианы между центроидом треугольника " стороной, соединяющую узлы г и j, имеющая направление, соответствующее обходу Ї2; против часовой стрелки; J І значение суммарного потока через составную границу S{.
Локальная нумерация узлов и обозначения на конечном элементе, (6) -- Обозначения диффузионные потоков, рассчитываемых через границы КО.
Дія построения локальной ММКО-аппроксимации волнового уравнения на симплексе Т достаточно построить дискретные аналоги интегралов вида / rjdil и (і. / h{i)xdy - Tjydx), і = 1, 2, 3 на каждой части конечного объёма. ж, Построепие локальной матрицы массы
Так как для представления функции возвышения свободной поверхности на Т используется кусочно-линейная интерполяция, т.е. V HviLi, (3.8) где f}i - значення возвышения свободной поверхности в точке і, Li базисные функции вида (2.1), введенные на конечном элементе Тк, соответствующие г-ому узлу, дискретный аналог интеграла / rjdQ, і — 1, 2, 3, на симплексе Т , можно (і представить в виде: = 1,2,3, (3.9) Mir)1 - Щ \ f Lj&l + / LJdQ J=l U nj где M локальная матрица массы, имеющая размерность 3x3, г) - гипорвек-тор, составленный из компонент возвышения свободной поверхности, заданных в соответствующих им узлах симплекса. Для определения элементов этой матрицы необходимо рассчитать значения интегралов в фигурных скобках выражения (3.9), например с помощью формул точного интегрировании (2.2), с учіп-им преобразования (2.4). Рассчитанные таким образом значения интегралов приведены в ПРИЛОЖЕНИИ А.
Построение локальной матрицы жесткости
Обозначим J{ - суммарный диффузионный поток в направлении внешней нормали через часть границы конечного объёма її;, принадлежащей симплексу Т . Очевидно, что достаточно рассмотреть построение аппроксимаций трех определяющих потоков (например, Jv2, hz) -hi) ич возникающих шести, поскольку J\ — —Ji2, hi = —J-ІЗІ J\n = —-Тзі (3, b). Исходя из этих соображений, можно записать: J\ = J\2 — ./зі» h = Лз - 12, h — -h\ — - 23- (3.10)
Для интерполяции функции возвышения свободной поверхности и» конечном элементе Т и представления глубины Л, используются кусочно-лшкійішо базис 38 ные функции нила (2.1), следовательно определяющие потоки имеют вид [13, 14]: з_з . hi = Yl h Y,fli LAdy LjCidx, (3.11) i=» -i si где к = 1, 2, З, I = 1, 2, 3, 6,-,с,- - коэффициенты линейной комбинации (2.1).
Зная выражения для определяющих диффузионных потоков па конечном а -менте, можно записать локальную матрицу жесткости J симплекса Tfc, имеющую размерность 3x3: ( J ЛтГ = (3.12) Интегралы в (3.12) рассчитываются аналитически по формулам, подобным (2.3) согласно преобразованию (2.4). Значения этих интегралов приведены н ПРИЛОЖЕНИИ А.
Дискретизация по времени
Для построения дискретизации по времени применим технику, подобную использованной для волнового уравнения. Для этого представим разложение величины р в ряд Тейлора по т, по временным слоям вплоть до второго порядка: 7}k+l + rVk + ч + 0{т\ (4.6) где т] = r)(t ). Испатьзуя для аппроксимации производной второго порядка следующую конечноразностную аппроксимацию: где { т, если 0 = 1, 1, если 9 ф1у выражение (4.6), можно переписать в следующем виде: +I - тгі+І = 9 + r(l - Ы 4- О(т о). (4-7)
При параметре в = 1 выражение (4.7) имеет второй порядок точности, для других значений в - только первый порядок точности. При значении 9 = 0 предложенная схема соответствует явной, полностью неявная схема может быть получена при значении 0 = 2.
Подставляя в (4.7) значения (ри выраженные из (4.5), можно записать схему временной дискретизации для нелинейной модели мелкой воды: pk+l + r[v ( pk+lVrk+i) + дИк+1Р (Нк+Х)] = рк + \тдНкР {h) + + ( - l) r[V ( pkVTk) + дНкР (Нк) - дНкР (Л)]. Линеаризуя полученное выражение, для чего члены Urfc+1 и Нк+1 замеяя значениями U и Я с предыдущего тага по времени [11], можно получить: pk+l + 2r[V ( +1 U7 ) + уНкР (tf +l)] = рк + т9НкР (А) + + ( - і) (y U) + дНкР (ff ) - дНкР {h)\. В результате была получена временная дискретизация системы уравнений мелкой воды (4.5), соответствующая по идеологии временной дискретизации волнового уравнения (3.1).
Предположим, что в двумерной области U задана некоторая триангуляция, такая, что два ее симплекса (конечных атемента), могут иметь общими лишь узел или ребро. При решении системы (4.8) будем считать, что поле скоростей и значение полной глубины рассчитываются в одних и тех же узлах и узлы эти совпадают с вершинами конечных элементов (9, а), за исключением дополнительных узлоп серединах рёбер симплекса, соответствующих квадратичным функциям формы для аппроксимации поля скоростей (рис. 9, Ь). Тогда двойственная, но отношению к конечноэлементной, сетка, для аппроксимации уравнения неразрывности на каждом симплексе, может быть образована центроидами треугольников и центрами их сторон [58]. Таким образом, произвольная вершина триангуляции і па конечном элементе окружена отрезками медиан (рис. 9, а), образующими часть конечного объёма Ї7(, соответствующего этому узлу. В случае аппроксимации уравнений движения, двойственная сетка на каждом конечном элементе может быть получена путём использования медиан четырех малых треугольников, образованных вершинами симплекса и центрами его сторон, так что каждый из узлов имеет конечный объём, подобный объёмам, соответствующим кусочно-линейной интерполяции (рис. 9, Ь) [58]. Такая конструкция была впервые предложена Пракашем и Патапкаром в [53) в качестве аналога структурированных разнесённых сеток.
Таким образом, для і-го узла триангуляции, ММКО-дискретизапия заключается в построении аппроксимаций конвективных потоков через совокупность граней дІ\\ соответствующих узлу конечного объёма Qj, а также аппроксимации объёмных интегралов на этом конечном объёме. В результате представления решения и коэффициентов системы уравнений на конечных элементах кусочно-полиномиальными интерполянтами с конечным носителем, из (4.9) получаем дискретный аналог интегро-балансного соотношения, которое связывает значение искомой функции в і-ом узле со значениями в смежных с ним узлах триангуляции. Следовательно, разработанная технология ММКО использует информацию только о базовой (первичной) триангуляции.
Подобно технологии метода конечных элементов, генерация матриц, соответствующих дискретному аналогу, выполняется обходом не по узлам триангуляции, а по конечным элементам [58}, что позволяет избежать многократных повторных вычислений геометрических характеристик на каждом из симплексов,
Введем локальную нумерацию узлов па симплексе Т (рис. 10, а). Базовые конечный обт ёмы ІЦ введем так, что П,- = Ufi;. Пусть .9,- граница конечного объёма Qi; з ассоциированная с г-ым узлом, тогда Sy - часть границы 5,-, образованная отрезком медианы между барицетром и стороной, соединяющей узлы г и j, имеющая направление, соответствующее обходу Q; против часовой стрелки; .7; значение суммарного потока через составную границу S .
Организация входных и выходных данных пакетов
Этот тестовый расчет предназначен для иллюстрации работоспособности предложенных вычислительных методов в расчётной области, геометрическими характеристиками приближенной к реальной. В качестве расчётной области О. была взята реальная акватория Японского моря с достаточно детальным описанием геометрии границ и распределением глубин.
На этапе подготовки данных к вычислительному эксперименту были сделаны следующие упрощения расчётной области. Так как используемые модели волновой гидродинамики не предназначены для описания наката волн, в прибрежной зоне глубины, меньшие значения минимальной глубины Л0 заменялись этим значением.
На рисунке 32 показана часть Японского моря, в котором производился вычислительный эксперимент. Слева сверху располагается побережье России, снизу справа - побережье о. Хонсю, справа - побережье о. Хоккайдо; все эти побережья моделировались с помощью граничных условий вида "жесткая стенка". В расчётной области заданы о. Окушири и о. Ребун, а также часть о. Садо. В левом нижнем и в правом верхнем углу расчётной области П были заданы "открытые" краевые условия, позволяющие волне покидать расчётный регион без отражения. На рисунке 33 приведено модифицированное распределение глубин в расчётной области по описанной выше методике. Выбор минимальной глубины ho = 75 м. обеспечивал устойчивый расчёт по классической модели мелкой воды смешанным ММКО. Начальное возмущение, описывающее очаг цунами 1940 года в Японском море, задавалось в виде верхней половины эллипса: f (х - Jo)2 (У - Уо)2 где хо = 902898 м., у0 = 137196 м., г = 37080 м., Я = 74160 м.( т?0 = 4 м.
В расчётной области Q была построена треугольная сетка, состоящая из NP = 13624 узлов и NT = 26191 треугольников, характерный размер грани которых !)() был задан 9270 метров. Часть этой триангуляции приведена на рис. 34.
На рисунках 35 41 показано развитие волнового процесса в Японском море с течением времени. В силу иллюстративности этого геста на этих рисунках приведены только расчёты по классической модели мелкой воды по вычислительной схеме сметанного ММКО. Численное моделирование этой задачи производилось с параметром схемы Тейлора Галёркина в равным 1, временной шаг т был использован равным 0.625 с.
На основании результатов этого численного эксперимента можно сформулировать следующие выводы:
1. Предложенная технология расчёта на неструктурированных сетках позволяет производить вычислительный эксперимент в областях с реальными очертаниями границ и распределением глубин.
2. Предложенная реализация граничных условий работоспособна при задании условий на геометрически сложных границах.
Настоящая работа посвящена созданию вычислительных технологий, соответствующих МЛІКО-аппроксимациям задач волновой гидродинамики, описываемых вол-нововым уравнением и системой уравнений нелинейной теории мелкой воды. Характерной особенностью разработанных алгоритмов является использование неструктурированных симшшциалъных разбиений расчётной области, а также использование конечиозлементных пространств и барицентрических разбиений в качестве двойственных. В диссертационной работе получены следующие основные результаты, выносимые на защиту:
1. Предложена технология кусочно-полиномиального представления решения задач волновой гидродинамики методами конечных объёмов. Технология основана па разложения по базису конечноэлементных пространств в терминах барицентрических координат, с дальнейшим тонным интегрированием их одночленов.
2. Разработан и реализован вычислительный алгоритм на базе модифицированного метода конечных объёмов с использонанием локальных базисных функций второго порядка, что позволяет удовлетворить условиям ЛВБ для системы уравнений нелинейной теории мелкой воды.
3. Разработан и реализован алгоритм для решения нестационарной задачи волновой гидродинамики с использованием технологии метода Тейлора Галёркина с параметром, управляющим вносимой диссипативностью схемы.
4. На основании разработанных вычислительных схем созданы комплексы программ для моделирования распространения волн для волнового уравнения на интерполяционных пшіиномах для свободной поверхности первого и второго порядков; и для системы уравнений мелкой воды на кусочно-линейных интерполяционных полиномах для полной глубины и кусочно-квадратичных интерполяционных полиномах для поля скоростей.
5. Выполнено тестирование разработанных вычислительных схем на задачах, имеющих аналитическое решение; проведены серии расчётов для задач, для которых имеются результаты моделирования других авторов.
