CC BY
21
Вычислительные технологии
Том 24, № 1, 2019
Оценка численной диффузии метода конечных объемов при моделировании поверхностных волн*
Е.С. Тятюшкинл1,2, А. С. Козелков1,2, А. А. Куркин2^, В. В. Курулин1, В. Р. Ефремов3, Д. А. Уткин1
1 Российский федеральный ядерный центр — Всероссийский научно-исследовательский институт экспериментальной физики, Саров, Россия
2Нижегородский государственный технический университет им. Р.Е. Алексеева, Нижний Новгород, Россия
3Конструкторское бюро приборостроения им. акад. А.Г. Шипунова, Тула, Россия ^Контактный e-mail: aakurkin@gmail.com
Обсуждается применение метода конечных объемов при решении уравнений Навье — Стокса для моделирования поверхностных волн. Сформулирована задача о распространении поверхностных волн, которая используется для оценки численной диффузии в решении уравнений Навье — Стокса. Предлагается методика оценки численной диффузии, выражаемой коэффициентом уменьшения амплитуды волны при прохождении ею одной своей длины (коэффициентом затухания). Дана оценка размеров сетки и шага по времени, выраженных в безразмерных величинах относительно параметров волны, необходимых для обеспечения приемлемого значения коэффициента затухания. Показана степень влияния каждого из сеточных параметров на увеличение коэффициента затухания.
Ключевые слова: уравнения Навье — Стокса, поверхностная волна, volume of fluid, метод конечных объемов, разностная схема, сеточное разрешение.
Библиографическая ссылка: Тятюшкина Е.С., Козелков А.С., Куркин А.А., Курулин В.В., Ефремов В.Р., Уткин Д.А. Оценка численной диффузии метода конечных объемов при моделировании поверхностных волн // Вычислительные технологии. 2019. Т. 24, № 1. С. 106-119. DOI: 10.25743/ICT.2019.24.1.008.
Введение
Применение методов численного моделирования, основанных на решении трехмерных уравнений Навье — Стокса [1-3], для исследования распространения волн на поверхности водной среды, предполагает построение сеточной модели, содержащей счетные узлы во всей толще водной среды [4, 5]. Любая численная модель должна иметь разумные ограничения на используемый шаг по времени. Сеточное разрешение и шаг по времени напрямую определяют точность полученных результатов и устойчивость итерационного процесса. Низкое сеточное разрешение приводит к недостаточной детализации полей скорости, давления, объемной доли жидкости, что повышает численную диффузию метода и, в конечном счете, приводит к занижению амплитуд колебаний водной среды
* Title translation and abstract in English can be found on page 119.
© ИВТ СО РАН, 2019.
в зоне расположения мареографов. Большой шаг по времени также вызывает "размытие" решения и существенно снижает его устойчивость, особенно при использовании схем, сжимающих фронт раздела сред (МСГСБАМ, НШС) [6, 7]. Оценка необходимого сеточного разрешения при использовании уравнений Навье — Стокса наиболее часто встречается при рассмотрении турбулентных течений однородной жидкости [8-10].
В указанных работах для определения сеточных размеров используются критерии, сформулированные относительно линейного масштаба пограничного слоя, а в качестве критерия для определения шага по времени применяется критерий Куранта. Показано, что использование недостаточного сеточного разрешения и схем с высокой численной диффузией приводит к значительным погрешностям в итоговых результатах. Вопрос сеточного разрешения и параметров численных схем актуален также для моделирования распространения поверхностных волн.
В настоящей статье приведены результаты оценки приемлемых размеров сетки и шага по времени, выраженных в безразмерных величинах относительно параметров волны, необходимых для обеспечения точности решения, достаточной для геофизических приложений. Представленная методика расчета трехмерных многофазных течений со свободной поверхностью основывается на численном решении уравнений Навье — Стокса в односкоростном приближении на базе полностью неявной связи скорости и давления. Методика расчетов детально изложена в [4, 11], она основана на неявной связи уравнения неразрывности и уравнений сохранения количества движения за счет слагаемых градиента давления и массового потока. В указанных работах приведены основные формулы дискретизации уравнений, вид коэффициентов, которые суммируются в общую связанную матрицу, и основные этапы вычислительной процедуры. Верификация методики осуществляется на задачах, имеющих экспериментальные данные (задача об обрушении плотины, гидравлический прыжок и падение параллелепипеда в воду).
Для определения приемлемых размеров сетки и шага по времени с использованием представленной методики введены параметры, на основе которых предложена методика оценки численной диффузии расчетного метода, выражаемая коэффициентом уменьшения амплитуды волны при прохождении ею одной своей длины (далее — коэффициентом затухания). С использованием описанной методики получены результаты оценки сеточного разрешения в горизонтальном и вертикальном направлении, оценка шага по времени, а также оценка влияния схемы дискретизации конвективного слагаемого.
1. Базовая задача
Для оценки численной диффузии метода при численном моделировании распространения поверхностных волн сформулируем базовую задачу. Рассмотрим распространение синусоидальных волн вдоль протяженного бассейна (рис. 1). Волны генерируются у левого края бассейна и движутся до правой границы, где свободно, без отражения, покидают расчетную область.
Введем обозначения: А — амплитуда волны; Л — длина волны; к = 2ж/X — волновое число; к — глубина бассейна; ш = 2п/Т — круговая частота, где Т — период, связанный с длиной волны дисперсионным соотношением [9, 10]; д — ускорение свободного падения:
Рис. 1. Общий вид задачи
Для генерации волн на левой границе используются следующие граничные условия:
ц = А
их = Аи $т(иг)е-ку, иу = Аш со$(^)е-ку,
где ц — уровень воды на левой границе; их, иу — компоненты скорости на левой границе. Такие граничные условия соответствуют решению задачи о распространении синусоидальных волн при А << к [12], что обеспечивает генерацию пакета синусоидальных волн с наименьшими искажениями на входной границе.
Будем рассматривать невязкую жидкость, по поверхности которой распространяются волны, сохраняя свою амплитуду1. В ходе численного моделирования задачи при конечном сеточном и временном разрешении амплитуда волн будет постепенно уменьшаться по причине влияния численной диффузии на решение задачи. По степени уменьшения амплитуды волн при прохождении ими определенного расстояния можно судить о величине численной диффузии, а следовательно, о качестве сеточного разрешения и шаге по времени.
На рис. 2 показана типичная расчетная сетка со сгущением к поверхности раздела воздух —вода, которая строится для решения задач данного класса. Сетка построена средствами препостпроцессора пакета программ ЛОГОС [13]. Подробное описание алгоритма генерации сеточных моделей для распространения поверхностных волн в произвольных бассейнах с учетом батиметрических данных приводится в [4, 5].
Вдоль горизонтальной оси (Ах) сетка имеет равномерное разбиение. Наименьший вертикальный размер Ау ячейки — у свободной поверхности, вверх и вниз от которой вертикальный размер ячеек увеличивается с коэффициентом 7 = Ауп/Ауп-\, пока не достигнет заданного размера Аут вблизи дна бассейна и верхней границы расчетной области.
При аппроксимации решения уравнений методом конечных объемов на расчетной сетке такого типа (особенности такой аппроксимации подробно изложены в [11, 14]) можно выделить следующие ключевые параметры, которые определяют численную диффузию расчетного метода:
1) вертикальный размер ячеек вблизи свободной поверхности, Ау;
2) вертикальный размер ячеек вблизи дна бассейна, Аут;
3) горизонтальный размер ячеек, Ах.
1 Для чистоты численного эксперимента физическая вязкость жидкости убирается сознательно, поскольку она вносит свой вклад в процесс затухания.
Первый параметр определяет количество ячеек, приходящееся на амплитуду волны А (А/Ду), второй параметр — количество ячеек, приходящееся на глубину бассейна к, третий — количество ячеек, приходящееся на длину волны Л. Эти параметры должны быть дополнены величиной используемого шага по времени Д£. Большой шаг по времени порождает большую численную диффузию [14], поэтому в расчетах он должен быть ограничен сверху значением, при котором численная диффузия имеет еще приемлемую величину. Наряду с этим при проведении практических расчетов шаг по времени ограничивается сверху дополнительным критерием, который зависит от устойчивости численного метода. Оценка того, какой из этих критериев для шага по времени наступает раньше в рассматриваемых задачах, — одна из целей данного исследования.
Методика исследования влияния вышеперечисленных параметров на численную диффузию расчетного метода заключается в следующем. Для определения влияния на решение одного из параметров проводится ряд расчетов с варьированием значения этого параметра. Значения остальных параметров являются фиксированными и избыточно малыми, так что гарантированно не вносят большого вклада в численную диффузию, влияющую на решение. Влияние на решение выражается в коэффициенте уменьшения амплитуды волны при прохождении ею одной своей длины:
=1 - ^ (1)
(п — номер волны, начиная от левого края бассейна). Будем рассматривать среднюю величину 8 =< 8П >. Очевидно, что чем мельче сетка и меньше шаг по времени, тем коэффициент затухания 8 ближе к нулю. Значение коэффициента затухания будем считать приемлемым, если вследствие влияния численной диффузии амплитуда волны при прохождении одной собственной длины
8П < 0.5 %. (2)
Используем далее этот факт для оценки всех четырех критериев, указанных выше.
2. Оценка сеточного разрешения в горизонтальном направлении
Для определения приемлемого горизонтального размера ячеек проведем серию расчетов с одними и теми же параметрами волны: А = 2.5 • 10-2 м, Л = 20 м, к = 5 м, длина бассейна Ь = 2 • 103 м (см. рис. 1). В каждом расчете остаются постоянными шаг по времени и вертикальный размер ячеек, горизонтальный размер ячеек изменяется. В табл. 1 приведены параметры расчетных сеток и шаг по времени (столбец означает количество расчетных ячеек, приходящееся на одну длину волны), а также итоговые значения 8, полученные по описанной методике, для всех расчетных случаев.
Во всех девяти случаях расчетная сетка имеет постоянные значения Ду, , причем такие, что по отношению к основному параметру Дх их влияние на численную диффузию пренебрежимо мало2. Это касается и шага по времени Д£, который взят малым по отношению к периоду волны. Размер ячеек вдоль горизонтальной оси изменялся от 0.25 до 2 м, что означает 80 и 10 расчетных ячеек на одну длину волны соответственно.
В расчетах на длину бассейна Ь приходятся 100 длин волн. Для получения статистически установившейся картины колебаний расчет производился до момента времени т = 200 Т. По окончании расчета строился график положения свободной поверхности вдоль горизонтальной оси. Положение свободной поверхности определялось путем построения изоповерхности по значению объемной доли жидкости, равной 0.5.
В качестве примера проанализируем расчет 1 (табл. 1). На рис. 3 представлены график свободной поверхности п(х) и график огибающей о(х), которая с учетом (1) может быть представлена через параметр 8:
а(х) = А(1 - 8)х/х. (3)
Здесь огибающая а(х) наилучшим образом повторяет падение амплитуды волн вдоль длины бассейна при значении 8 = 0.15%, это и будет являться оценкой значения параметра 8 для фиксированной величины Дх. Значение 8 = 0.15% заметно меньше допустимого, поскольку получено для самого мелкого шага Дх = 0.25, что соответствует разрешению в 80 ячеек на одну длину волны.
Таблица 1. Параметры и результаты расчетов первого этапа
Номер расчета Дх, м Дх/Х Ду/А Ду-и, /к ДЪ/Т 5, %
1 0.25 0.0125 80 0.1 0.1 0.007 0.15
2 0.33 0.0165 60 0.1 0.1 0.007 0.18
3 0.50 0.0250 40 0.1 0.1 0.007 0.25
4 0.75 0.0375 27 0.1 0.1 0.007 0.35
5 1.00 0.0500 20 0.1 0.1 0.007 0.8
6 1.25 0.0625 16 0.1 0.1 0.007 1.1
7 1.50 0.0750 13 0.1 0.1 0.007 1.5
8 1.75 0.0875 11 0.1 0.1 0.007 1.7
9 2.00 0.1000 10 0.1 0.1 0.007 2.2
2Данное утверждение сделано исходя из серии предварительных численных экспериментов, которые показали, что такие значения параметров минимально влияют на общую численную диффузию метода.
О 500 1000 1500 2000
Рис. 3. Графики свободной поверхности п(х) — сплошные линии и огибающей а(х) — штрихи
Результаты показывают, что при увеличении значения Ах и, следовательно, уменьшении Уд величина 5 монотонно возрастает, что свидетельствует о возрастании численной диффузии. Расчетный случай с наибольшим значением Ах, при котором 5 еще удовлетворяет условию (2), выделен в таблице курсивом. На рис. 4 приведены зависимость 5(Ах/Х) (точки соответствуют значениям из табл. 1), аппроксимирующая функция 5н(Ах/Х) и прямая 5 = 0.5, отвечающая верхней границе условия (2).
Аппроксимирующая функция 5ь(Ах/Х) взята в следующем виде (к = 250):
= кх2.
Используя функцию 6ь(Ах/Х), можно определить предельное значение горизонтального разрешения:
Ах/Х < 0.045, Уа > 22.
При таких значениях горизонтального разрешения численная диффузия не превышает 5 = 0.5. Этот критерий получен для длины волны Л = 20 м, но неизвестно, будет ли он универсален и справедлив для других длин волн. Для проверки этого увеличим в два раза длину волны: А1 = 2Л = 40 м. На рис. 5 представлены график свободной поверхности п(х) для проведенного расчета и график огибающей а(х). Коэффициент затухания
0 0.02 0.04 0.06 0.08 0.1 Ах/Л
Рис. 4. Графики ¿(Ах/А) — значки ■, ¿^(Ах/А) — сплошная кривая и 5 = 0.5 — пунктир
Рис. 5. Графики свободной поверхности п(х) — сплошные линии и огибающей а(х) — штрихи (расчет 3)
составил 5 = 0.14%, что близко к значению, которое получено ранее для другой длины волны, но при таких же значениях безразмерных параметров: Ах/Х = 0.0125, = 80 (расчетный случай 1, табл. 1).
3. Оценка сеточного разрешения в вертикальном направлении
Определение приемлемого вертикального размера ячеек вблизи свободной поверхности Ау и Ау,ю проводится с использованием той же методики. Параметры волны и бассейна взяты такие же, как и в предыдущем случае. Вначале при неизменных шаге по времени А£, горизонтальном размере Ах и вертикальном размере у стенки Ауы варьировалось значение параметра Ау. В табл. 2 приведены параметры расчетных сеток и результаты расчетов по значению параметра 5 = А/Ау). Как видно из таблицы, в расчетах Ау меняется от 0.002 до 0.04 м. В самой мелкой расчетной сетке на амплитуду волны приходятся более 12 ячеек, в самой крупной сетке на двойную амплитуду волны — чуть более одной расчетной ячейки.
С ростом Ау наблюдается уменьшение числа Куранта при постоянном шаге по времени. Это объясняется тем, что в расчетной сетке Ау — минимальный линейный размер ячеек, которые располагаются вблизи свободной поверхности, где волны распространя-
Таблица 2. Параметры анализа вертикального разрешения
Номер расчета Ау, м Ау/А Ах/Х Ауш/Ь АЬ, с СРЬ(Аг)
1 0.002 0.08 12.5 0.025 0.1 0.025 0.50 0.25
2 0.004 0.16 6.3 0.025 0.1 0.025 0.27 0.24
3 0.008 0.32 3.1 0.025 0.1 0.025 0.14 0.23
4 0.010 0.40 2.5 0.025 0.1 0.025 0.10 0.25
5 0.012 0.48 2.1 0.025 0.1 0.025 0.07 0.33
6 0.016 0.64 1.6 0.025 0.1 0.025 0.06 0.38
7 0.020 0.80 1.2 0.025 0.1 0.025 0.05 0.40
8 0.040 1.60 0.6 0.025 0.1 0.025 0.045 0.42
ются с наибольшими скоростями. Поэтому Ду значительно влияет на максимальное число Куранта, которое наблюдается в расчете.
Результаты показывают, что значение Ду существенно меньше влияет на численную диффузию, чем горизонтальный размер Дх. До значений Ма = 0.21 величина 8 держится на одном и том же уровне, диктуемом, по всей видимости, другими расчетными факторами. Отмечено незначительное колебание 8 при увеличении Ду, которое наблюдается в расчетах 2-4 (табл. 2), что, вероятно, связано с изменением числа Куранта в расчетах.
Чувствительность численной диффузии к значению Ду наблюдается лишь после того, как на амплитуду волны начинает приходиться менее двух ячеек расчетной сетки (расчеты 6-8, табл. 2). С учетом этого для вертикального сеточного размера вблизи свободной поверхности выберем следующий критерий:
Ду/А < 0.5, МА > 2.
Такой вертикальный сеточный размер вблизи свободной поверхности позволяет сохранять численную диффузию метода в пределах выбранного допустимого значения.
Для оценки приемлемого вертикального размера ячеек вблизи дна бассейна проведены расчеты при неизменных шаге по времени Д£, горизонтальном размере Дх и вертикальном размере у свободной поверхности Ду. Размер варьировался путем изменения коэффициента роста вертикального размера ячеек 7 = Дуп/Дуп-\. В табл. 3 приведены параметры расчетных сеток и результаты расчетов параметра 8 (^ — количество расчетных ячеек, приходящееся на всю глубину бассейна). Как видно из таблицы, коэффициент скорости роста 7 в расчетах варьировался от значений 1.1 до крайне больших — 3.36, которые обычно не применяются на практике. При этом количество ячеек, приходящееся на глубину бассейна, изменялось от 58 до 9. Такие сильные изменения в сеточном разрешении водной толщи бассейна не привели к заметным изменениям в численной диффузии. Коэффициент затухания для всех расчетов составил около 8 = 0.25, и это значение диктовалось горизонтальным размером ячеек Дх/Х = 0.025.
Третья серия расчетов показала, что при выборе сеточного размера критерий, связанный с численной диффузией, не являлся существенным. При решении практических задач главным критерием здесь будет достаточное количество ячеек для описания рельефа дна, а также сохранение в разумных пределах скорости роста размеров ячеек, который, как правило, не превышает 7 = 1.2.
Таблица 3. Параметры расчетов третьего этапа
Номер расчета ДуУ], м 7 Дх/Х Дут/Ъ Д£, с СТЬД) 5,%
1 0.50 1.10 58 0.025 0.1 0.025 0.5 0.25
2 0.66 1.15 43 0.025 0.1 0.025 0.5 0.24
3 0.83 1.20 35 0.025 0.1 0.025 0.5 0.25
4 1.00 1.25 30 0.025 0.1 0.025 0.5 0.25
5 1.35 1.35 23 0.025 0.1 0.025 0.5 0.25
6 1.70 1.50 18 0.025 0.1 0.025 0.5 0.26
7 2.37 1.90 14 0.025 0.1 0.025 0.5 0.24
8 3.26 2.87 9 0.025 0.1 0.025 0.5 0.26
4. Оценка значения шага по времени
Учитывая полученные ранее результаты, выберем оптимальное сеточное разрешение для задачи: Ах/Х = 0.375 в соответствии с табл. 1, Ау/А = 0.48 в соответствии с табл. 2 и 7 = 1.2. На расчетной сетке с такими параметрами оценим влияние на численную диффузию шага по времени АЬ путем проведения серии расчетов с его варьированием. В табл. 4 представлены параметры расчетов четвертого этапа и итоговые значения коэффициента затухания 8.
Число Куранта в задаче определим как максимальное значение выражения для всех ячеек расчетной сетки:
А1 £
СРЬ(Аг) = —-.
Здесь суммирование проводится по граням ячейки с положительным объемным потоком Рк+; V — объем ячейки.
Результаты показывают сильную зависимость коэффициента затухания 8 от шага по времени А£. Приемлемые результаты согласно (3) получаются лишь при АЬ < 0.037, что соответствует АЬ/Т = 0.01, т.е. примерно Ыт = 100 расчетным шагам на один период волны. При таком шаге по времени число Куранта составило СРЬ(АЬ) = 0.18. В последнем расчетном случае 6 с максимальным шагом по времени на начальном этапе генерирования волн возникают численные осцилляции, которые приводят к расхождению решения.
Пятая серия расчетов показала граничные значения АЬ/Т и СРЬ(АЬ), однако не ответила на вопрос, какое именно из этих чисел определяет уровень численной диффузии при расчете распространения волны. Первый параметр зависит только от физических свойств волны, второй — от физических свойств волны и величины ячеек в расчетной сетке. Для ответа на вопрос повторим расчеты на сетке с двукратно уменьшенным вертикальным разбиением вблизи свободной поверхности: Ау/А = 0.24. Это сохранит на прежнем уровне соотношение АЬ/Т, но в два раза увеличит СРЬ(АЬ). Параметры пятой серии расчетов представлены в табл. 5. Как видно, двукратное увеличение вертикального разбиения вызывает двукратное увеличение числа Куранта. Однако рост числа Куранта приводит не к увеличению численной диффузии, а к ее уменьшению примерно на четверть. Это связано лишь с уменьшением сеточного разбиения, что хорошо согласуется с результатами табл. 2, где двукратное уменьшение Ау также приводит к уменьшению на четверть численной диффузии.
Таблица 4. Параметры расчетов четвертого этапа
Номер расчета А* Аг/т сть(Аг)
1 0.0125 0.003 300 0.05 0.30 +0.05
2 0.025 0.007 150 0.11 0.32
3 0.037 0.010 101 0.18 0.45
4 0.050 0.013 75 0.22 0.75
5 0.060 0.016 62 0.29 0.90
6 0.075 0.020 50 ^0.35 —
Таблица 5. Параметры расчетов пятого этапа
Номер расчета At At/T NT CFL(At)
1 0.0125 0.003 300 0.1 0.25
2 0.025 0.007 150 0.22 0.27
3 0.037 0.010 101 0.36 0.37
4 0.050 0.013 75 0.44 0.58
5 0.060 0.016 62 0.60 —
6 0.075 0.020 50 ^0.70 —
Указанные факты означают, что критерием шага по времени для обеспечения приемлемого уровня численной диффузии является именно соотношение At/T либо количество шагов, приходящееся на период волны Nt, а не число Куранта. Число Куранта в задачах такого класса определяет только меру устойчивости расчетного алгоритма, о чем свидетельствует получение расходящегося решения уже не в одном расчетном случае (табл. 4, расчет 6), а в двух (табл. 5, расчеты 5 и 6).
Из представленных результатов следует, что наибольшее влияние на численную диффузию оказывают горизонтальное сеточное разбиение по отношению к длине волны Ax/X и шаг по времени по отношению к периоду волны At/T. Вертикальное сеточное разбиение вблизи свободной поверхности по отношению к амплитуде волн не оказывает большого влияния на численную диффузию при Ay < 0.5А. Количество расчетных точек, приходящееся на глубину бассейна, также не оказывает заметного влияния на численную диффузию и, скорее всего, будет диктоваться другими факторами (корректным учетом рельефа дна, трения и пр.).
Выше были установлены количественные критерии сеточного разрешения для коэффициента затухания 5 < 0.5 %. С использованием аналогичной методики получены оценки необходимого сеточного разрешения и для других значений коэффициента затухания, которые приведены в табл. 6.
5. Оценка влияния схемы дискретизации конвективного слагаемого
В [14, 15] показано, что численная диффузия метода напрямую зависит от схемы дискретизации конвективного слагаемого. Так, использование схемы GAMMA вместо схемы UD позволяет в несколько раз уменьшить численную диффузию при расчете турбулентных течений в однородной несжимаемой жидкости [15].
В табл. 6 приведены критерии необходимого сеточного разрешения, полученные для схемы UD [8, 12]. С применением аналогичной методики получены оценки необходимого
Таблица 6. Результаты расчетов для различных значений 5
Номер расчета 5,% Ax/X Nx At/T NT Ay/A Na
1 0.5 0.045 22 0.010 100 0.5 2
2 1.0 0.075 13 0.015 67 0.5 2
3 2.0 0.090 11 0.017 59 0.5 2
Таблица 7. Результаты расчетов для различных схем дискретизации конвективного слагаемого
Схема Ax/X Nx At/T NT Ay/A na
UD 0.045 22 0.010 100 0.5 2
LUD 0.054 19 0.010 100 0.5 2
QUICK 0.063 16 0.010 100 0.5 2
BCD 0.085 12 0.010 100 0.5 2
сеточного разрешения для схем LUD, QUICK [8, 12] и BCD [11] при коэффициенте затухания 5 = 0.5% (3). Результаты расчетов представлены в табл. 7.
При использовании различных схем дискретизации конвективного слагаемого численная диффузия, связанная с горизонтальным разрешением Ах/Х, заметно изменяется. Для схемы LUD горизонтальный размер ячеек может быть увеличен в 1.2 раза, для схемы QUICK — в 1.4 раза, для схемы BCD — в 1.9 раза. Для трехмерных задач такое увеличение горизонтального размера ячеек приводит к существенному снижению требуемого числа ячеек; в частности, для схемы BCD оно относительно схемы UD может быть снижено в 3.5 раза. Чувствительность численной диффузии к вертикальному размеру ячеек вблизи свободной поверхности и шагу по времени остается прежней.
Оценки по используемым схемам, сеточному разрешению и выбору шага по времени, представленные в настоящей работе, являются ценной информацией для моделирования волн цунами космогенного и оползневого происхождения на основе уравнений Навье — Стокса.
Таким образом, представлены результаты оценки сеточных размеров, шага по времени, а также влияния конвективного слагаемого для применения метода конечных объемов при моделировании распространения поверхностных волн. Введены параметры, на основе которых предложена методика оценки численной диффузии расчетного метода, выражаемая коэффициентом уменьшения амплитуды волны при прохождении ею одной своей длины.
Благодарности. Работа выполнена при финансовой поддержке гранта Президента Российской Федерации по государственной поддержке научных исследований ведущих научных школ Российской Федерации НШ-2685.2018.5, молодых российских ученых-докторов наук МД-4874.2018.9, а также при финансовой поддержке РФФИ — проекты N 16-01-00267 и 17-05-00067.
Список литературы / References
[1] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1973. 847 с.
Loytsyanskiy, L.G. Mechanics of liquid and gas. Moscow: Nauka, 1973. 847 p. (In Russ.)
[2] Ландау Л.Д., Лифшиц Е.М. Гидродинамика. М.: Наука, 1986. 736 с.
Landau, L.D., Lifshits, E.M. Hydrodynamics. Moscow: Nauka, 1986. 736 p. (In Russ.)
[3] Ferziger, J.H., Peric, M. Computational method for fluid dynamics. New York: SpringerVerlag, 2002. 423 p.
[4] Козелков А.С. Методика численного моделирования цунами оползневого типа на основе уравнений Навье — Стокса // Вычисл. механика сплошных сред. 2016. Т. 9, № 2. C. 218-236.
Kozelkov, A.S. Numerical technology for landslide tsunami simulations based on Navier — Stokes equations // Comput. Continuum Mechanics. 2016. Vol. 9, No. 2. P. 218-236. (In Russ.)
[5] Kozelkov, A.S., Kurkin, A.A., Pelinovsky, E.N. et al. Landslide-type tsunami modelling based on the Navier —Stokes equations // Sci. of Tsunami Hazards. 2016. Vol. 35, No. 3. P. 106-144.
[6] Ubbink, O. Numerical prediction of two fluid systems with sharp interfaces: PhD Thesis. London, Imperial College of Science, Technology & Medicine, 1997. 138 p.
[7] Храбрый А.И. Численное моделирование нестационарных турбулентных течений жидкости со свободной поверхностью: Дис. ... канд. физ.-мат. наук. СПб., 2014. 154 с. Khrabryy, A.I. Numerical simulation of unsteady turbulent fluid flow with free surface: PhD Thesis. Sankt-Peterburg, 2014. 154 p. (In Russ.)
[8] Козелков А.С., Курулин В.В., Пучкова О.Л., Тятюшкина Е.С. Моделирование турбулентных течений вязкой несжимаемой жидкости на неструктурированных сетках с использованием модели отсоединенных вихрей // Матем. моделирование. 2014. Т. 26, № 8. С. 81-96.
Kozelkov, A.S., Kurulin, V.V., Tyatyushkina, E.S., Puchkova, O.L. Application of the detached eddy simulation model for viscous incompressible turbulent flow simulations on unstructured grids // Mathematical Models and Computer Simulations. 2014. Vol. 26, No. 8. P. 81-96. (In Russ.)
[9] Kozelkov, A., Kurulin, V., Emelyanov, V. et al. Comparison of convective flux discretization schemes in detached-eddy simulation of turbulent flows on unstructured meshes // J. of Sci. Computing. 2016. Vol. 67. P. 176-191.
[10] Козелков А.С., Ш^агалиев Р.М., Курулин В.В. и др. Исследование потенциала суперкомпьютеров для масштабируемого численного моделирования задач гидродинамики в индустриальных приложениях // Вычисл. математика и матем. физика. 2016. Т. 56, № 8. С. 1524-1535.
Kozelkov, A.S., Shagaliev, R.M., Kurulin, V.V. et al. Investigation of supercomputer capabilities for the scalable numerical simulation of computational fluid dynamics problems in industrial applications // Comput. Mathematics and Math. Physics. 2016. Vol. 56, No. 8. P. 1506-1516.
[11] Козелков А.С., Мелешкина Д.П., Куркин А.А. и др. Полностью неявный метод решения уравнений Навье — Стокса для расчета многофазных течений со свободной поверхностью // Вычисл. технологии. 2016. Т. 21, № 5. С. 54-76.
Kozelkov, A.S., Meleshkina, D.P., Kurkin, A.A. et al. Fully implicit method for solution of Navier-Stokes equations for simulation of multiphase flows with free surface // Comput. Technologies. 2016. Vol. 21, No. 5. P. 54-76. (In Russ.)
[12] Fenton, J.D. A fifth-order Stokes theory for steady waves // Coastal and Ocean Eng. 1985. Vol. 111, No. 2. P. 216-234.
[13] Сафронов А.В., Дерюгин Ю.Н., Жучков Р.Н. и др. Результаты валидации многофункционального пакета программ ЛОГОС при решении задач аэрогазодинамики старта и полета ракет-носителей // Матем. моделирование. 2014. Т. 26, № 9. С. 83-95. Safronov, A.V., Deryugin, Yu.N., Zhuchkov, R.N. et al. Validation results for the LOGOS multifunction software package in solving problems of aerodynamics and gas dynamics for the lift-off and injection of launch vehicles. // Math. Models and Comput. Simulations. 2015. Vol. 7, No. 2. P. 144-153.
[14] Weinman, K.A., Valentino, M. Comparison of hybrid RANS-LES calculations within the framework of compressible and incompressible unstructured solvers / Eds Sh.-H. Peng, P. Doerffer, W. Haase. Progress in Hybrid RANS-LES Modelling. Springer, 2010. P. 329-338.
[15] Козелков А.С., Крутикова О.Л., Курулин В.В. и др. Применение численных схем с выделением пограничного слоя для расчета турбулентных течений с использованием вихреразрешающих подходов на неструктурированных расчетных сетках // Вычисл. математика и матем. физика. 2017. Т. 57, № 6. С. 1048-1060.
Kozelkov, A.S., Krutyakova, O.L., Kurulin, V.V. et al. Application of numerical schemes with singling out the boundary layer for the computation of turbulent flows using eddy-resolving approaches on unstructured grids // Comput. Mathematics and Math. Physics. 2017. Vol. 57, No. 6. P. 1036-1047.
Поступила в 'редакцию 17 декабря 2018 г.
Evaluation of numerical diffusion of the finite volume method when modelling surface waves
Tyatyushkina, Elena S.1'2, Kozelkov, Andrey S.1'2, Kurkin, Andrey A.2'*, Kurulin, Vadim V.1, Efremov, Valentin R.3, Utkin, Dmitriy A.1
1 Russian Federal Nuclear Center — All-Russian Research Institute of Experimental Physics, Sarov, 607189, Russia 2Nizhny Novgorod State Technical University n.a. R.E. Alekseev, Nizhny Novgorod, 603950, Russia
3Joint-Stock Company "KBP n.a. Academician A. Shipunov", Tula, 300001, Russia * Corresponding author: Kurkin, Andrey A., e-mail: aakurkin@gmail.com
The application of numerical simulation methods based on the solution of the full three-dimensional Navier —Stokes equations for modelling of wave propagation on the water surface requires the construction of a grid model containing countable nodes throughout the entire volume of water medium. Insufficient grid resolution leads to insufficient detailing of the fields of velocity and pressure, as well as volume fraction of the liquid, which increases the numerical diffusion of the method and, ultimately, leads to an underestimation of oscillation amplitudes of the medium. A large time step also results in a "blurring" of the solution and significantly reduces its stability, especially when using the schemes which compress the front of the media interface.
This paper presents the results of an assessment of acceptable grid sizes and time steps expressed in dimensionless parameters with respect to the wave parameters necessary to ensure accuracy of the solution sufficient for geophysical applications. The estimate is given for the method of calculating three-dimensional multiphase flows with a free surface based on solving the Navier —Stokes equations in a one-velocity approximation based on a completely implicit connection between velocity and pressure using the finite volume method. The finite volume method for the numerical solution of the Navier-Stokes equations is implemented for use on arbitrary unstructured grids. The methodology for estimation of numerical diffusion of the calculation method is proposed. This estimation is expressed as a percentage of the wave amplitude decrease at the distance equal to the one wavelength. In turn the methodology is based on the parameters entered to estimate the acceptable grid sizes and time step for the
© ICT SB RAS, 2019
calculation method. Based on the described methodology, the results of the estimation of the grid resolution in the horizontal and vertical directions, the estimation of the time step, and the evaluation of the influence of the discretization scheme of the convective term are presented.
Keywords: surface waves, volume of fluid, method of finite volume, difference scheme, grid resolution.
Cite: Tyatyushkina, E.S., Kozelkov, A.S., Kurkin, A.A., Kurulin, V.V., Efremov, V.R., Utkin, D.A. Evaluation of numerical diffusion of the finite volume method when modelling surface waves // Computational Technologies. 2019. Vol. 24, No. 1. P. 106-119. (In Russ.) DOI: 10.25743/ICT.2019.24.1.008.
Acknowledgements. The work was supported by a grant from the President of the Russian Federation for state support of scientific research of leading scientific schools of the Russian Federation NSh-2685.2018.5, young Russian doctors of science MD-4874.2018.9 and also with financial support from the RFBR (projects No. 16-0100267 and 17-05-00067).
Received 17 December 2018
