[ «Численное моделирование течений в приближении мелкой воды на основе регуляризованных уравнений ...»
Московский государственный университет им. М.В. Ломоносова
Физический факультет
На правах рукописи
Булатов Олег Витальевич
Численное моделирование течений в приближении
мелкой воды на основе регуляризованных уравнений
Специальность 05.13.18 –
математическое моделирование,
численные методы и комплексы программ
Диссертация на соискание ученой степени кандидата физико-математических наук
Научный руководитель:
доктор физ.-мат. наук, профессор Елизарова Татьяна Геннадьевна Москва – Оглавление Page Введение.................................... Глава 1. Регуляризованные уравнения мелкой воды и метод численного решения.................................. 1.1 Уравнения мелкой воды...................... 1.2 Регуляризованные уравнения мелкой воды............ 1.3 Регуляризованные уравнения для течений при малых числах Фруда................................. 1.4 Численный метод для одномерных течений........... 1.5 Задача Римана............................ 1.5.1 Построение автомодельного решения........... 1.5.2 Результаты численного расчета.............. 1.6 Задача о транскритическом течении над неровностями дна.. 1.7 Задача об отражении поверхностных волн от подводной возвышенности............................... 1.8 Автомодельное решение и численное моделирование задач Римана при наличии уступов дна.................... 1.8.1 Введение........................... 1.8.2 Тест 1............................. 1.8.3 Тест 2............................. 1.8.4 Тест 3............................. 1.8.5 Тест 4............................. 1.8.6 Тест 5............................. 1.8.7 Оценка точности численного метода............ Глава 2. Условие сухого дна для одномерных задач............ 2.1 Постановка условия для сухого дна на примере водоема с холмом и сухим верхом............................ 2.2 Одномерный разрыв........................ 2.3 Задача Римана с разбегающейся жидкостью........... 2.4 Сравнение точного решения с численными расчетами для случая постоянного наклона дна...................... 2.5 Набегание цунами на наклонный берег.............. Глава 3. Обобщение алгоритма на пространственные течения...... 3.1 Численный алгоритм для двумерной прямоугольной сетки... 3.2 Условие покоящейся жидкости для неровного дна........ 3.3 Задача о разрушении несимметричной дамбы.......... 3.4 Набегание цунами на берег сложной формы........... 3.5 Расчеты волны прорыва в расширяющемся канале....... Глава 4. Численный метод для неструктурированных сеток....... 4.1 Разностная аппроксимация уравнений.............. 4.2 Эффективная реализация численного алгоритма........ 4.3 Условие покоящейся жидкости................... 4.4 Задача о разрушении столба жидкости.............. 4.5 Задача о разрушении плотины и затоплении поверхности с тремя конусами............................... Заключение................................... Литература................................... Приложение А. Усовершенствования алгоритма для задач сухого дна.. Приложение Б. Численное моделирование течений газа на основе квазигидродинамических уравнений.................... Б.1 Система уравнений и численный алгоритм............ Б.2 Задачи Римана о распаде разрыва................. Введение Движение несжимаемой жидкости со свободной поверхностью в поле сил тяжести может быть описано в приближении мелкой воды. Уравнения мелкой воды (МВ) представляют собой упрощенную модель полных уравнений Навье-Стокса, описывающих пространственные нестационарные течения вязкого сжимаемого газа [1], [2], [3], [4], [5]. При выводе системы уравнений МВ предполагается, что среда представляет собой достаточно тонкий слой, глубина которого много меньше его продольного размера, поэтому вертикальной составляющей скорости в слое можно пренебречь и полагать, что продольные скорости постоянны по толщине слоя. Дополнительно предполагается, что жидкость несжимаема, находится в поле сил тяжести, и ее температура постоянна [6], [7] [8].
Математическая модель мелкой воды широко используется для решения задач, представляющих как академический, так и практический интерес. К последним относится моделирование течений в относительно неглубоких водоемах, реках, водохранилищах, течений вблизи побережья морей и океанов, расчет волн цунами и сброса вод вблизи гидроэлектростанций, а также множество других задач, непосредственно связанных с проблемами экологии. Это приближение применяется для описания береговых течений, течений в реках и озерах, гидравлических течений в водозаборниках, технических сужениях и лотках, распространения волн прорыва при разрушении гидротехнических сооружений, для численного моделирования волн цунами и приливных бор в реках.
Подробную информацию о разнообразных задачах физики, где применяется модель мелкой воды, можно найти в книгах [8], [9]. Приближение мелкой воды применяется к атмосферным течениям и используются для задач прогноза погоды [10]. Уравнения мелкой воды используется при численном моделировании крупномасштабных атмосферных и океанических течений ( [11], [12], [13]), где существенны ускорения Кориолиса и его широтные вариации.
Если жидкость расслаивается по причине разной солености или температуры, то полученный в результате слоистый поток по своей структуре похож на течение мелкой воды. Примеры использования многослойной модели мелкой воды приведены, например, в [8], [14].
В последние десятилетия был разработан целый ряд численных алгоритмов для моделирования задач в приближении мелкой воды. Среди способов дискретного расчета уравнений мелкой воды получили наибольшее распространения три численных метода. К ним относятся метод конечных разностей ( [20], [21]), метод конечных элементов ( [22], [23], [24]) и метод конечного объема ( [25], [26], [27], [28], [29], [30], [31], [32], [33], [34]). В русскоязычной литературе, следуя терминологии А.А. Самарского [62], называют интегро-интерполяционным методом. Основным преимуществом метода конечного объема является его понятная физическая интерпретация, локальное и глобальное сохранение массы жидкости, а также простота, с которой метод конечного объема расширяется и обобщается для неструктурированных сеток.
Уравнения МВ выведены для случая, когда высота уровня жидкости много меньше характерных размеров задачи, и форма дна является достаточно гладкой функцией. Тем не менее задачи с разрывным профилем дна интенсивно изучаются в рамках приближения МВ как аналитически, так и численно. Построение аналитических решений для таких задач достаточно трудоемко даже для плоских одномерных течений, и этим вопросам посвящена обширная литература (см., например, [8], [15], [16], [17], [18]). Аналитические решения служат как для оценки ситуаций, возникающих в ряде практических случаев, так и для тестирования численных алгоритмов, разработанных для расчета таких течений. К подобным течениям относятся, в частности, течения на порогах шлюзов, течения при разрушении шлюзовых ворот или при переливах воды через гребень плотины (см. [19]), течения в узких морских проливах со сложной формой дна (например, течение в проливе Гибралтар) и ряд других задач.
Трудности при численном моделировании задач с разрывным профилем дна вызваны возникновением сложной конфигурации разрывов в решении, обусловленных как нелинейностью самих уравнений, так и разрывным профилем подстилающей поверхности. В [17], [18], [19] предложены способы преодоления этих проблем путем выделения линии разрыва, связанной с положением границы уступа или ступеньки дна и модификации системы уравнений МВ. Использование такого подхода делает численный алгоритм более точным, но лишает его однородности. Последнее не всегда удобно при расчетах практических задач. Другим способом решения поставленной задачи является использование двухслойных уравнений мелкой воды (см. [8]).
Таким образом, построение удобного однородного численного алгоритма для решения задач с разрывами дна представляется актуальным.
При численном моделировании течений жидкости со свободной поверхностью часто возникаю ситуации, когда высота уровня жидкости становится малой, то есть возникают так называемые зоны сухого дна. Например, такие ситуации возникают при расчетах течений рек, затоплении и осушении низменностей, набегании волн на береговую линию, расчет прибрежных волн цунами и в других случаях. Трудности в численном моделировании течения жидкости с областями сухого дна связаны с появлением движущейся границы, разделяющей сухую область и область, занятую жидкостью. Это довольно сложная проблема, которой посвящено множество работ ( [31], [32], [33], [36]).
В связи с перечисленными выше задачами усовершенствование и разработки новых эффективных алгоритмов для математического моделирования течений в приближении МВ является актуальной. Новые численные алгоритмы, удобные для практических применений, должны обладать следующими особенностями:
алгоритм должен быть достаточно универсальным и однородным для моделирования течений с неизвестными заранее особенностями, такими как гидравлические скачки и волны разрежения допускать возможность расчета течений с подвижными областями сухого дна адаптироваться к сложным неструктурированным расчетным сеткам, которые требуются для описания течений в сложных пространственных областях - например, в задачах затопления в поймах и руслах рек представлять возможность расчета течений в зонах со сложной формой подстилающей поверхности, включая ступеньки и уступы дна допускать возможность распраллеливания алгоритма на большое число процессоров для ускорения счета Данная диссертационная работа посвящена созданию, программной реализации и верификации нового численного алгоритма решения уравнений МВ, удовлетворяющего перечисленным свойствам.
По своей природе уравнения мелкой воды тесно связаны с уравнениями динамики газа. В отсутствие внешних сил и дополнительных усложняющих задачу факторов уравнения МВ формально можно выписать на основе уравнения Эйлера в баротропном приближении. Эта аналогия между уравнениями МВ и системой уравнений Эйлера для невязкого сжимаемого газа хорошо известна и является причиной того, что численные алгоритмы, используемые для решения уравненй МВ, как правило, основываются на методах, развитых для уравнений Эйлера.
В работах [37], [38] и [39] предложен способ построения сглаженных, или регуляризованных уравнений Навье-Стокса и Эйлера. На этом пути были выписаны квазигазодинамические и квазигидродинамические (КГД) уравнения, которые показали свою эффективность при численном моделировании широкого круга течений вязкого сжимаемого газа и несжимаемой жидкости.
КГД-уравнения отличаются от уравнений Навье-Стокса и Эйлера дополнительными дивергентными слагаемыми, появление которых связано с использованием сглаживания исходных уравнений по малому интервалу времени, что приводит к появлению дополнительных нелинейных слагаемых с малым параметром. Эти дополнительные -слагаемые выполняют роль регуляризаторов и обеспечивают устойчивость и точность численных алгоритмов, построенных на основе КГД подхода.
Эта идея легла в основу диссертационной работы, где впервые построены регуляризованные, или сглаженные, уравнения мелкой воды (РУМВ), на основе которых разработаны новые однородные и эффективные численные схемы для математического моделирования течений со свободной поверхностью.
В первой главе выписаны уравнения МВ, предложены два способа построения регуляризованных уравнений мелкой воды и показана их прямая связь с КГД уравнениями. Первый способ является более общим и применим к течениям с произвольным числом Фруда. Второй способ удобен для расчета течений с малыми скоростями. Численный метод строится для более универсального первого варианта регуляризованных уравнений (параграф 1.4).
Алгоритм явный, используется метод конечных объемов, потоковые величины аппроксимируются центральными разностями Устойчивость численного алгоритма обеспечивают дополнительные -слагаемыми. Шаг по времени и пространственный шаг связаны условием Кураната.
Важным свойством алгоритмов для численного моделирования течений в приближении МВ является выполнение условий покоящейся жидкости для течений над сложной формой подстилающей поверхности. Последнее означает, что в изначально покоящейся жидкости не должны возникать возмущения, обусловленные неровностями дна. В англоязычной литературе численный алгоритм, обладающий этим свойством, называют "well-balanced scheme"(см. [36] или [41], [42]). В построенном автором алгоритме выполняется условие покоящейся жидкости, то есть этот алгоритм относится к классу "well-balanced scheme".
В параграфе 1.5 описанный выше алгоритм тестируется на задаче Римана о распаде разрыва. Для уравнений МВ данные задачи носят название задач о разрушении плотины. В первой части параграфа построены аналитические решения для задачи Римана для уравнений МВ, далее приведено сравнение численного решения в рамках РУМВ с аналитическим, показана сходимость по сетке и влияние параметра регуляризации на устойчивость и точность численного решения. В параграфах 1.6 и 1.7 алгоритм тестируется на двух известных задачах о течении жидкости над неровным дном.
В последнем и наиболее объемном параграфе 1.8 изучается задача Римана, описывающая распад разрыва над подстилающей поверхностью в виде ступеньки и уступа дна. Этот тип задач значительно более сложен, чем задача о распаде разрыва над гладкой поверхностью. Решение любой задачи Римана можно представить в виде волн расширения, ударных волн и стационарных разрывов. При наличии ступеньки или уступа возникает дополнительный стационарный разрыв, который располагается над границей уступа или ступеньки. Есть целый ряд работ (см. [16]), которые посвящены аналитическому решению задачи Римана над ступенькой. Но полностью аналитически данная задача не решена, в отличии от известной задачи Римана для плоского дна. В данном параграфе 1.8 построены аналитические решения для пяти вариантов задачи. Показана сходимость численных решений РУМВ к аналитическим при сгущении пространственной сетки.
Во второй главе выполнено расширение численного алгоритма для моделирования течений, в которых возможно появление зон с нулевым уровнем жидкости - так называемых зон сухого дна. Условие сухого дно первоначально построено для покоящейся жидкости (параграф 2.1).
В параграфе 2.3 алгоритм определения границы сухого дна используется для расчета движущейся границы жидкости на примере распада одномерного разрыва (задача разрушения плотины), где изначально справа расположена зона с сухим дном. Оценки точности численного решения выполнены путем его сравнения с полученным автором точным решением задачи. В параграфах 2.4 и 2.5 рассмотрен класс одномерных задач с постоянным наклоном дна, которые моделируют профиль береговой зоны. Задача из параграфа 2. состоит в моделировании многократного набегания и сбегания волны с наклонного берега. Одно из аналитических решений, полученных в работе [43], автор сопоставяет с численными расчетами. Вторая задача (параграф 2.5) используется для моделирования характерных особенностей набегания одиночной волны цунами на берег с постоянным наклоном. В обоих примерах проведено сравнение точного и численного решений и показана монотонная сходимость численного решения к эталону. В главах 3 и 4 полученный алгоритм для расчета течений с сухим дном обобщается для прямоугольных и неструктурированных сеток.
В третье главе автор проводит обобщение построенного им алгоритма численного решения РУМВ на случай пространственных течений с использованием двумерных прямоугольных сеток в декартовой системе координат.
Численный алгоритм для расчета двумерных течений строится по аналогии с алгоритмом расчета одномерных течений. Система РУМВ для двумерного течения аппроксимируется с помощью метода конечного объема, причем все пространственные производные аппкросимируются центральными разностями со вторым порядком точности. Метод в целом аналогичен методу решения КГД-уравнений для двумерных течений (см., например, [38] и [44]).
В этом алгоритме учитывается выполнение условия покоящейся жидкости и приводится способ решения задачи для случая появления зон сухого дна. В параграфе 3.3 полученный РУМВ-алгоритм тестируется на известном примере о разрушении несимметричной дамбы. В параграфах 3.4 и 3.5 проводится численное моделирование двух экспериментов, выполненных в лабораторных условиях и моделирующих реальные физические явления. В параграфе 3. рассмотрена задача о набегании цунами на берег сложной формы. Численные расчеты выполнены в соответствии с данными натурного эксперимента.
Для постановки эксперимента строилась модель, в основе которой лежал реальный ландшафт береговой линии в соотношении 1 : 400. В этом эксперименте моделировалось цунами Окушири (яп. Okushiri tsunami), которое произошло в 1993 году в долине Монай (Monai Valley). Его характерной особенностью стало необычно большой размер береговых волн, размер которых на пике составил 31,7 метров. Соответствующий эксперимент был проведены в Научно-исследовательском институте электроэнергетики города Абико, Япония (Research Institute for Electric Power Industry in Abiko, Japan). Постановку задачи и результаты эксперимента можно посмотреть на ресурсах [45], [46].
В параграфе 3.5 проведено численное моделирование распространения волны прорыва в расширяющемся канале. Данный расчет выполнен в целях верификации алгоритмов для численного моделировании течений, возникающих при разрушении реальных шлюзовых камер и других гидротехнических сооружений, проводимых в вычислительном отделе Центра гидравлических исследований ОАО НИИЭС РусГидро. Для оценки точности численного метода использовались данные натурного эксперимента, выполненные в лабораторных условиях [47], а также численные расчеты задачи методом Годунова I и II порядка точности.
В четвертой главе диссертации построена разностная аппроксимация ре- гуляризованных уравнений МВ для неструктурированных сеток. Для построения численного алгоритма используется метод конечного объема, аналогичный описанному в Главе 3. Подробно изложен способ эффективной программной реализации численного алгоритма (параграфа 4.2), особенностью которой является подход, в котором для каждого узла сетки вычисляется и сохраняется суммарное значение всех потоков, втекающих в рассматриваемый контрольный объем. В параграфе 4.3 проведена модификации численного алгоритма, обеспечивающая выполнение условия покоящейся жидкости.
В качестве тестовых задач рассмотрено течение, возникающее при распаде цилиндрического столба жидкости (параграф 4.4) и течения, возникающего при разрушении плотины и затоплении поверхности с тремя конусами разных размеров (параграф 4.5).
На основе разработанных автором одномерных и двумерных программ другими авторами был выполнен расчет формирования уединенной волны в ветровом гидроканале [51] и решена задача о расчете нагрузок на стенки бака при колебаниях находящегося в нем слоя жидкости [52].
Основные результаты по теме диссертации изложены в 5 печатных изданиях [79], [80], [81], [82], [83] 4 из которых изданы в журналах, рекомендованных ВАК [79], [80], [81], [82]. Кроме этого научные результаты изложены в трудах международных конференций [84], [85] и одном препринте [83].
Результаты, представленные в диссертации, докладывались автором на следующих международных научных конференциях:
Численные методы в динамике жидкости (ICFD 2010), Университет Рединга, Великобритания, 12–15 апреля, 2010;
9-ая международная конференция по городскому сейсмостойкому строительству (9CUEE) и 4-ая азиатская конференция по сейсмостойким строениям (4ACEE), Токийский технологический институт, Токио, Япония, 6–8 марта, 2012;
XIX Международная конференция студентов, аспирантов и молодых учёных Ломоносов-2012, МГУ им. М.В.Ломонсова, Москва, 9–13 апреля, 2012;
6-ая европейская конференция по численным методам в прикладной науке и технике (ECCOMAS 2012), Венский университет, Австрия, 10– 14 сентября, 2012;
XX Международная конференция студентов, аспирантов и молодых учёных Ломоносов-2013, МГУ им. М.В.Ломонсова, Москва, 8–12 апреля, 2013;
Международная конференция Потоки и структуры в жидкостях, Российский государственный гидрометеорологический университет, СанктПетербург, 25–28 июня, 2013;
Суперкомпьютерные технологии математического моделирования (SCTEMM 2013), Якутск, 8–11 июля, 2013.
Результаты научной работы также докладывались автором и обсуждались на научном семинаре кафедры математики физического факультета МГУ им. М.В.Ломоносова под руководством профессора А.Н. Боголюбова, научном семинаре Санкт-Петербургского государственного морского технического университета, 30 января 2013 г., и на научном семинаре ИПМ им.
М.В.Келдыша под руководством проф. В.Ф. Тишкина и А.А. Кулешова.
Работа поддержана грантами РФФИ 10-01-00136, 13-01-00703а. Научная работа автора также стала победителем в конкурсе работ талантливых студентов, аспирантов и молодых ученых МГУ имени М.В.Ломоносова, учрежденный О.В. Дерипаска 2012 г.
Глава 1. Регуляризованные уравнения мелкой воды и метод численного решения 1.1. Уравнения мелкой воды Необходимо отметить, что уравнения мелкой воды можно получить несколькими способами. В книге [54] уравнение мелкой воды для ровного дна получается методом возмущений в первом приближении, когда рассматривается задача о течении идеальной, несжимаемой жидкости в канале малой глубины. В книге [6] уравнения мелкой воды получаются из простых предположений:
1) Пренебрегают компонентой вертикального ускорения 2) Пренебрегают вертикальными силами, кроме силы тяжести 3) Амплитуда колебаний жидкости мала по сравнению с глубиной жидкости К получившимся уравнениям нужно еще добавить уравнение, которое играет роль уравнения неразрывности для мелкой воды. Оно выводится из общих физических соображений. Например, его вывод этим способом можно найти в [6] и [1].
Рис. 1.1. Рисунок для иллюстрации используемых обозначений высоты уровня жидкости h(x, y) и профиля дна b(x, y) Также уравнения мелкой воды без учета вязкости и неровностей дна можно получить из двумерных уравнений Эйлера где - плотность жидкости, p - давление, u - скорость течения жидкости.
Проведем следующую постановку, используя баротропное приближение Если через h обозначить высоту уровня жидкости, то получим систему уравнений мелкой воды для случая плоского дна b = const. Интересно отметить тот факт, что в основе вывода КГД уравнений способом, который применяется в данной работе, лежат уравнения Эйлера (см. [38]). Поэтому можно получить регуляризованные уравнения мелкой воды с помощью баротропного приближения, использую указанную подстановку для квазигазодинамических уравнений. Однако это способ обладает очевидными недостатками: игнорируется случай неровного дна.
Касательно модели мелкой воды, следует отметить, что существуют модели с кинематической вязкостью [53]. Эти модели отличаются сложностью и разнообразием, и для широкого круга задач модели мелкой воды с членами, отвечающими за вязкость, не используются. Поэтому будем использовать уравнения мелкой воды в том виде, в каком они даны в книге [5].
где h(x, y, t) - высота уровня жидкости над уровнем дна, b(x, y) - функция, описывающая форму дна, ux (x, y, t), uy (x, y, t) - скорость течения жидкости (см., например, [5], [7]). На рис. 1.1 показано, что каждая из переменных h и b призвана обозначать. Например, в этом случае возмущение свободной поверхности находится как сумма = h + b. В практических задачах обязательно присутствуют внешние силы fi, которые учитывают, например, следующие факторы: силу Кориолиса, воздействие ветра, шероховатость и трение о дно.
В модели мелкой воды аналогом числа Маха в газовой динамике Ma = |u|/c, где c = RT скорость звука в газе, является число Фруда F r = |u|/c.
При этом скорость распространения малых возмущений вычисляется как c = 1.2. Регуляризованные уравнения мелкой воды Регуляризованные уравнения мелкой воды были получены по аналогии c КГД уравнениями. По поводу КГД системы следует обратиться к источникам [37], [38] и [74]. Интересно, что под КГД уравнениями можно понимать квазигазодинамические уравнения или квазигидродинамические уравнения.
В книге [38] на примере уравнений Навье-Стокса был показан формальный способ получения КГД уравнений. Однако он не является единственным способом их получения. В том же источнике можно найти вывод КГД уравнений на основе уравнения Больцмана. Сам метод аналогичен способу получения уравнений Эйлера, и его модификации - уравнения Навье-Стокса, из уравнений для функции распределения.
В качестве исходной системы уравнений возьмем систему уравнений мелкой воды для общего случая (1.2), (1.3), (1.4). Перепишем систему в других обозначениях и используем усреднение по времени Предполагается, что сила fi не успевает сильно измениться за промежуток времени (t, t +t). Это может быть справедливо для силы трения о дна водоема или воздействия ветра. В записи уравнений для краткости использовали тензор ij с компонентами Будем предполагать, что интервал усреднения (t, t+t) настолько мал, что в качестве средних значений можно взять величины с одного временного слоя, отвечающего моменту времени t (t < t < t + t). Эти величины также обозначим звездочками, то есть h (xi, t) = h(xi, t ), и ui (xi, t) = ui (xi, t ).
Далее используем разложение по времени Здесь для величин h и ui имеет одно и то же значение. Есть определенное удобство в том, чтобы использовать те же обозначения, которые использовались ранее для КГД систем. Поэтому обозначим поток h ux переменной jmx. Для квазигазодинамической система данная величина обозначала поток массы. Следует учитывать, что после проведения выкладок, отбрасываются все члены, которые содержат 2. В качестве примера рассмотрим следующие выкладки Используем систему уравнений (1.2), (1.3), (1.4), чтобы переписать производную по времени через пространственные производные. Тогда величина jmx принимает вид где В дальнейшем отбрасываем члены O( 2 ), и вместо h ux используем jmx. Проводим аналогичную процедуру, чтобы получить jmy где Таким же образом поступаем с величиной ij. При этом выражение для ij расписываем таким образом, чтобы была полная аналогия, если перейти от КГД уравнений к регуляризованным уравнениям мелкой воды (в случае плоского дна) с помощью баротропного приближения.
Выражения для ij имеют вид первых пространственных производных с коэффициентом Для КГД уравнений аналогичную конструкцию ij записывают в несколько другом виде. Это упрощает их запись, т.к. длинное выражение разбивается таким образом на части. ij записывается в виде Здесь переменными wx, wy и R обозначим следующие выражения Сразу отметим, что конструкция ij получилась несимметричная. Поэтому данная величина не несет такого же физического смысла, как тензор напряжения. Но противоречия в этом нет. Исходный тензор ij был симметричным, поэтому в данном случае симметричной остается величина ij = uj jm,i ij +ij 2 gh2. При этом компоненты по отдельности становятся несимметричными.
Проведя вышеизложенную процедуру, получаем новые величины, которые подставляем в усредненные уравнения. Как финальный штрих, отбрасываем слагаемые порядка O( 2 ). Выпишем получившуюся систему дифференциальный уравнений регуляризованной модели мелкой воды где Очевидно, что выписанные уравнения справедливы для декартовой системы координат. Регуляризованные уравнения мелкой воды путем замены переменных можно переписать для сферических и полярных координат. Тогда их можно использовать для моделирования течений на сфере.
Выписанная система тесно связана с системой уравнения мелкой воды (1.2), (1.3), (1.4). Почти очевиден тот факт, что стационарные решения уравнений мелкой воды является также стационарными решениями регуляризованных уравнений мелкой воды. Дело в том, что слагаемые с коэффициентом появляются, когда мы выражаем производные по времени с помощью уравнений (1.2), (1.3), (1.4) через пространственные производные. Соответственно, в стационарном случае эти выражения обращаются в нуль, а значит выражение при также обнуляются. Обратное утверждение не является очевидным, и для него также нужно специальное рассмотрение. Самым простым стационарным решение является задача о покоящемся озере : для стационарной задачи при ux = uy = 0 в отсутствии внешних сил fx = fy = 0 система уравнений МВ удовлетворяет условию гидростатического баланса то есть в этом случае уровень воды в озере горизонтален h(x, y) + b(x, y) = const. Для системы уравнений МВ с добавками указанное решение также удовлетворяется, что может быть проверено его непосредственной подстановкой.
Уравнения мелкой воды в случае плоского дна b = const и при отсутствии внешних сил fx = fy = 0 называются уравнениями Сен-Венана.
Для регуляризованных уравнений Сен-Венана доказан ряд теорем. В работах [71], [58] доказана теорема о неубывании специфической энтропии, что указывает на диссипативный характер регуляризованных уравнений мелкой воды. Также в работе [58] доказана теорема единственности. Согласно этой теореме единственности классическое решение начально-краевой задачи для квазигидродинамических уравнений в приближении мелкой воды (см. параграф 1.3) при отсутствии потоков через границу рассматриваемой области является единственным при любом выборе положительных параметров кинетической вязкости и малого параметра.
1.3. Регуляризованные уравнения для течений при малых числах Рассмотрим другой способ усреднения по времени исходных уравнений мелкой воды (1.5) и (1.6). Предположим теперь, что временной промежуток настолько мал, что толщина жидкости за это время не успевает сильно измениться, но успевает измениться только скорость. В результате для осредненных по времени величин 1.15 можно использовать следующие выражения Используем разложение по для постановки в усредненную систему уравнений 1.5, 1.6. Тогда для h ui получим выражение Как и в первом случае используем систему уравнений (1.2), (1.3), (1.4), чтобы выразить производные по времени через пространственные производные.
Тогда величины jmx, jmx принимают вид где Аналогичный подход используем для остальных членов уравнений. После преобразований регуляризованная система уравнений мелкой воды принимает вид:
В качестве компонент ij обозначены следующие выражения Выполняется условие для покоящейся жидкости над неровной поверхностью.
Как уже было сказано, при таком подходе мы считаем, что величина h при усреднении не изменяется. В этом случае градиент от величины b можно перенести в левую часть уравнения и сгруппировать его с произдводной от h.
При подстановке условий h + b = const и ux = uy = 0, fx = fy = 0 все члены уравнений будут равны нулю. Численные расчеты на основе полученных в этом параграфе уравнений можно найти в статье [56].
По результатам численного моделирования для расчетов реальных течений и задач предпочтительнее использовать квазигазодинамические уравнения и полученные на их основе регуляризованные уравнений мелкой воды.
Тем не менее квазигидродинамические уравнения представляют собой интерес с точки зрения аналитической теории, так как их особенности и свойства легче изучать. Например, доказана единственность классического решения основной начально-краевой задачи для квазигидродинамических уравнений [57]. Аналогичное утверждение доказано для квазигидродинамических уравнений в приближении мелкой воды [58]. В приложении Б будет подробно рассказано про квазигидродинамические уравнения и одномерные численные тесты.
1.4. Численный метод для одномерных течений В этом разделе подробно рассмотрим построение численного метода на основе регуляризованных уравнений мелкой воды в одномерном случае. Выпишем систему уравнений где Численный алгоритм также строится по аналогии с алгоритмом для квазигазодинамических уравнений. Получается явная по времени разностная схема, в которой используется интегро-интерполяционный метод с аппроксимацией пространственных производных в потоках центральными разностями.
Значения искомых переменных h(x, t) и u(x, t) задаются в узлах пространственной сетки i. Значения переменных в полуцелых пространственных точках i+1/2 вычисляются как среднее арифметическое значений в соседних точках hi+1/2 = 0.5(hi +hi+1 ), ui+1/2 = 0.5(ui +ui+1 ), bi+1/2 = 0.5(bi +bi+1 ), (1.21) Используя эти значения, вычисляют потоки в полуцелых узлах где Таким же образом определяют В качестве следующего этапа численного метода аппроксимируют первое уравнение (1.17) где индекс k соответствует текущему слою по времени, а шаг по времени составляет t. Все пространственные производные вычисляются на временном шаге k. Именно уравнения (1.17), (1.18) аппроксимируются интегроинтерполяционным методом. Это может показаться незаметным в одномерным случае. Но интегро-интерполяционный метод лежит в основе построения численного алгоритма на пространственных сетках.
Уравнение для определения расхода (1.18) аппроксимируется следующим образом:
Здесь Выпишем аппроксимированное уравнение для расхода (1.26), полагая ui = 0, fi = 0. Тогда получаем простое выражение Рассмотрим задачу о "покоящемся озере". Для нее выполняется hi + bi = const. Очевидно, что для точного выполнения нужно hi = 0.5(hi1/2 + hi+1/2 ).
Данное равенство не всегда выполняется, поэтому мы получим более аккуратные результаты, если вместе hi в выражении (1.27) будем использовать другое Устойчивость численного алгоритма обеспечивается слагаемыми с коэффициентом, величина которого связывается с шагом пространственной сетки x и вычисляется в виде где 0 < < 1 числовой коэффициент, выбираемый из условий точности и устойчивости счета. Условие устойчивости имеет вид условия Куранта, где шаг по времени выбирается по формуле Здесь число Куранта 0 < () < 1 зависит от величины параметра регуляризации и подбирается в процессе вычислений для обеспечения устойчивости численного решения [70].
1.5. Задача Римана 1.5.1. Построение автомодельного решения В параграфе 1.4 будем разбирать, как построить аналитическое решение для задачи Римана над уступом. Поэтому разберемcz более подробно, как решается задача Римана без уступа. Решение данной задачи можно найти в книге [6]. Постановка задачи выглядит довольно просто. Используем приближение мелкой воды для одномерного случая, дно плоское. Уравнения принимают вид Имеется область длиной L, где в начальный момент времени находится разрыв уровня жидкости, среда покоится. Слева от разрыва h1, справа h2. Для таких задач можно построить аналитическое решение. Обратимся к рис. 1.2, Рис. 1.2. Схематическое изображение распределения чтобы продемонстрировать поведение решения со временем. На нем схематично показано распределение h в какой-то момент времени. Есть три области I, III, IV, в которых величина h постоянна. В области II решение имеет вид u = U ( x ) и h = Z( x ). Точки A, B и C движутся с разными скоростями.
Сначала находим общий вид решения, который справедлив для области II, для этого делаем простую замену X = x. Система уравнений в частных проt изводных превращается в простую систему дифференциальных уравнений, общее решение которой имеет вид В общем решении фигурирует константа c. Известно, что в точке, жидкость еще не начала двигаться, значит в этой точке X = c или x = ct. Значит величина c является скоростью движения точки. Но нам также известно значение h в точке, которое равно h1. Получаем соотношение h1 = 9g (2c+c)2, из которого находим скорость c = gh1. В области III значения величин постоянны, их обозначают h и u. Эти величины связаны с h1 следующим образом. Найдя константу c = gh1, запишем вид решения в области II Избавляясь от переменной x, получаем связь Нужно также использовать условия на разрыве (аналог условий Гюгонио из газовой динамики для ударной волны в приближении мелкой воды, см. [6]).
N обозначает скорость разрыва, другими словами, скорость точки. Величины слева и справа от разрыва обозначаются и +. Соответственно, слева h = h и u = u, справа - h+ = h2 и u+ = 0. Формулы приводим к удобному для решения виду используем в дальнейшем вместо (1.35). Прежде чем находить неизвестные величины величины N, h, u и т.д., задачу обезразмериваем по следующей схеме. Вводим характерные скорости и безразмерные переменные c, u и N.
Получаем систему уравнений Используя систему, исключаем c2 и получаем квадратное уравнение для u. В результате величины c и u выражаются через N и b.
Уравнение для N получаем, используя связь между h, u и h1. Выражение, связывающее эти величины (1.34), в безразмерном виде 2 + u = 2. Выпишем получившееся нелинейное уравнение, к решению которого сводится данная задача.
Еще нужно упомянуть про один неочевидный момент: скорость точки B. Есть выражение u = 2 x + gh1. Из него находим скорость uB = 3 u gh1.
Наконец из уравнение (1.37) можно найти N даже для случая сухого дна, когда h2 = 0 и b = 0. Если b устремить к нулю, тогда N = 2.
1.5.2. Результаты численного расчета В большинстве случаев численный метод для газовой динамики или гидродинамики начинают тестировать с одномерной задачи распада разрыва. В предыдущем разделе описан способ получения аналитического решения для такой задачи, которая еще называется задачей о разрушении плотины.
Эта задача использовалась в [5] и [42] для тестирования численных методов решения уравнений мелкой воды на основе разностных алгоритмов, с использованием схем Годунова первого и второго порядков точности. Эта же задача исследовалась в других работах, среди которых есть [59], где тестируются схемы с хорошей балансировкой, построенные на основе кинетического подхода, и [60], посвященной построению конечно-разностных алгоритмов произвольного порядка точности. В [59] также отмечено, что большинство вычислительных алгоритмов, основанных на аппроксимации потоков направленными разностями, не позволяют рассчитывать указанный тест.
Рассматривается плоское одномерное течение жидкости в канале длины L с плоским дном b = const. В начальный момент в центре области задается разрыв уровня воды, разделяющий два однородных состояния с высотой уровня h = hl слева от разрыва и h = hr справа от разрыва. В начальный момент времени справа и слева от разрыва жидкость неподвижна, ul = ur = 0.
Были рассчитаны варианты распада разрыва для значений hr /hl от 0.5 до 0.0001, однако далее представлен лишь расчет, соответствующий статье [42], где L = 2000 м, hl = 10м, hr = 0.1м, g=9.8м/с2. Результаты приведены для момента времени t = 50 с.
На рис. 1.3 представлены профили высоты жидкости h(x), полученные на последовательности равномерных пространственных сеток с шагами x = 0.5, 1 и 2м. (сплошная линия - автомодельное решение). Фрагмент общей картины представлен на рис. 1.3, где узлы сетки отмечены маркерами – x = 2м. – треугольники, x = 1м – квадраты, x = 0.5 м – кружки. Вычисления Рис. 1.3. Задача о распаде разрыва. Распределения уровня жидкости и его фрагмент для x = 0.5, 1, 2 и = 0.1.
проведены для параметра регуляризации (1.30), вычисленного при = 0.1.
Графики демонстрируют монотонную сходимость численного решения к автомодельному при сгущении пространственной сетки.
Зависимость численного решения от величины параметра регуляризации показана на рис. 1.4 (фрагмент). Приведены профили высоты жидкости h(x), полученные на сетке с шагом x = 2м для последовательности параметров регуляризации (1.30), вычисленных при = 0.05, 0.1, 0.2 и 0.3. Шаг по времени (1.31) был равен t = 0.02сек, что соответствует числу Куранта = 0.1. С ростом разрыв решения сглаживается, при уменьшении в численном решении появляются осцилляции ( = 0.05), и при дальнейшем уменьшении величины регуляризатора решение теряет устойчивость. Опыт расчетов показывает, что имеется область оптимальных значений параметра, числа Куранта и шага сетки, обеспечивающих требуемую точность решения задачи.
Аналогичная задача изучалась в статье [56]. В ней использовался другой численный метод, основанный на другом виде уравнений мелкой воды с Рис. 1.4. Задача о распаде разрыва. Фрагмент распределения уровня жидкости для x = 2; = 0.05, 0.1, 0.2 и 0.3.
регуляризаторами. Их получают способом, который упоминается в приложении Б.
Для = 0.1 ширина скачка занимает 4 точек разностной сетки (рис.
1.3) и при уменьшении величины ширину ударного перехода можно уменьшить до 3 точек разностной сетки. В других работах (которые используют более сложные численные алгоритмы) удается аппроксимировать разрыв на меньшем числе точек сетки; так, в [59] разрыв аппроксимируется на 3 – 4 точках разностной сетки, а использование обобщенного метода Римана (см. [42]) позволяет описать разрыв на 1 – 2 точках сетки. Однако указанные методы существенно превосходят предложенный алгоритм по числу вычислений для одного шага по времени в расчете на одну точку пространственной сетки.
1.6. Задача о транскритическом течении над неровностями дна Этот тест - классическая задача о транскритическом течении над неровным дном, где наблюдается изменение числа Фруда от значений, меньших единицы, до значений, больших единицы, включая сверхкритические значения. В зависимости от начальных и граничных условий в течении может или формироваться, или не формироваться неподвижный разрыв – гидродинамический скачок. Для течения без разрыва система уравнений МВ имеет точное решение вида Система уравнений МВ с регуляризатором также удовлетворяет этому точному решению, поскольку на этом решении все слагаемые с коэффициентом обращаются в ноль.
Согласно [59], рассмотрим два наиболее сложных случая транскритического течения с образованием гидродинамического разрыва и без него. Рассмотрим канал длиной 25м, форма дна которого описывается функцией Случай 1: течение без разрыва. Слева налагается граничное условие hu = 1.53 м2 /с, h/x = 0. Справа ставим условия сноса потока h/x = 0, u/x = 0. В качестве начального условия выбрано h + b = 0.4м, u = 0.
Рассчитанные значения hu и jm с высокой точностью соответствуют аналитическому решению этой задачи Q0 = jm = hu = 1.53 м2 /с. В зоне над неровностью дна имеется отличие величины hu от точного значения, которое составляет 0.001 м2 /с. Расчеты выполнены для значений параметра регуляризации = 0.6, числа Куранта = 0.05 и шагах сетки x = 0.0625м и 0.125м. Последнее значение совпадает со значением из [59]. При этом величина ошибки в вычислении значения hu над неровностью дна оказывается меньше, чем в [59], использующей громоздкий вычислительный газокинетический алгоритм. Известно, что аккуратное определение расхода над препятствием представляет большие вычислительные сложности, чем определение скорости или уровня жидкости в этой зоне.
Распределение уровня жидкости h + b и числа Фруда показано на рис.
1.5 для шагов сетки x = 0.0625м (сплошная линия) и 0.125м(маркерыквадраты). Видно малое отличие численных решений при сгущении сетки, что свидетельствует о достигнутой точности, а также транскритический характер течения: при переходе через холмик, который образует поверхность дна, наблюдается рост скорости, соответствующий переходу числа Фруда через единицу. (рис. 1.5).
Рис. 1.5. Течение без разрыва. Распределение уровня жидкости и числа Фруда. x = 0.125 и 0.0625.
Случай 2: течение с неподвижным разрывом. Формирование гидродинамического разрыва обусловлено следующими граничными и начальными условиями: слева на границе полагаем hu = 0.18м2 /с, h/x = 0, справа на границе требуем выполнения условия h = 0.33м и условия сноса для скорости u/x = 0. В качестве начального условия выберем h + b = 0.33м, u = 0.
При численном интегрировании этой задачи в зоне за разрывом появляется вычислительная неустойчивость - сеточные осцилляции решения. Для сглаживания этих нефизичных колебаний по аналогии с КГД алгоритмом для моделирования сверхзвуковых течений газа, выражение для тензора вязких напряжений (1.20) следует дополнить слагаемым Это слагаемое представляет собой дополнительный регуляризатор типа вязкости Навье-Стокса. Его наличие действительно приводит к тому, что численное решение сглаживается.
Расчеты выполнены для значений параметра регуляризации = 0.6, числа Куранта = 0.1 и шагов сетки x = 0.125м(штриховая линия) и 0.0625м(сплошная линия), см. рис. 1.6 и 1.7.
На рис. 1.7 справа представлены рассчитанные значения jm (маркерыквадраты) и hu. Точное решение задачи соответствует величине расхода Q0 = hu = 0.18 м2 /сек. Во всей области течения для обоих сеток jm = Q0. Расход hu также хорошо совпадает с аналитическим значением везде, за исключением узкой зоны вблизи разрыва, где наблюдаются осцилляции, область которых занимает 8–10 шагов разностной сетки и уменьшается при ее сгущении.
Распределение уровня жидкости и числа Фруда для указанных сеток показаны на рис. 1.7. Видно очень близкое совпадение результатов расчета уровня жидкости для обоих сеток (рис. 1.7, справа) и монотонная сходимость численных результатов к аналитическому решению для числа Фруда (рис.
1.7, слева). Для сетки x =0.125 м. максимальное значение числа Фруда составляет F rmax = 2.35, для сетки x= 0.0625м. будет F rmax = 2.48, аналитическое значение равно F rmax =2.78. Заметим, что, в отличие от результатов из [59] и цитируемых в ней работ, в представленных расчетах не наблюдается превышения величины числа Фруда над его аналитическим значением.
Согласно [59], соответствие численного и аналитического решений задаh+b Рис. 1.6. Течение с разрывом. Распределение уровня h + b. x = 0.125 и 0.0625.
чи в зоне над препятствием представляет собой наиболее сложную проблему практически для всех численных алгоритмов.
1.7. Задача об отражении поверхностных волн от подводной возвышенности В иностранных публикациях данную задачу также называют тестом Le Veque. В этой тестовой задаче (работа [42], а также [59], [61]) изучается поведение со временем слабого возмущения в покоящемся канале, дно которого описывается функцией b(x) = 0.25[cos(10(x 0.5) + 1], если |x 0.5| < 0. Рис. 1.7. Течение с разрывом. Слева - Распределение числа Фруда F r. Справа - распределение потоков jm и hu. x = 0.125 и 0.0625.
В начальный момент времени жидкость покоилась (u = 0). Поверхность задается следующим распределением, в которое включено малое возмущение.
Согласно [42] система (1.17) – (1.20) решается численно (при этом полагают g = 1) для двух значений возмущения = 0.01 и 0.2. Область расчета составляет [0, 1], шаг по времени определяется числом Куранта = 0.2. Параметр регуляризации = 0.2. Численный расчет идет до момента времени t = 0.7.
Результаты для = 0.01 и 0.2 и пространственным шагом сетки x = 0.001 и 0.0001 показаны на рис. 1.8 – 1.9. На них видна сходимость по сетке. Результаты, полученные данным методом для x = 0.001, соответствуют результатам из [42] и [59, 61], в которых x = 0.01.
Неровное дно задается условием (1.38). В данном случае должны просто получить h(x) + b(x) = 1 и u(x) = 0. Если использовать аппроксимацию, указанную в (1.29), то из-за численных погрешностей получившийся результат и точное решение будут совпадать с точностью 1012.
Рис. 1.8. Тест Le Veque = 0.01, слева – уровень жидкости, справа – скорость жидкости.
Рис. 1.9. Тест Le Veque = 0.2, слева – уровень жидкости, справа – скорость жидкости.
1.8. Автомодельное решение и численное моделирование задач Римана при наличии уступов дна 1.8.1. Введение Рассмотрим решение одномерной задачи о распаде разрыва над ступенькой или уступом дна. Пусть в начальный момент времени жидкость покоится.
Таблица 1. Данные для различных вариантов тестовых задач Любое решение задачи такого типа представляет собой комбинацию волн разрежения, ударных волн и стационарного разрыва (или скачка), расположенного на ступеньке. Известно, что таких комбинаций может быть множество. В частности в [8] приводится 56 конфигураций решения для одномерного разрыва на ступеньке, для которых возможны, в частности, следующие конфигурации. Если жидкость стекает со ступеньки, то в левой области может располагаться волна разрежения, а в правой волна разрежения с ударной волной или две ударные волны. Причем в последнем случае ударные волны могут быть правыми (распространяться направо), или может формироваться левая и правая ударные волны. Конкретный вид этих конфигураций решения зависит от начальных данных. Информация о конкретном виде конфигурации решения существенно упрощает его аналитическое построение. Нужную информацию можно получить либо опираясь на уже изученные задачи, либо используя численное решение задачи по приведенному выше алгоритму.
Для получения численного решения задачи используется описанная выше разностная схема, основанная на РУМВ. Длина рассматриваемого канала составляет 100м. Разрыв расположен в точке с координатой x = 50м. Начальные данные слева и справа от разрыва для рассматриваемых далее тестовых случаев приведены в Табл. 1.1. Для разных вариантов величины, выбираются индивидуально. Величина выбирается из диапазона 0.1 < < 0.9.
В первом варианте = 0.01, в остальных случаях число Куранта составляло = 0.1. Практика расчетов показывает, что увеличение коэффициента приводит к сглаживанию численного решения, а его уменьшение до определенных значений сначала приводит к появлению осцилляций вблизи разрывов, а в дальнейшем и к разрушению численного решения. Ослабление этих негативных эффектов достигается путем уменьшения числа Куранта.
Расчеты всех вариантов проведены на четырех сгущающихся сетках x = 0.1м, 0.05м, 0.025м, 0.0125м с числом узлов, соответственно, 1000, 2000, 4000 и 8000. На всех приведенных далее рисунках представлены расчеты на сетке 2000 ячеек с шагом x = 0.05м. Результаты всех расчетов представлены на момент времени tout = 2с.
Далее для пяти характерных тестов приведено построение аналитических решений и их сопоставление с результатами численного моделирования.
Первый вариант построения аналитического решения изложен достаточно подробно, остальные варианты строятся аналогично, и их изложение приведено кратко.
1.8.2. Тест Рассматриваемый нами случай (hL = 7м, hR = 0.01м, bL = 3м, u = 0м/c) относится к изученному в [7], [5] типу решений. Схематичное изображение конфигурации решения представлено на рис. 1.10, где разрыв помещен в точку x = 50м. Здесь и далее на схемах = h + b - высота уровня жидкости, включая отметку дна. Через N обозначена скорость правой ударной волны.
Обозначим через h,u значения величин слева от стационарного разрыва на ступеньке, а через h1, u1 значения справа от разрыва на ступеньке.
Для этих величин на разрыве над ступенькой выполняются законы сохранеhL, u= Рис. 1.10. Схематическое изображение решения первого варианта задачи.
Распределение высоты жидкости для разрыва над ступенькой высотой bL = ния массы и полной энергии В области [x1, x2 ] располагается волна разрежения, которая описывается автомодельным решением (см. [6]) Из свойств решения (1.41) получаем, что для любой пары значений (h, u ) в этой области выполняется условие Волна разрежения в области [x3, x4 ] также описывается автомодельным решением Из (1.43) для нахождения неизвестной величины c2 получаем условие Слева от ударной волны величины равны h2, u2, а справа равны hR, u0, где u0 = 0 - начальная скорость течения. Здесь N - скорость правой ударной волны. Для ударной волны в точке x5 выполняется следующие условия на разрыве:
Однако условий (1.40), (1.42), (1.44) и (1.45) не достаточно для решения данной задачи, включающей в себя 7 неизвестных, а именно: h, u, h1, u1, h2, u2, N, для которых выписано 6 условий. Чтобы число неизвестных не превышало число уравнений, поставим еще одно дополнительное условие. Зная значения толщины слоя жидкости и скорости h, u слева от разрыва, из (1.41) определяем скорость движения точки x2 :
В данной задаче предполагаем, что точка x2 = 0 не движется. Это эквивалентно требованию, чтобы слева от точки x2 число Фруда имело вид F r = |u|/ gh = 1. Такое предположение соответствует аналитическому решению, полученному в [7], [5] на основе метода характеристик и численному решению на основе метода РУМВ. С учетом предположения о неподвижности разрыва в точке x2 = 0 из уравнений (1.41) получаем Решив кубическое уравнение, найдем h1 и u1. Запишем полную систему нелинейных уравнений для неизвестных h2, u2, N :
которую можно упростить до одного уравнения. Из второго и третьего уравнения (1.47) выразим h2, u2 через N, после чего произведем подстановку в первое уравнение и получим Нам удобнее решать численно систему (1.47) без упрощений. Результаты численного расчета исходной системы (1.2 – 1.4) методом РУМВ удобно использовать в качестве начального приближения. Окончательная запись точного решения для теста 1 имеет вид h(x), u(x) = Отметим, что точное решение центрировано относительно начала координат, которое следует сдвинуть для его сопоставления с численным решением системы уравнений (1.2 – 1.4).
На рис. 1.11 приведено сравнение полученных точных и числовых результатов для распределения уровня жидкости и ее скорости. Видно очень хорошее согласие числового и аналитического решений во всей области течения. Незначительные отличия видны в области сопряжения волны разрежения и ударной волны над углом ступеньки в точке x2 и в области сопряжения волны разрежения и зоной постоянного течения в точке x3.
Рис. 1.11. Распределение толщины жидкости и скорости на момент времени tout = 2с для первого варианта задачи. Сплошная линия – точное решение, квадратные символы – результаты численного расчета с параметрами = 0.5, = 0.01, x = 0. 1.8.3. Тест Аналогичным образом построено автомодельное решение для второго варианта задачи. Для него начальные значения имеют вид bL = 3м, hL = 7м, bR = 0м, hR = 1м. Конфигурация решения для данных начальных значений приведена на рис. 1.12.
Жидкость стекается со ступеньки. В левой области находится волна разрежения, которая, как и в первом варианте, примыкает к стационарному разрыву в точке x2. В правой области находятся две правые ударные волны N1, N2 – скорости этих ударных волн. В первую очередь находим h1,u1, решив Рис. 1.12. Схематическое изображение решения второго варианта задачи.
Распределение высоты жидкости для разрыва над ступенькой высотой bL = уравнение (1.46).
Для нахождения неизвестных h2, u2, N1, N2 составляем систему нелинейных уравнений, которые получим из условий на разрыве для двух правых ударных волн в точках x3 и x4.
Тогда точное решение примет вид h(x), u(x) = h = h1, u = u1, если x2 < x < x3 = tN1, На рис. 1.13 приведено сравнение точного и числового решения, которые отвечают второму варианту задачи распада разрыва на ступеньке.
Рис. 1.13. Распределение толщины жидкости и скорости на момент времени tout = 2с для второго варианта задачи. Сплошная линия – точное решение, квадратные символы – результаты численного расчета с параметрами = 0.9, = 0.1, x = 0. 1.8.4. Тест Третий случай является одномерным разрывом на ступеньке для начальных значений bL = 3м, hL = 7м, bR = 0м, hR = 1м. Ступенька расположена слева, конфигурация решения показана на рис. 1.14.
Рис. 1.14. Схематическое изображение решения третьего варианта задачи.
Распределение высоты жидкости для разрыва над ступенькой высотой bL = Как и в двух предыдущих случаях, жидкость стекает со ступеньки.
Волна разрежения или центрированная волна Римана примыкает к стационарному разрыву, что означает неподвижность точки x2. В правой области за ступенькой находятся две ударные волны. Одна из них распространяется направо, а другая налево, но из-за наличия ступеньки слева она останавливается в точке x2, то ее скорость N1 = 0.
Из рассмотренных ранее примеров мы уже знаем значения h,u слева от разрыва. Также было выписано кубическое уравнение (1.46) для h1. Все выкладки, которые были выписаны для первого и второго варианта, остаются справедливыми. Нам также нужно найти h2,u2,N1,N2, для чего нужно решить систему (1.48). Главное отличие третьего варианта задачи от второго заключается в том, что скорость распространения ударной волны N1 меняет знак и при построении аналитического решения системы может обращаться в ноль.
Необходимо сделать важное замечание к этому варианту задачи. В [17], Таблица 1. Зависимость значений точного решения от величины диссипации K [18] указывается, что на ступеньке может происходить потеря энергии потока, например, за счет турбулизации течения перед ступенькой. В этом случае полная энергия (1.40) на ступеньке не сохраняется. Пусть K - некий коэффициент, который должен описывать потерю энергии на ступеньке за счет диссипации. Согласно [17], [18] добавим его в условие сохранения полной энергии на ступеньке (1.40):
Коэффициент K имеет размерность в метрах, а его добавление в уравнение (1.40) эквивалентно изменению высоты ступеньки. Значения коэффициента K заранее не известно. Покажем влияние коэффициента на аналитическое решение задачи. Для этого решим систему уравнений (1.48), (1.49), учитывая потерю энергии на ступеньке. Результаты решения представлены в Табл.
1.2, из которых следует, что подбором значения можно уменьшить скорость левой ударной волны вплоть до ее полной остановки. Таким образом, диссипация на ступеньке приводит к тому, что скорость левой ударной волны, набегающей на ступеньку, уменьшается и может обратиться в ноль. Получаем дополнительное условие, что величина K должна отвечать такой диссипации на ступеньке, при которой скорость левой ударной волны будет N1 = 0 и точка оказывается неподвижной.
Сведем вместе условия (1.48) и (1.49), учитывая N1 = 0. Исключая скорости, получаем систему из четырех уравнений:
где для удобства записи введен коэффициент При численном решении системы (1.50) требуется взять в качестве начального приближения величины, обеспечивающие выполнение неравенства h1 < h2.
Получим следующие значения: K = 1.8099м, h1 = 1.3622м, u1 = 12.6112м/с, h2 = 6.003м, u2 = 2.8617м/с, N2 = 8.5766м/с. Из полученного решения следует, что значение K = 1.8099м отвечает диссипации на ступеньке, при которой N1 = 0.
Заметим, что для этого варианта точка x = x2 является тройной точкой, в которой имеется три значения для высоты и скорости потока. Аналитическое решение задачи записывается в виде Сравнение точного решения с учетом диссипации, соответствующей условию N1 = 0, и числовых результатов приведены на рис. 1.15.
Рис. 1.15. Распределение толщины жидкости и скорости на момент времени tout = 2с для третьего варианта задачи. Сплошная линия – точное решение, квадратные символы – результаты численного расчета с параметрами = 0.9, = 0.1, x = 0. Из приведенных рис. 1.15 видно, что численное решение однозначно дает вариант, соответствующий диссипации энергии на ступеньке с величиной этой диссипации, соответствующей значению N1 = 0. Последнее определяется физичным характером введенного в систему уравнений мелкой воды регуляризатора и наличием для РУМВ теоремы о диссипации энергии (см. [70], [71], [72]). На рис. 1.15 ясно видна тройная точка, расположенная над краем ступеньки в точке x2. Естественно, что в числовом решении значения высоты и скорости жидкости, соответствующие этой тройной точке, смещены друг относительно друга на один шаг пространственной сетки. При этом их значения очень хорошо соответствуют величинам, полученным при аналитическом решении задачи (см. фрагмент на рис. 1.16).
Рис. 1.16. Распределение толщины жидкости в момент времени tout = 2с на отрезке [49.8, 50.2]. Сплошная линия – точное решение, толстая линия с квадратными символами – результаты численного расчета с параметрами = 0.9, = 0.1, x = 0. 1.8.5. Тест В четвертом варианте задачи ступенька располагается справа, разрыв на ступеньке задается значениями hL = 10м, bL = 0м, hR = 0.2м, bR = 3м.
Теперь жидкость натекает на ступеньку. В левой области образуется волна разрежения или волна понижения. В правой области ударная волна и волна разрежения, которая примыкает к стационарному скачку на ступеньке в точке x3. На рис. 1.17 приведена конфигурация решения.
Здесь на участке [x3, x4 ] сохраняется величина Из соотношений (1.41) и условия x3 = 0 получаем связь между h,u и h2,u2 :
Собирая вместе условия для ударной волны в точке x5, стационарного скачка в точке x3, а также условие сохранения инварианта для волн разрежения, Рис. 1.17. Схематическое изображение решения четвертого варианта задачи.
Распределение высоты жидкости для разрыва над ступенькой высотой bL = 3м для четвертого варианта задачи.
получаем систему Численно решив приведенную выше систему, находим аналитическое решение h(x), u(x) = h = h1, u = u1, если x2 < x < x3 = 0, На рис. 1.18 приведено сравнение аналитического и числового решения четвертого варианта задачи о распаде разрыва. Здесь также видно очень хорошее согласие полученных решений.
Рис. 1.18. Распределение толщины жидкости и скорости на момент времени tout = 2м для четвертого варианта задачи. Сплошная линия – точное решение, квадратные символы – результаты численного расчета с параметрами = 0.3, = 0.1, µ = 0, x = 0. 1.8.6. Тест В этом варианте жидкость натекает на ступеньку. Заданы начальные значения слева и справа от разрыва hL = 10м, bL = 0м, hR = 2м, bR = 3м. В левой области располагается волна разрежения, а в правой – ударная волна. На ступеньке расположен стационарный скачок. Схема конфигурации течения показана на рис. 1.19.
Рис. 1.19. Схематическое изображение решения пятого варианта задачи. Распределение высоты жидкости для разрыва над ступенькой высотой bL = 3м При составлении системы уравнений воспользуемся сохранением инварианта (1.42) на участке [x1, x2 ]. К этому добавим условия на разрыве x4 и условия для стационарного скачка x3.
Для решения системы (1.51), как и в предыдущих случаях, удобно использовать результаты численного расчета для поиска начального приближения.
В результате проведения выкладок, аналогично предыдущим задачам, для описания аналитического решения получаем функцию Сравнение точного решения с аналитическим приведены на рис. 1.20, где видно их близкое соответствие.
1.8.7. Оценка точности численного метода На рис. 1.21 приведена зависимость норм L1 (h) и L1 (hu) от шага пространственной сетки. Значения норм определены по следующим формулам:
где hex,uex - точные решения соответствующих задач, выписанные выше. Из графиков видно, что для всех рассмотренных тестов численные результаты монотонно сходятся к точным решениям при сгущении шага пространственной сетки. Заметим, что наиболее быстрая сходимость к точному решению наблюдается для третьего варианта, который является наиболее сложным и обладает тройной точкой в области неподвижного разрыва над ступенькой.
Скорость сходимости для всех вариантов пропорциональна шагу пространственной сетки, что соответствует теоретическим оценкам для схем первого Рис. 1.20. Распределение толщины жидкости и скорости на момент времени tout = 2с для пятого варианта задачи. Сплошная линия – точное решение, квадратные символы – результаты численного расчета с параметрами = 0.7, = 0.1, µ = 0, x = 0. порядка точности на гладких решениях. Однако следует заметить, что здесь рассматриваются разрывные решения, для которых, как правило, точность разностных схем, включая схемы высокого порядка, не достигает первого порядка.
0. 0. 0. 0. 0. 0. Рис. 1.21. Зависимость норм L1 (h) и L1 (hu) от шага сетки x для пяти тестов.
Цифрами на графиках обозначены номера соответствующих тестов.
Глава 2. Условие сухого дна для одномерных задач 2.1. Постановка условия для сухого дна на примере водоема с холмом и сухим верхом В очень многих практических задачах (разливы рек из-за наводнений, обрушение волн цунами на прибрежную территорию) присутствует сухое дно.
Это области, где толщина жидкости считается нулевой h = 0. Для взаимодействия с сухими областями разные численные методы используются разные подходы (см. [41] или [36], [42]). Для нашей работы мы стали использовать метод, изложенный в статье [41].
Исходим из того, что в области сухого дна жидкость покоится. Но из уравнений (!ссылка на рег. урав.) следует, что наличие неровного дна b() = 0, градиент которого входит в правую часть в качестве внешней силы, приводит жидкость в движение, которого быть не должно. Поэтому поставим условие, чтобы для малых h скорость равнялась нулю u = 0.
В нашем численном методе фигурирует параметр, который входит во все регуляризующие добавки. Из естественных соображений понятно, что в области сухого дна нам регуляризующие добавки не нужны, так как жидкость покоится. Ставим такое же ограничение для при малых h.
Параметр может быть сколь угодно большим. Поэтому нам нужно задать его нижнюю границу. Рассмотрим случай покоящейся жидкости на рис. 2.1.
Граница с сухим дном находится в точке i + 1. Наложим такие условие на, чтобы в точке i было u = 0. Для этого находим высоту жидкости в этой точке и ставим ее нижней границей для В случае прямоугольной или неструктурированной сетки получаем условие В случае прямоугольной сетки шаблоном является прямоугольник, а для неструктурированной сетки он представлен в виде выпуклого многоугольника. Мы укажем конкретный вид для каждой задаче в соответствующей главе.
Рис. 2.1. Выбор оптимального для задач с сухим дном Перейдем к самой простой задаче о водоеме с холмом и сухим верхом.
Эта задача рассматривалась в работах [65] и [67]. С помощью этой задачи проверяется выполнения условия покоящейся жидкости (в англоязычной литературе численный алгоритм в этом случае называют well-balanced scheme).
Под этим условием подразумевается, что изначально покоящаяся жидкость остается в покое, даже при наличии неровного дна. Это условие является достаточно простым и естественным, но нуждается в проверке для численного метода.
В области размером L = 1м форма дна задается выражением Рельеф дна показан на рис. 2.2. На нем видно, что холмик разделяет два резервуара с водой. В начальный момент t = 0 жидкость покоилась u = 0.
Рис. 2.2. Холм с сухим верхом Стационарный уровень жидкости находится на высоте 0.1м. В области сухого дна толщина слоя жидкости равняется h = 0 (h(x) = 0 если b(x) > 0.1m).
Начальные данные были поставлены так, чтобы выполнялось условие h(x) + b(x) = const. Граничные условия записывают в следующем виде Расчеты проводились для параметров = 0.5, = 0.5. Расчеты велись до момента времени t = 200с. На рис. 2.2 изображены результаты для пространственного шага сетки x = 0.001м. Было установлено, что жидкость остается в покое при выполнении двух условий. Во-первых, для всех точек должно выполняться При выполнении этих условия получаем, что для задачи о водоеме с холмом и сухим верхом полученная точность составляет 106.
2.2. Одномерный разрыв В дальнейшем мы будем рассматривать различные задачи, в которых в силу их природы появляется "сухое дно". Самой простой из подобного рода задач является одномерный разрыв, где с одной из сторон нет жидкости (т.е.
h = 0). Рассмотрим область [0, 50]м. В центре отрезка (точка с координатой x = 25м) расположен разрыв. В начальный момент времени t = 0с жидкость покоится (ul = 0м/с, ur = 0м/с).
В левой области высота уровня жидкости hl = 1.0м, справа сухое дно hr = 0м. Расчет ведется до момента времени t = 5с. Зная общее решение для задачи разрыва, легко записать аналитическое решение для данного случая.
Рис. 2.3. Схема решения. В правой области сухое дно, точки x1, x2 движутся в разных направлениях. Жидкость на отрезке [x1, x2 ] описывается автомодельным решением уравнений мелкой воды.
Хотя для этой задачи легко найти скорость передвижения береговой точки на границе с сухим дном x2 = 2t ghL, но на практике при численном расчете мы не можем определить точное расположение береговой точки, так как условие для сухого дна убирают из рассмотрения область в которой h <.
Расчеты проводились на трех различных сетках x = 0.05м (N = 1000), x = 0.025м (N = 2000), x = 0.0125м (N = 4000), x = 0.00625м (N = 8000). Остальные параметры расчета имели значения: параметр регуляризации = 0.2, число Куранта = 0.1. В задаче с сухим дном необходимо рациональным образом выбирать, для данной задачи подходит = 0.0001м.
Можно взять = 0.001м и большую величину, но тогда результаты расчета вблизи береговой точки будут еще хуже. Если выбрать = 0.00001м тогда численный алгоритм будет неустойчив и придется уменьшать до 0.01. Поэтому из всех вариантов мы выбираем самый оптимальный.
В численных расчетах мы можем модифицировать параметр, который входит в добавочные члены. Для модельной задачи о распаде разрыва с сухим дном мы взяли два различных выражениях для. Это позволяет нам сделать вывод о виде, который больше подходит для задач с сухим дном. Выбор можно сделать следующим образом. В первом варианте в качестве характерной скорости берется скорость распространения длинных волн c = gh, тогда Во втором случае в качестве характерной скорости берется c = gh + |u|, тогда На рис. 2.8 показана сходимость по сетке. Численное решение на самом деле совпадает с аналитическим на всем протяжении графиков, но заметное отличие видно только на границе с сухим дном. Из графиков также становится понятно, что уменьшение шага сетки сдвигает эту границу. Кроме шага сетки на положение береговой линии может влиять выбранный параметр (при наличии неровного дна мы привязываем его к размеру сетки) и величины вблизи этой граничной точки.
На рис. 2.5 приведены численные расчеты на одной сетке x = 0.05м (N = 1000) с различными (см. выражения (2.5), (2.6)). В первом случае график получается более сглаженным из-за большой величины для тех же параметров. Но в любом случае выбор не влияет существенным образом на положение границы с сухой области. Хотя в дальнейшем предпочтительней выбирать (2.6).
В задаче об одномерном разрыве с сухим дном жидкость натекает на сухую область. Более сложной является задача, где жидкость стекает с какойлибо поверхности, образуя сухую область. Эти гораздо более предпочтительные задачи (неровная поверхность дна, изменяющийся наклон берегового склона, стекание жидкости) мы рассматриваем, когда говорим об эффективном применении численного алгоритма для задач сухого дна.
Рис. 2.4. Сходимость по сетке. Слева приведены профили жидкости на отрезке [40, 44]. Справа высота жидкости на отрезке [40, 44.5]. Для всех вариантов выбранный момент времени t = 3с, остальные параметры = 0.0001, = 0.2, = 0.1.
2.3. Задача Римана с разбегающейся жидкостью Рассмотрим задачу Римана с разбегающейся в разные стороны жидкостью. Пусть задана область [0, 100]м. В центре отрезка x = 25м расположен разрыв. В начальный момент времени t = 0с высота уровня жидкости везде постоянна (hL = hR = 10м), но скорости различны (uR = uL = 25м/с). На рис. 2.6 схематически изображено решение рассматриваемой задачи.
В области [x2, x3 ] будет появляться сухое дно только в случае выполнения условия uR > 2 ghR. Отметим, что в точках x2 и x3 также расположен разрыв для скорости u. Значения скорости слева и справа от точки x различаются. Это вызвано тем, что именно в этих точка h переходит к нуt = 3s Рис. 2.5. Сравнение профиля скорости, полученного с помощью численного расчета, с аналитическим результатом для момента времени t = 3с. Параметры = 0.0001, = 0.2, = 0.1. Пунктирная линия - аналитическое решение, сплошная линия - численный расчет с формулой (2.5), пунктир с точками численный расчет с формулой (2.6) левому значению. Аналитическое решение состоит из двух волн разрежения и области с сухим дном. На отрезке [x1, x2 ] волна разрежения описывается автомодельным решением (см. 1.41) Аналогичным образом в области [x3, x4 ] решение описывается функциями Из начальных условий находим константы c1 и c Рис. 2.6. Схема решения. В области [x2, x3 ] образуется сухое дно. Точки x2, x3 движутся в разных направлениях.
Запишем аналитическое решение для данной задачи.
Время расчет t = 2c, параметр регуляризации = 0.5, число Куранта = 0.01. Шаг сетки составил x = 0.025м (N = 2000). Расчеты проведены для двух параметров отсечения = 0.1м и = 0.01м, чтобы показать влияние выбора на результаты численного расчета.
Из-за наличия численных погрешность величина не может быть нулевой. При нулевом значении любое численное возмущение для задач, в которых есть сухие области h = 0, приведет к отрицательным значениям h. С физической точки зрения при рассмотрении течения глубиной 10м, мы пренебрегаем движением жидкости в слое толщиной соответственно 10см ( = 0.1м) и 1см ( = 0.01м). Таким образом, параметр нужно также выбирать согласно целесообразности каждой отдельной задачи. Например, если требуется смоделировать набегание 10-метровой волны на берег, то мы можем пренебречь движением жидкости в 1см слое вблизи сухого дна. Десятиметровые волны могут разрушить здания и инфраструктуру, при таких масштабах течения можно пренебречь движением в лужах глубиной 1см.
Рис. 2.7. Распределение высоты уровня жидкости и скорости на момент времени tout = 1с. Сплошная линия – точное решение, квадратные символы – результаты численного расчета с параметром отсечения = 10см, окружности – результаты численного расчета с параметром отсечения = 1см Слева на рис. 2.7 показано распределение высоты уровня жидкости для момента времени t = 1c. Аналитическое решение обозначено сплошной линией. Квадратиками обозначено численное решение, полученное для параметра отсечения = 10см. Видно, что почти на всем протяжении численное решение и аналитическое решение хорошо согласуются. Справа на рис. 2.7 показано распределение скорости в момент времени t = 1c. Треугольниками обозначено численное решение, полученное при тех же параметрах расчета, но для параметра отсечения = 1см. Из рисунка видно, как величина отсечения влияет на положение границы жидкости с сухой областью.
Рис. 2.8. Сходимость по сетке. Слева приведены профили жидкости на отрезке [40, 44]. Справа высота жидкости на отрезке [40, 44.5]. Для всех вариантов выбранный момент времени t = 3с, остальные параметры = 0.0001, = 0.2, = 0.1.
2.4. Сравнение точного решения с численными расчетами для случая постоянного наклона дна В качестве тестовой задачи для численных расчетов данная задача рассматривалась в [69]. Для этого проводится сравнение с известным точным решением периодической волны Кэри и Гринспана (Carrier and Greenspan).
Это аналитическое решения часто применяют, чтобы проверить способность X_shore, m Рис. 2.9. Зависимость координаты "береговой точки"от времени для трех шагов сетки x = 0.1м, x = 0.05м, x = 0.025м численного алгоритма моделировать набегание и сбегание волны с наклонного берега. В частности, оно позволяет решить, как поставить граничное условие для береговой линии. Вид решения взят из работы [69], в которой приводится ссылка на статью с решением [43]. Аналитическое решение удобно записать с помощью безразмерных переменных.
Решение записывается в неявном виде Из вида решения следует, что в момент времени t = 0 жидкость покоилась.
Движение периодическое, период которого T = L/(g · tg). Величина является безразмерным аналогом скорости длинных волн c = gh. Поэтому подстановка = 0 позволяет получить информацию о координате x и скорости u в береговой точке. Для справедливость выражения необходимо, чтобы величина A была ограничена 0 A 1. Рассмотренная задача соответствует случаю A = 0.6, g = 9.8м/с2, длине области L = 20м, углу наклона tg = 1/30. Областью расчета является отрезок [100м, 10м]. Дно имеет постоянный наклон b(x) = tg · x. В качестве начальных условий взято решение при t = 0. Слева на границе области поставлено периодическое по времени граничное условие для h и u, которое получим из точного решения при фиксированном x0 = 100м.
Расчет проводился для трех шагов сетки x = 0.1м, x = 0.05м, x = 0.025м. Параметры численного алгоритма имели значения = 0.2, = 0.1 для всех шагов сетки. Так как дно с постоянным наклоном, то для трех разных шагов брали три разные значения = x · tg, которые соответственно равны = 1/300м, = 1/600м, = 1/1200м.
Далее приведена информация о полученных результатах численного расчета. На рис. 2.9 показана зависимость координаты "береговой точки"(moving shoreline), которая обозначает разделение жидкости и сухого дна, в зависимости от времени. Взято изменение координаты береговой точки для трех периодов 3T, значение периода нам известно T = L/(g · tg).
Точное решение обозначено сплошной линией, результаты численного расчета для трех шагов сетки помечены различными прерывистыми линиями. На графике шаг сетки x мы обозначаем с помощью hx. Соответствующим образом обозначены линии для пространственных шагов x = 0.1м, x = 0.05м и x = 0.025м.
Изменение шага сетки, которое влечет изменение = x·tg, не сильно влияет на положение береговой точки. На левом рис. 2.9 изменение, связанные с уменьшением шага сетки, заметны только около пиковых значений. На правом рис. 2.9 показано, что уменьшение шага сетки приближает результаты численных расчетов к аналитическому решению. Хотя сам численной метод не позволяет исследовать течения в непосредственной близости от береговой точки, потому что из рассмотрения выкидывается область {x : 0 < h(x) < }.
На рис. 2.10 показано распределение скорости u для момента времени t = 5с. Здесь также приведено аналитическое решение и результаты численного расчета для трех шагов сетки x = 0.1м, 0.05м, 0.025м. Остальные параметры численного алгоритма остаются неизменными = 0.2, = 0.1.
На левом рис. 2.10 приведено распределение по скорости для всего отрезка [100м, 10м]. В таком масштабе численные результаты совпадают с аналитическим решением. Разница заметна только в граничной точке. На правом рис. 2.10 показан увеличенный фрагмент распределения скорости в окрестности "береговой точки". Из вида аналитического решения следует, что скорость терпит разрыв около этой береговой точки. Используемый численный метод не позволяет определить точное положение разрыва, которое совпадает с положением "береговой точки". Также в данном численном методе нет возможности специальным образом отслеживать перемещение разрыва.
В результате, скачок по скорости оказывается сглажен. При этом на правом рис. 2.10 продемонстрировано, что численное решение приближается к аналитическому при уменьшении шага сетки x и соответствующем уменьшении.
Из представленных результатов уже можно сделать некоторые выводы о рассматриваемом численном методе. Во-первых, условия для сухого дна, которое имеет вид u = 0 при h < является самым адекватным граничным условием. Потому что значение u = 0 просто отвечает значению справа от разрыва скорости около береговой точки. Во-вторых, наблюдается расхождение в скорости между аналитическим решением и численным расчетом в окрестности той же береговой точки. Это обусловлено тем, что в нашем численном расчете всякий разрыв получается сглаженным, даже если в задаче нет физических составляющих, отвечающих за диссипацию. Кроме того, мы не можем точно определить положение береговой точки, так как не рассматривается движение в области {x : 0 < hx < }.
только для целых, но и полуцелых точек. К тому же величину оцениваем из предположения о стационарности жидкости. Тогда глубина жидкости в окрестности береговой точки x0 (L = b(x0 )) приближенно равна h(x) = L b(x) = L b(x0 ) x x. Получим оценку для, которая в расчетах может быть больше указанного значения, но того же порядка.
Рис. 2.10. Распределение скорости u в момент времени t = 5с для трех пространственных шагов x = 0.1м, x = 0.05м, x = 0.025м 2.5. Набегание цунами на наклонный берег Рассматриваемая задача взята из The third international workshop on long-wave runup models (Benchmark Problem №1). Постановка задачи и аналитические результаты находятся на электронном ресурсе [45].
На этом ресурсе можно найти данные о начальном распределении жидкости. Там же берем аналитические результаты по распределению возмущения свободной поверхности и скорости жидкости u для трех моментов времени t1 = 160с, t2 = 175с и t3 = 220с, а также информацию о движении береговой линии на промежутке времени [0, 355с]. Распределение жидкости в нулевой момент времени t = 0 нам известно из начальных условий. Они не записываются в аналитической форме, поэтому обратимся к рисунку 2.11.
Справа на рис. 2.11 изображена свободная поверхность жидкости и наклонный берег с постоянным углом наклона (tg = 0.1). На рис. 2.11 слева показана зависимость возмущения свободной поверхности жидкости от координаты. Область расчета занимает отрезок [500м, 50000м]. В начальный момент времени t = 0 жидкость покоится, и задано только возмущение свободной поверхности. На отрезке [0, 50000м] задаем распределение жидкость.
Область [500м, 0] отведена для сухого дна, в которой задано h = 0. Слева ставим условия для сухого дна. Справа держим постоянный уровень жидкости h(x = 50000м) = 5000м.
Рис. 2.11. Начальное возмущение свободной поверхности жидкости. На графиках разные масштабы по горизонтали и вертикали. Жидкость занимает область размером 50км, а разница между максимумом и минимум свободный поверхности, как видно из левого рисунка, составляет 10м.
Расчеты проводились для трех шагов сетки x = 5м, x = 2м и x = 1м. Для демонстрации влияния параметров численного алгоритма на результаты расчетов взято два значения параметров = 0.3 и = 0.4. При этом численные результаты для и u приведены для параметра = 0.4.
Другой параметр = 0.5 оставался неизменным. На рис. 2.12 приведена зависимость координаты береговой точки для различных параметров расчета x и по сравнению с точным решением. Графики получаются ступенчатые, поэтому нужно использовать процедуру сглаживания для получения информации о скорости движения береговой точки. Установлено, что значение = 0.3 является оптимальным для расчетов в данной задаче.
На группе рисунков (рис. 2.13, 2.14, 2.15) представлено сопоставление результатов численных расчетов и точного решения для трех моментов времени t = 160с, 175с и 220с. На каждом графике приведены результаты расчетов для уменьшающейся серии пространственных шагов x = 5м, x = 2м и x = 1м. Но разницы в таком масштабе между известным точным решением и численными расчетами незаметна.
Обратимся к другой серии графиков (рис. 2.16, 2.17, 2.18), которые представляют сравнение результатов расчета и известного точного решения распределения скорости для тех же трех моментов времени t = 160с, 175с и 220с. Численные расчеты для скорости также приведены в случае трех пространственных шагов x = 5м, x = 2м и x = 1м. Трем разным шагам сетки отвечают разные значения согласно условию = x · tg. Только здесь - угол наклона. Численные результаты для и u получены для параметров численного алгоритма = 0.4, = 0.5. Для графиков скорости u видно значительное расхождение в правой части отрезка полученных численным методом значений от точных. Это вызвано тем, что на рис. 2.16, 2.17, 2.18 изображен отрезок, примыкающей к береговой точке. Поэтому такие результаты являются допустимыми. Обратимся к предыдущему параграфу 2.4, в котором было показано, что величина скорости имеет разрыв/скачок в береговой точке. А результаты численного расчета для скорости u будут сглаживать этот скачок.
Рис. 2.12. Зависимость от времени координаты точки, которая отделяет жидкость от сухой области. Точное решение обозначено сплошной линией. Результаты расчетов для сеток с шагом x = 5м до x = 1м обозначены пунктирными линиями.
Рис. 2.13. Возмущение свободной поверхности жидкости в момент времени t = 160с Рис. 2.14. Возмущение свободной поверхности жидкости в момент времени t = 175с Рис. 2.15. Возмущение свободной поверхности жидкости в момент времени t = 220с Рис. 2.16. Профиль скорости в момент времени t = 160с Рис. 2.17. Профиль скорости в момент времени t = 175с Рис. 2.18. Профиль скорости в момент времени t = 220с Глава 3. Обобщение алгоритма на пространственные течения 3.1. Численный алгоритм для двумерной прямоугольной сетки Уже рассматривался численный алгоритм для одномерного случая. Теперь приведем построение для двумерного случая. Используем прямоугольную сетку, соответствующий шаблон показан на рис. 3.1. Величины в центре ячеек находим, использую значения в целых узлах. Выберем переменную h для демонстрации расчетных формул.
Рис. 3.1. Шаблон для прямоугольной сетки hi+1/2,j+1/2 = 0.25(hi,j + hi+1,j + hi,j+1 + hi+1,j+1 ) hi+1/2,j1/2 = 0.25(hi,j + hi+1,j + hi,j1 + hi+1,j1 ) hi1/2,j+1/2 = 0.25(hi,j + hi1,j + hi,j+1 + hi1,j+1 ) hi1/2,j1/2 = 0.25(hi,j + hi1,j + hi,j1 + hi1,j1 ) Зная значения в центре ячеек, находим значения на ребрах hi+1/2,j = 0.5(hi+1/2,j+1/2 + hi+1/2,j1/2 ) hi1/2,j = 0.5(hi1/2,j+1/2 + hi1/2,j1/2 ) hi,j+1/2 = 0.5(hi+1/2,j+1/2 + hi1/2,j+1/2 ) hi,j1/2 = 0.5(hi+1/2,j1/2 + hi1/2,j1/2 ) Аналогичная аппроксимация используется для всех остальных параметров кроме. Значение на ребре вычисляется как среднее значение от точек, которые это ребро соединяет.
Хотя величина и меняется от точки к точке, но остается малой величиной.
Поэтому способ ее аппроксимации играет мало существенную роль.
В первую очередь нужно найти значения потоковых величин на ребрах.
Так как выражение для j записывается с помощью поправки к скорости, то ее выражаем простым образом ji+1/2,j = hi+1/2,j (ux ji,j+1/2 = hi,j+1/2 (uy Здесь для удобства использован верхний индекс для обозначения x и y компоy y нент. Таким образом, нам нужны значения wi+1/2,j, wi1/2,j и wi,j+1/2, wi,j1/2.
В выражения для них входят производные, которые аппроксимируются ценx x тральными разностями. Например, приведем пример для wi+1/2,j, wi1/2,j wi+1/2,j = wi1/2,j = Аналогично для wi,j+1/2, wi,j1/ i,j+1/2 hi+1/2,j+1/2 ui+1/2,j+1/2 ui+1/2,j+1/2 hi1/2,j+1/2 ui1/2,j+1/2 ui1/2,j+1/ wi,j+1/2 = i,j1/2 hi+1/2,j1/2 ui+1/2,j1/2 ui+1/2,j1/2 hi1/2,j1/2 ui1/2,j1/2 ui1/2,j1/ wi,j1/2 = Используем интегро-интерполяционный метод для аппроксимации уравнения (1.10) Когда будем записывать аппроксимацию для ij, то будем использовать вид (1.9). Чтобы не загромождать формулу, вместо переменных wx, wy и R будем использовать wsx, wsy и Rs. Суть от небольшой перемены в обозначениях не меняется.
i+1/2,j = ui+1/2,j wsi+1/2,j + Rsi+1/2,j i1/2,j = ui1/2,j wsi1/2,j + Rsi1/2,j i,j+1/2 = ui,j+1/2 wsi,j+1/2 + Rsi,j+1/2 i,j1/2 = ui,j1/2 wsi,j1/2 + Rsi,j1/ Далее аппроксимируем выражения, входящие в wsx, wsy и Rs. Все формулы можно получить, глядя на шаблон из рис. 3.1. По сути выражения, которые здесь приводятся, являются промежуточными выкладками для нахождения h, ux, uy на следующем временном слое. Поэтому все выписанное выражения просто облегчают программную реализацию, иначе окончательные выражения получились бы еще более громоздкими. Рассмотрим величину wsx Теперь также аккуратно аппроксимируем выражение для wsy i+1/2,j i1/2,j i,j+1/ i,j1/ Осталось разобраться с последней величиной Rs Rsi+1/2,j = gi+1/2,j ux 0.5hi+1/2,j+1/2 hi+1/2,j+1/2 0.5hi+1/2,j1/2 hi+1/2,j1/ Rsi1/2,j = gi1/2,j ux 0.5hi1/2,j+1/2 hi1/2,j+1/2 0.5hi1/2,j1/2 hi1/2,j1/ Rsi,j+1/2 = gi,j+1/2 ux Rsi,j1/2 = gi,j1/2 ux На этом этапе для рассматриваемой точки можно найти ux и uy Отдельно выпишем выражение для h (Пока полагаем, что hi,j = hi,j = hi,j ) Все выписанные выражения являются оптимальной формой для записи разностных уравнений. На самом деле, описанный способ расчета можно оптимизировать. Это распространенная практика, когда сначала вычисляют все потоки на ребрах прямоугольников, а затем находят значения величин на следующем временном слое. Действительно, число вычислений сократится.
Здесь, рассматривая точку (i, j), находили значения ji1/2,j и ji1/2,j. Но величина ji1/2,j с точностью до знака совпадает с jix+1/2,j, где i = i 1. И приведенный способ приводит к тому, что некоторые вычисления проводятся по два раза, когда достаточно одного. Это один из известных способов оптимизации, который сократит время работы, но увеличит память, которую для работы будет использовать программа, так как нужно где-то хранить промежуточные величины.
3.2. Условие покоящейся жидкости для неровного дна Проведем дополнительные модификации численного алгоритма, чтобы для него выполнялось условие покоящейся жидкости. Поставим условие, что жидкость в любой момент времени остается в покое, если для начальных условий выполняется fx = fy = 0, ux = uy = 0. К этому нужно добавить обязательное условие h(x, y) + b(x, y) = const. Используя эти условия, преобразуем выражения (3.2) и (3.3) к следующему виду После этих упрощений легко заметить, каким образом нужно поменять hi,j, hi,j, чтобы выражения (3.2), (3.3) были равны нулю. В результате мы полуx) (y) чим, что hi,j, hi,j должны вычисляться как средние значения двух соседних полуцелых точек на соответствующих осях x и y. Соответствующий результат запишем для выражения (3.4) В таком виде условие покоящейся жидкости для неровного дна будет выполняться автоматически. В дальнейшем условие для покоящейся жидкости будет расширено на численный алгоритм для неструктурированных сеток.
3.3. Задача о разрушении несимметричной дамбы Построение численного алгоритма для расчета двумерных течений было приведено в предыдущем разделе. В нем также был явно указан вид аппроксимированных величин. Возможности двумерного алгоритма иллюстрируются примером расчета нестационарного течения, возникающего при разрушении несимметричной дамбы. Эту задачу также используют для тестирования численного алгоритма, например, в [41], [60] и в приведенных в этих работах ссылках. В этих расчетах внешние силы не учитываются, fx = fy = 0, и дно водоема полагается плоским b(x, y) = 0.
В соответствии с [41] рассматривается задача о течении, возникающем при мгновенном разрушении дамбы, разделяющей два бассейна с водой. Высота уровня воды в левом бассейне составляет 10м, в правом 5м. Длина разрыва равна 75м, начало разрыва расположено в точке с координатой y = 95м.
(рис. 3.2). Толщина стенки дамбы равна 10м, и ее левая сторона расположена в точке с координатой x = 95м. На всех границах дамбы ставятся отражающие граничные условия. g = 9.8м/с2.
Рис. 3.2. Разрушении несимметричной плотины. Толщина слоя жидкости и линии тока.
Высота уровня жидкости и распределение линий тока (касательные к скорости), а также изолинии числа Фруда представлены на рис. 3.2 и 3.3 для момента времени t = 7.2 с. и равномерной пространственной сетки с шагами x = y = 1м. Вычисления проведены для параметра регуляризации = 0.2 и числа Куранта = 0.2. На рисунке видны характерные особенности возникающего к указанному моменту времени течения, а именно, сглаженный немонотонный профиль слева от разрыва и резкий но монотонный профиль в правом бассейне, а также отражение волны от верхней стенки бассейна.
На рис. 3.4 продемонстрированы одномерные распределения уровня Рис. 3.3. Разрушении несимметричной плотины. Распределение числа Фруда.
жидкости h(x) и числа Фруда F r(x), соответственно, вдоль линии y = 160м для последовательности сеток с шагами x = y= 1 (линия 1), 2 (линия 2) и 4 (линия 3)м. Из этих рисунков следует сходимость численного решения при сгущении сетки к эталонному решению, в качестве которого выбрано решение этой задачи, приведенное в [41]. В указанной работе задача о дамбе решается на неструктурированной пространственной сетке с помощью двух численных алгоритмов высокого порядка точности, результаты которых оказываются близкими. Результаты расчетов на основе алгоритма с регуляризацией для x = y = 1м хорошо совпадают с решением [41], полученным на сетке с шагами 2м.
Двумерный алгоритм был адаптирован для проведения расчетов на многопроцессорной вычислительной системе. Это возможно благодаря тому, Water height level Рис. 3.4. Разрушении несимметричной плотины. Сходимость по сетке. Момент времени t = 7.2с, срез по линни y = 160 м.
что схема является явной по времени, и тому, что расчеты проводили на прямоугольной сетке. В этом случае область делится на несколько одинаковых частей. На рис. 3.5 область разделена на четыре части. При реализации создается массив S(i, j), в котором записана информация о типе точке: является ли точка внутренней или граничной, выпадает ли рассматриваемая точка из области. Каждый тип обозначен собственным номером. Поэтому здесь мы используем прямоугольные матрицы для хранения величин, даже если граница имеет сложную форму. На рис. 3.5 также видно, что области граничат по вертикальным линиям, на которые ложатся точки, принадлежащие одному столбцу матрицы. Каждый отдельный процессор обсчитывает свою часть области. Но на каждом новом временном слое они должны обменяться значениями, которые принадлежат граничной линии. А это составляет один столбец. Стоит заметить, что данная задача и данный метод просто созданы для распараллеливания. Играет свою роль и то, что используются прямоугольные сетки. Если вести расчет на неструктурированных сетках, то возникает проблема разделения области на части и нахождения точек, которые являются граничными для соседних областей.
«