CC BY
23
ISSN 1992-6502 (P ri nt)_
2016. Т. 20, № 3 (73). С. 153-163
Ъьомт, QjrAQnQj
ISSN 2225-2789 (Online) http://journal.ugatu.ac.ru
УДК 532.543
Разработка и примеры приложения специализированного
параллельного кода для численного моделирования турбулентных нестационарных течений со свободной поверхностью
А. И. Храбрый , Д. К. Зайцев , Е. М. Смирнов
aero@phmf.spbstu.ru
ФГАОУ ВО «Санкт-Петербургский политехнический университет Петра Великого» (СПбПУ)
Статья рекомендована к публикации программным комитетом международной научной конференции
«Параллельные вычислительные технологии 2016»
Поступила в редакцию 06.07.2016
Аннотация. Излагается вычислительная методика для моделирования течений жидкости со свободной поверхностью, с акцентом на оригинальные элементы, повышающие точность решения и обеспечивающие корректную работу метода VOF (Volume-Of-Fluid) в особо сложных с вычислительной точки зрения случаях. Методика реализована в специализированном трехмерном коде, оперирующем с неструктурированными расчетными сетками. Параллелизация вычислений осуществлена по технологии «domain decomposition». Даются примеры приложения разработанного кода к решению двумерных и трехмерных задач о нестационарном натекании турбулентного потока воды на различные препятствия. Результаты расчетов сопоставляются с экспериментальными данными.
Ключевые слова: численное моделирование; параллельный код; domain decomposition; течение со свободной поверхностью; натекание на препятствие; турбулентные течения; метод VOF; вязкий отрыв.
ВВЕДЕНИЕ
Многие из интересных с практической точки зрения течений со свободной поверхностью являются нестационарными, сопровождаются натеканием потока на те или иные препятствия, что может приводить к сильной деформации свободной поверхности. Задачи, стоящие перед исследователями, включают как изучение детальной структуры потока, так и определение величин, действующих на препятствия нагрузок.
К настоящему времени доступны экспериментальные данные для ряда модельных задач о натекании потока на препятствия различной формы (треугольник [1], трапеция [2], вертикальная стенка [3] и др.). Как правило, в экспериментальных работах приводятся фотографии свободной поверхности в различные моменты времени.
В последнее время численное моделирование становится применимым ко все более широкому классу течений, заменяя собой эксперимен-
тальные методы исследования. Те или иные численные модели характеризуются различными областями применимости и затратами вычислительных ресурсов. Эксперимент при этом во все большей степени выступает в роли источника данных для валидации математических моделей и программных кодов.
До недавнего времени численное моделирование нестационарного набегания потока воды на препятствия, как правило, выполнялось в рамках приближения «мелкой воды» (см., например, [1, 4]), позволяющего удовлетворительно предсказывать эволюцию подъема жидкости перед препятствием. Однако такие расчеты не способны воспроизвести многие детали взаимодействия потока с препятствием, например, опрокидывание отраженной от препятствия волны, с невысокой точностью предсказывают распределение давления на стенках препятствия (особенно в случае «крутых» препятствий). Кроме того, как отмечено в некоторых из вышеупомянутых работ, при использовании приближения «мелкой
Работа поддержана грантом РФФИ 14-07-00065.
воды» недостаточно точно предсказывается скорость фронта основной волны и отраженных волн.
Альтернативой использованию приближения «мелкой воды» является решение полных (трехмерных или двумерных) уравнений Навье-Стокса, осредненных по Рейнольдсу и замкнутых с привлечением подходящей модели турбулентности. При решении уравнений Навье-Стокса наиболее популярным методом для отслеживания положения свободной поверхности в настоящее время является метод УоЫте-о^ Fluid (VOF) [5], позволяющий проводить расчеты течений с сильной деформацией свободной поверхности (включая опрокидывание волны), обеспечивая при этом высокую детализацию структуры эволюционирующего потока. Сопоставление методик, основанных на приближении «мелкой воды» и на полных уравнениях гидродинамики, проводится, например, в работе [2], где рассматривается натекание потока на трапециевидное препятствие. Результаты данной работы, полученные на основе полных уравнений, свидетельствуют, в частности, о том, что с наветренной стороны препятствия может формироваться отрывная зона, вызванная ростом давления вдоль основного потока при его натекании на препятствие. Модель «мелкой воды» в принципе не позволяет воспроизвести этот эффект.
Применение численного моделирования для детального анализа столь сложного течения, как набегание потока на препятствие, предполагает проверку влияния схемных факторов на получаемое решение. Причем это относится как к задаче расчета формы свободной поверхности, так и к воспроизведению эффектов пристеночного трения. Для рассматриваемых течений характерны высокие числа Рейнольдса, и обычно для определения величины трения на стенках применяется метод пристенных функций (вкупе с «вы-сокорейнольдсовыми» расчетными сетками), а нередко расчеты проводятся и вовсе с заданием условия проскальзывания на стенках. Известно, однако, что пристенные функции могут плохо работать в условиях отрыва пограничного слоя, вызванного «неблагоприятным» градиентом давления.
Получение решения, свободного от влияния схемных факторов, требует весьма густых расчетных сеток в целом по расчетной области с дополнительным сильным сгущением к твердым границам (так называемых «низкорейнольдсовых» сеток). Применительно к нестационарным трехмерным расчетам, это требование сопряжено с особо высокой затратностью вычислительных ресурсов,
что до недавнего времени делало такого рода расчеты невозможными. Сегодня, благодаря распространению мощных многопроцессорных систем, проведение аккуратных расчетов трехмерных нестационарных турбулентных течений со свободной поверхностью становится все более реализуемым. Соответственно, возникает потребность в разработке и апробации эффективных вычислительных средств, позволяющих достичь высокого качества пространственного и временного разрешения течения.
Настоящая работа направлена на разработку основанного на методе VOF численного алгоритма для расчета нестационарных трехмерных течений жидкости со свободной поверхностью, его программную реализацию в параллельном коде и применение к решению модельных задач практической направленности о нестационарном натекании потока воды на различные препятствия. Излагается вычислительная методика с акцентом на оригинальные элементы, повышающие точность решения. Дается описание вычислительных приемов, обеспечивающих корректную работу метода в особо сложных с вычислительной точки зрения случаях.
ЧИСЛЕННАЯ МЕТОДИКА И ПРОГРАММНАЯ РЕАЛИЗАЦИЯ
Реализация и тестирование численных схем для решения системы уравнений метода VOF проводились в программной среде гидро-динамического кода внутреннего пользования Flag-S, разрабатываемого сотрудникам кафедры гидроаэродинамики СПбПУ и исходно не предназначавшегося для расчета течений со свободной поверхностью. В итоге была создана версия кода для расчета течений данного класса, названная Flag-FS (Flag - Free Surface). Новая версия унаследовала от кода Flag-S следующие основные свойства и элементы. Код оперирует неструктурированными расчетными сетками, состоящими из полиэдральных ячеек, что позволяет проводить расчеты течений в областях сложной геометрии. Дискретизация решаемых уравнений выполняется по методу конечных объемов (МКО), связывающему изменение той или иной физической величины в центре ячейки расчетной сетки с потоками этой величины через грани ячейки. Вычисление конвективных и диффузионных потоков основывается на схемах второго порядка, оперирующих значениями самой величины и ее градиента в центрах смежных по грани ячеек. «Перевязывание» скорости и давления в итерационном процессе, реализуемом на каждом шаге по времени, осуществляется по методу
SIMPLEC [6]. Для предотвращения пространственных осцилляций в поле течения применяется коррекция Рхи-Чоу [7]. Реализованы несколько моделей турбулентности, включая используемую в представленных ниже расчетах SST модель Ментера [8] с обобщенными пристенными функциями [9].
Определяющие уравнения и их дискретные аналоги. В методе VOF для определения положения свободной поверхности вводится специальная маркер-функция С, представляющая собой объемную долю жидкости (газу соответствует значение С=0, жидкости - С=1). При расчетах по данному методу движущаяся по сетке межфазная граница становится «размытой», приобретая вид переходного слоя толщиной в несколько расчетных ячеек, в котором значений маркер-функции меняются от 0 до 1. Положение свободной поверхности обычно определяется изоповерхностью С=0,5.
В так называемой одножидкостной формулировке метода VOF уравнения гидродинамики решаются «насквозь», как для единой среды с переменными материальными свойствами, выражаемыми через величину С по следующим соотношениям:
Р = С?1 + (! - С К ( = С( + (1- С ,
(1) (2)
где р - плотность среды, ( - динамическая вязкость среды, индексы I и g относятся к жидкости и газу соответственно.
Поскольку объемная доля жидкости сохраняется вдоль траекторий, ее динамика может быть описана уравнением конвективного переноса (V - скорость потока):
йС йг
С дг
+ v - ус = 0.
Однако, это уравнение не является консервативным и не подходит для использования в рамках концепции МКО. Для реализации в коде Flag-FS уравнение было переписано в консервативной форме (3), одновременно с записью уравнения неразрывности в форме (4), предполагающей несжимаемость жидкости и газа:
— + У - С ) = 0, дг у ' -
У-V = 0.
(3)
(4)
дС
дг
V + £// = 0,
/
Е Г = 0 ,
/
(5)
(6)
где индекс Р относится к центру ячейки, V -ее объем, суммирование проводится по всем граням ячейки, индекс f относится к центру грани, Г = 8/Пу - - объемный расход через грань, Б/ -
площадь грани, и - внешняя нормаль к ней. Отметим, что форма (5) гарантирует сохранение общего количества маркер-функции в расчетной области в целом, таким образом, обеспечивая сохранение объемов, занимаемых жидкостью и газом.
По причинам, детально изложенным в работе [10], уравнение баланса импульса было тождественно преобразовано из традиционной консервативной формы к форме (7) с оригинальным представлением конвективных членов:
р — + У-(руу)-V У-(ру ) = -Ур + У- т + р£, (7)
дг
где р - давление,
ускорение силы тяжести,
Дискретные конечно-объемные аналоги данных уравнений имеют следующий вид:
т = (( + ()(Уг' + (Уу)т) - тензор вязких напряжений, ( - турбулентная вязкость (дается моделью турбулентности).
Форма (7) может быть обобщена на случай уравнений переноса (8) любой физической величины ф - в ее роли может выступать компонента скорости или параметр турбулентности:
р^ + У-^)-фУ-^) = Ш8, (8) дг
где КИБ - это сумма источниковых и диффузионных членов.
Дискретный аналог уравнения (8) имеет следующий вид:
рр%V + ЕР;Ф/Г/ -ФРЪ/Г/ = ш8 (9)
дг / /
Данный способ дискретизации конвективных слагаемых позволяет избежать нефизичных осцилляций в решении, вызываемых большим перепадом значений плотности вблизи межфазной границы.
Фигурирующие в дискретных аналогах определяющих уравнений значения величин на гранях ячеек находятся с помощью процедуры реконструкции, которая задействует значения из центров ближайших ячеек. Для компонент скорости и параметров турбулентности в коде
Flag-FS используются противопоточные схемы второго порядка. Однако такие схемы не подходят для вычисления значений маркер-функции Cf на гранях ячеек, так как вносимые ими эффекты численной диффузии приводят к размытию межфазной границы на множество ячеек расчетной сетки. В литературе предложен ряд специализированных, т.н. «сжимающих» численных схем для нахождения значений Cf. Большинство таких схем удовлетворяют т.н. CBC критерию [11], который гарантирует локальную ограниченность получаемого решения в случае одномерной постановки и равномерной расчетной сетки. Однако свойства схем в многомерном случае и при расчетах на неравномерных/скошенных сетках каждый раз требуют отдельного изучения. В частности, результаты детального сопоставления возможностей двух популярных схем HRIC [12] и CICSAM [13], а также предложенной позднее схемы M-CICSAM [14], представлены в работе [15]. В проведенных тестах схема M-CICSAM продемонстрировала явное превосходство над другими двумя схемами: меньшую зависимость качества результатов расчетов от густоты и скошенности расчетной сетки, величины шага по времени. С учетом результатов работы [15], для реализации в коде Flag-FS была выбрана схема M-CICSAM в сочетании со схемой Кранка-Николсон для аппроксимации по времени. Такая комбинация хорошо работает в широком диапазоне чисел Куранта, практически не искажая межфазную границу и не допуская сколько-нибудь заметных «вылетов» значений маркер-функции за пределы диапазона [0, 1]. Во избежание получения нефизичных (отрицательных) значений плотности и вязкости, в выражениях (1) и (2) упомянутые «вылеты» (которые, как правило, не превосходили 0,001) обрезались.
Кроме того, для устранения небольшого остаточного размытия межфазной границы (которое, как показано в работе [16], в отдельных случаях может приводить к существенному искажению поля течения), в коде Flag-FS была реализована оригинальная методика «обострения» фронта [16]. Отметим, что схема Кранка-Ни-колсон использовалась для аппроксимации и остальных уравнений, что позволило обеспечить достаточно высокую точность решения при числах Куранта, достигающих 0,8. Именно с таким значением числа Куранта проводились представленные ниже расчеты (шаг по времени подирался автоматически). Вопросы чувствительности решения к густоте расчетной сетки исследовались в работе [17] на тестовой задаче натека-ния потока на квадратное препятствие.
Методические расчеты показали также, что при использовании «низкорейнольдсовых» сеток для решения задач о растекании жидкости вдоль сухой твердой стенки, вблизи последней формируется тонкая нефизичная газовая прослойка, приводящая к существенному искажению величины трения на стенке [18]. Для устранения описанного артефакта в коде Flag-FS реализована оригинальная методика, вводящая искусственную диффузию маркер-функции, которая действует только вблизи твердой стенки, по направлению к ней:
~dt
■ +
V- (С v )= max
д
кдпь
X
дС
\ Л
дп
(10)
b J J
где йь = Уй - орт направления, перпендикулярного ближайшей твердой границе, % - коэффициент искусственной диффузии, определяемый по формуле:
X =
0
d > d
cell
X
d „— d
-ceä- , d < d ,,
max i ' cell
d
cell
где й - расстояние до ближайшей стенки, йоеп и %тах - пользовательские параметры. Параметр йоеи определяет расстояние от стенки, на котором действует искусственная диффузия, %тах задает максимальное значение коэффициента искусственной диффузии, достигаемое на стенке.
Дискретный аналог уравнения (10) имеет следующий вид:
dt
V + I CfFf =
(11)
: maxW(xf (VVCf - Vd)(nf - Vd)Sf ), 0 I.
f
Алгоритм вычислений и вопросы распараллеливания. Совместное решение определяющих уравнений (3), (4), (8) в коде Flag-FS осуществляется посредством организации глобального итерационного процесса, на каждой итерации которого решаются линейные уравнения относительной приращений скорости, давления и параметров турбулентности для нахождения уточненных значений данных величин. «Перевязка» приращений скорости и давления осуществляется по методу SIMPLEC.
Поскольку для корректного решения уравнения переноса маркер-функции (3) требуется бездивергентное поле скорости, обновление поля
0
Блок 1 Блок 2
Рис. 1. Передача данн
маркер-функции проводится не на каждой глобальной итерации по решению уравнений гидродинамики (4) и (8), а с некоторым шагом, достаточным для получения поля скорости, близкого к дивергентному (как правило, шаг составлял десяток итераций, а общее количество глобальных итераций доходило до 100). Исходя из того, что ограничения на величину шага по времени для уравнений гидродинамики, как правило, мягче, чем для уравнения (3), была реализована методика, позволяющая, в целях экономии вычислительных ресурсов, выполнять дробные шаги по времени для уравнения (3) в пределах одного шага по времени для решения уравнений гидродинамики [10]. На каждом из этих дробных шагов дискретный аналог уравнения переноса маркер-функции - уравнение (11) - решается по методу установления (как правило, требовалось 35 итераций).
Для решения систем линейных уравнений в коде Flag-FS используются итерационные методы на основе подпространств Крылова. Уравнение Пуассона для приращения давления в методе SIMPLEC решается с помощью метода сопряженных градиентов (CG) с неполным разложением Холецкого в качестве предобуслав-ливателя (см. например, [19]). Для решения остальных уравнений используется вариант метода би-сопряженных градиентов (BiCGStab, [20]) с предобуславливателем SGS (Symmetric Gauss-Seidel).
Распараллеливание вычислений производится по технологии «domain decomposition», в рамках которой расчетная область разбивается на блоки примерно равного размера, и для выполнения вычислений внутри каждого блока создается отдельный процесс, оперирующий своим набором данных в оперативной памяти и выполняющийся на отдельном процессорном ядре. Для обеспечения связности численного решения производятся обмены данными между блоками (стыковка) с удовлетворением требованию «про-
Блок 1 «Заграничные» Блок 2 ячейки
при стыковке блоков
зрачности» границ, т.е. отсутствия влияния разбиения расчетной области на получаемое решение.
В коде Flag-S разбиение расчетной стеки проводится на соприкасающиеся (не накладывающиеся) блоки. Связанность боков обеспечивается созданием в каждом из них ряда «заграничных» ячеек, повторяющих приграничные ячейки смежного блока (рис. 1). Для прозрачности границ требуется, чтобы потоки величин через грани ячеек на межблочных границах были одинаковыми для смежных ячеек, находящихся в разных блоках. Поскольку в рамках реализованных в коде Flag-S численных схем второго порядка для вычисления потока величины через смежную грань двух ячеек требуются лишь значения самой величины и ее градиента в центрах этих ячеек, прозрачность границ обеспечивается копированием приграничных значений величин и их градиентов. Учитывая, что при вычислении градиента величины в ячейке сетки используются значения этой величины из соседних ячеек, стыковка блоков осуществляется в два этапа: сначала блоки обмениваются значениями величин, по которым в каждом из блоков вычисляются их градиенты, а затем передаются вычисленные градиенты.
Обмены данных между блоками производятся по технологии MPI. Для уменьшения числа обращений к сравнительно медленным функциям MPI, все данные по каждой стыковке передаются единым массивом; при этом для пересылки используются неблокирующие процедуры MPI_ISEND и MPI_IRECV, так что каждая стыковка обрабатывается по мере готовности данных в стыкуемых блоках, независимо от других стыковок.
РАСЧЕТЫ НАТЕКАНИЯ ПОТОКА НА ТРЕУГОЛЬНОЕ ПРЕПЯТСТВИЕ
В качестве одного из приложений разработанного кода выступала тестовая задача о нате-кании потока, вызванного обрушением дамбы, на треугольное препятствие. Постановка задачи
соответствовала эксперименту [1]. Схема экспериментальной установки приведена на рис. 2.
Торцевая стенка Боковые стенки
Рис. 2. Схема экспериментальной установки работы [1]. Все размеры даны в см
Перед началом эксперимента вода удерживалась вертикальной перегородкой («дамбой») в левой части длинного канала прямоугольного сечения. Затем перегородка резко убиралась, вода начинала распространяться вправо, натекая на препятствие и перетекая через него. Боковые стенки канала были прозрачными; положения свободной поверхности в различные моменты времени фиксировались при помощи фотокамеры, расположенной сбоку. Число Рейнольдса, определенное по начальной высоте столба жидкости Н и характерной скорости потока (2gH)l¡2, составляло 1,6 105.
Представляемые ниже расчеты, проведенные в двумерной и трехмерной постановках, были направлены на исследование влияния вязких эффектов вблизи дна и боковых стенок канала на эволюционирующую картину течения. В этих расчетах пользовательский параметр йоги, входящий в модель искусственной диффузии маркер-функции у стенок (10), полагался равным половине характерного размера ячеек вдоль стенки, а коэффициент искусственной диффузии %тах = а^Н3)112, где значение коэффициента а=10-4 подобрано эмпирически.
Результаты расчетов в двумерной постановке. Совместное решение определяющих уравнений (3), (4), (8) в коде Flag-FS осуществляется посредством организации глобального итерационного процесса, на каждой итерации которого решаются линейные уравнения относительной приращений скорости, давления и параметров турбулентности для нахождения уточненных значений данных величин. «Перевязка» приращений скорости и давления осуществляется по методу SIMPLEC.
Двумерные расчеты были проведены в двух вариантах: с постановкой условий прилипания на стенках (и, соответственно, с учетом пристеночных вязких эффектов) и с использованием условий проскальзывания. В первом случае сетка имела сгущение к ниж-
ней границе, обеспечивающее значения нормализованного расстояния от стенки до центров пристенной ячейки у+= у1у}т^р/ц, не
превосходящие единицы, и состояла из 25 000 прямоугольных ячеек (с небольшой скошенностью в области над препятствием). Во втором случае использовалась сетка из 16 000 ячеек, без сгущения.
Форма свободной поверхности, полученная в расчетах с разными граничными условиями, иллюстрируется на рис. 3. Результаты расчетов наложены на экспериментальные фотографии из работы [1] для двух моментов времени, отсчитываемых от устранения удерживавшей воду перегородки. Как видно из рисунка, форма свободной поверхности в расчете с условиями проскальзывания (рис. 3, а, в) весьма сильно отличается от наблюдаемой в эксперименте. Учет эффектов придонного трения приводит к существенному улучшению качества предсказания течения. Основное отличие между двумя решениями состоит в появлении отрывных зон в расчете с учетом придонного трения. Рис. 3, б, г показывает, что сначала одна ^ = 3,0 c), а затем две ^ = 3,7 о) крупные отрывные зоны возникают перед препятствием и на его наветренной стороне. Эти зоны приводят к возникновению «холмов» на свободной поверхности, видимых и на фотографиях.
Рис. 3. Форма свободной поверхности и отрывных зон для моментов времени 3,0 с (а, б) и 3,7 с (в, г). Результаты расчетов с условием проскальзывания (а, в) и с условием прилипания (б, г) наложены на фотографии эксперимента [1]
Более детально картина течения, полученная в расчете с учетом эффектов придонного трения, показана на рис. 4, где приведено векторное поле количества движения (отображены векторы количества движения в каждой третьей расчетной ячейке по обоим направле-
ниям). Из рисунка видно, что придонные градиенты скорости хорошо разрешены на сетке. Можно отметить также, что зоны возвратного (рециркуляционного) движения весьма велики, причем в момент времени 3,7 с внутри левой отрывной зоны виден дополнительный небольшой вихрь вблизи дна.
Рис. 4. Форма свободной поверхности и векторы количества движения при / = 3,0 с (а) и / = 3,7 с (б),
Трехмерные расчеты проводились с целью анализа влияния вязких эффектов у боковых стенок канала на структуру течения и форму свободной поверхности. Рассматриваемая конфигурация обладает симметрией относительно средней вертикальной плоскости. Это позволило использовать расчетную область, соответствующую половине экспериментальной. Расчетная сетка состояла из одного миллиона шестигранных ячеек и была сгущена к нижней и к боковой стенкам; степень сгущения к обеим границам была такой же, как в двумерном расчете.
В расчете были задействованы 62 вычислительных ядра кластера с процессорами типа AMD Opteron 6300 (ускорение примерно в 45 раз по сравнению с однопроцессорным вариантом); общее время счета - около недели (2000 шагов по времени).
На рис. 5 представлена расчетная форма свободной поверхности в те же моменты времени, что и на рис. 3, а на рис. 6 дана векторная картина скоростного поля в момент времени 3,0 с. Видно, что, как и в двумерном расчете, возникают сначала одна, а потом две отрывные зоны, приводящие к появлению «холмов» на свободной поверхности. Примерно на половине полуширины канала, со стороны плоскости симметрии, высота подъема свободной поверхности на «холмах» практически постоянна поперек канала. Однако по мере приближения к боковой стенке свободная поверхность существенно отклоняется от двумерной формы, а «холмы» практически вырождаются.
Рис. 5. Свободная поверхность по результатам трехмерных расчетов для моментов времени
(а) / = 3,0 с и (б) / = 3,7 с. Также показаны вихревые линии (пунктир), идущие из центров отрывных зон в плоскости симметрии
Рис. 6. Расчетное поле скорости на свободной поверхности и в плоскости симметрии, / = 3,0 с
Трехмерные эффекты еще в большей степени проявляются в конфигурации отрывных зон. Исходящие из центров отрывных зон в плоскости симметрии вихревые линии (обозначены пунктиром на рис. 5) на большей части своей длины отклоняются от прямой линии, которая имела бы место в двумерном течении, а вблизи боковой стенки канала «поворачивают» наверх и выходят на поверхность жидкости. Свидетельство выхода
«оси» отрывной зоны на свободную поверхность можно видеть и в картине течения на свободной поверхности для момента времени 3,0 с (рис. 6), где отчетливо видна область рециркуляционного движения вблизи боковой стенки (над наветренной стороной препятствия).
На рис. 7 для того же момента времени приведено поле поверхностного трения на дне канала, а именно: на наветренной стороне препятствия (ограниченной линиями А-А и Б-Б) и на примыкающем горизонтальном участке канала.
Т7: -100 -80 -60 -40 -20 0 20 40 60 Ра
Рис. 7. Расчетное распределение поверхностного
трения на наветренной стороне препятствия и на близлежащей части дна канала при t = 3,0 с: (а) поле продольной компоненты трения, (б) векторы трения вблизи передней кромки препятствия
Как видно из рис. 7, а, вблизи плоскости симметрии (где течение близко к двумерному), продольная компонента трения, т™, имеет два локальных минимума. Первый из них располагается вблизи кромки препятствия, а второй, с большим по модулю отрицательным значением, соответствует ядру рециркуляционной зоны, показанной на рис. 4, а.
Векторное поле трения на стенке (рис. 9, б) свидетельствует о наличии трехмерной отрывной зоны вблизи точки контакта передней кромки клина и боковой поверхности. Отметим также, что в этой области трение на стенке достигает максимального по модулю отрицательного значения.
На рис. 8 расчетная форма свободной поверхности (вид сбоку) сопоставляется с экспериментальными фотографиями [1]. Здесь свободная поверхность представляется не линией (как в двумерном случае), а имеет вид ленты, изменения в видимой толщине которой отражают зависимость высоты подъема жидкости от трансвер-сальной координаты. Примечательно, что сходные «ленты» видны и на экспериментальных фотографиях (темные области вдоль свободной поверхности). Сопоставление рис. 3 и 8 позволяет заключить, что результаты трехмерного расчета намного лучше согласуются с экспериментальными данными.
Рис. 8. Результаты трехмерных расчетов течения, соответствующего эксперименту [1] (форма свободной поверхности) в сопоставлении
с экспериментальными фотографиями для моментов времени 3,0 с (а) и 3,7 с (б)
РАСЧЕТЫ НАТЕКАНИЯ ПОТОКА НА МНОЖЕСТВЕННЫЕ (ОДНОРЯДНЫЕ И ДВУРЯДНЫЕ) ПРЕПЯТСТВИЯ
Другой пример приложения разработанного кода - это модельная задача о натекании потока воды, возникшего в результате разрушения дамбы, на препятствия в форме параллелепипедов, выстроенных поперек потока в один и два бесконечных ряда (во втором случае в шахматном порядке, рис. 9). Данная модель имитирует постройки в прибрежной зоне. Основания препятствий имели размеры 16x16 см, расстояние между рядами препятствий равнялось 16 см, шаг между препятствиями одного ряда -32 см. Число Рейнольдса Яе = рI ^Н3)1/2 /ц составляло 1,8 106. В расчетах с одним рядом отсутствовал перед-ний ряд (кубических) препятствий.
Рис. 9. Постановка задачи о натекании воды на ряд препятствий. Все размеры даны в см
С учетом присущей задаче симметрии, расчетная область в поперечном направлении располагалась от середины низкого препятствия до середины соседнего с ним высокого препятствия (на продольных вертикальных границах области ставились условия симметрии). Для обеих конфигураций расчетные сетки состояли из примерно полумиллиона шестигранных ячеек и
Рис. 10. Положения свободной поверхности и поля давления на лицевой стороне высокого препятствия в моменты времени 0,8 с (верхний ряд), 1,2 с (средний ряд) и 2,0 с (нижний ряд): результаты расчетов в отсутствие первого ряда препятствий (слева) и при его наличии (справа)
162
информатика, вычислительная техника и управление
имели одинаковое сгущение ко дну, обеспечивающее значения у+, не превосходящие единицы. Результаты расчетов для вариантов с одним и двумя рядами препятствий приведены на рис. 10 (для удобства восприятия визуализированных картин течения расчетные поля размножены с учетом симметрии, присущей задаче).
Как видно из рисунка, характер течений, рассчитанных для вариантов с одним и двумя рядами препятствий, существенно различен. В отсутствие первого ряда поток ударяется в нижнюю часть высоких препятствий, часть жидкости движется вверх вдоль стенок, а затем опрокидывается назад, навстречу потоку. При наличии первого ряда из кубических препятствий поток сначала ударяется о них. Часть потока «перелетает» через низкие препятсвия и ударяется о второй ряд препятствий, достигая середины их высоты. Различия в картине течений заметно проявляются в поле давления на наветренной стороне высоких препятствий. В отсутствие первого ряда область максимальных значения давления располагается вблизи дна. При наличии первого ряда зона высокого давления несколько уменьшается, однако перемещается выше, создавая, таким образом, больший опрокидывающий момент (12,1 Нм против 4,8 Нм при г = 0,8 о). Очевидно, последнее представляет большую опасность с точки зрения потенциального разрушения прибрежных зданий под действием мощной набегающей волны.
ВЫВОДЫ
На базе кода внутреннего пользования Flag-S создана специализированная версия, позволяющая на основе метода VOF проводить параллельные расчеты трехмерных и двумерных нестационарных течений жидкости со свободной поверхностью. Разработан (и реализован) ряд оригинальных составляющих метода VOF для обеспечения его корректной работы в сложных с вычислительной точки зрения случаях. Произведен поиск оптимальных численных схем в плане чувствительности результатов расчетов к густоте используемой расчетной сетки и шагам по времени. Разработанный программный код применен к решению тестовых и модельных задач о нестационарном натекании потока воды препятствия различной формы, в частности, на одиночное треугольное препятствие и на множественные (однорядные и двурядные) препятствия в форме параллелепипедов.
На примере задачи с треугольным препятствием показано, что пристенные вязкие эффекты могут оказывать весьма сильное влияние
на картину течения через формирование крупных отрывных зон, приводящих к возникновению «холмов» на свободной поверхности. Установлено, что вязкие эффекты вблизи боковых стенок канала делают картину течения существенно трехмерной. Результаты трехмерного расчета хорошо согласуются с экспериментальными данными.
Для натекания потока на множественные препятствия получен интересный для практики вывод о том, что наличие первого ряда препятствий может привести к увеличению опрокидывающего момента, действующего на препятствия второго ряда.
СПИСОК ЛИТЕРАТУРЫ
1. Soares-Frazao S. Experiments of dam-break wave over a triangular bottom sill // Journal of Hydraulic Research. 2007. Vol. 45, extra issue. P. 19-26.
2. Ozmen-Cagatay H., Kocaman S. Dam-break flow in the presence of obstacle: experiment and CFD simulation // Engineering Applications of Computational Fluid Mechanics. 2011. Vol. 5. P. 541-552.
3. Hu C., Sueyoshi M. Numerical simulation and experiment on dam break problem // Journal of Marine Science and Application. 2010. Vol. 9. P. 109-114.
4. Федотова З. И., Чубаров Л. Б. Численное моделирование наката цунами // Труды Международной конференции RDAMM-2001, 2001. № 6, ч. 2, спец. выпуск. С. 380-396. [ Z. I. Fedotova, L. B. Chubarov, "Numerical modeling of tsunami rolling" (in Russian), in Proc. Int. Conf. RDAMM-2001, 2001, vol. 6, spec. issue, part 2, pp. 380-396. ]
5. Hirt C. W., Nichols B. D. Volume of fluid (VOF) method for the dynamics of free boundaries // Journal of Computational Physics. 1981. Vol. 39. P. 201-226.
6. Vandoormaal J. P., Raithby G. D. Enhancements of the SIMPLE method for predicting incompressible fluid flows // Numerical Heat Transfer. 1984. Vol. 7, issue 2. P. 147-163.
7. Rhie C. M., Chow W. L. A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation // AIAA Journal. 1983. Vol. 21. P. 1525-1532.
8. Menter F. R. Two equation eddy-viscosity turbulence models for engineering applications // AIAA Journal. 1994. Vol. 32. P. 1598-1605.
9. Smirnov E. M., Zaitsev D. K. Modification of wall boundary conditions for low-Re k-e turbulence models aimed at grid sensitivity reduction // Proc. European Conf. for Aerospace Sciences (EUCASS 2005), 4-7 July 2005, Moscow. CD-ROM, ID 2.09.06. 7 p.
10. Храбрый А. И., Зайцев Д. К., Смирнов Е. М. Численное моделирование течений со свободной поверхностью на основе метода VOF // Труды Крыловского государственного научного центра. 2013. Вып. 78. С. 53-64. [ A. I. Khrabry, E. M. Smirnov, D. K. Zaytsev "Numerical simulation of free surface flows on the base of the VOF method" (in Russian), in Proc. of Krylov state research centre, vol. 78, pp. 53-64, 2013. ]
11. Gaskell P. H., Lau A. K. C. Curvature-compensated con-vective-transport: SMART, A new boundedness-preserving transport algorithm // Int. J. Numer. Methods Fluids. 1988. Vol. 8. Issue 6. P. 617-641.
12. A two-fluid Navier-Stokes solver to simulate water entry / Muzaferija S. [et al.] // Proc. 22th Symp. on Naval Hydrodynamics, Washington DC, August 9-14 1998. P. 638-651.
13. Ubbink O., Issa I. A method for capturing sharp fluid interfaces on arbitrary meshes // Journal of Computational Physics. 1999. Vol. 153. P. 26-50.
14. Waclawczyk T., Koronowicz T. Remarks on prediction of wave drag using VOF method with interface capturing approach // Archives of Civil and Mechanical Engineering. 2008. Vol. 8. P. 5-14.
15. Khrabry A. I., Smirnov E. M., Zaytsev D. K. Solving the Convective Transport Equation with Several High-Resolution Finite Volume Schemes: Test Computations / In: Computational Fluid Dynamics 2010. Ed.: A. Kuzmin. New-York: Springer, 2011. 954 p. P. 535-540.
16. Free surface flow computations using the M-CICSAM scheme added with a sharpening procedure / Khrabry A. I. [et al.] // Proc. 6th European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2012), Vienna, Austria, 10-14 September 2012. CD-ROM. ISBN: 978-39502481-9-7. 2 p.
17. Расчет течений со свободными поверхностями: влияние схемных факторов и модели турбулентности / Зайцев Д. К. [и др.] // Тр. 14-й междунар. конф. «Супервычисления и математическое моделирование», Саров, 1-5 октября 2012. С. 282-292. [ Zaytsev D. K. et al., "Numerical simulation of free surface flows: influence of scheme factors and turbulence model" (in Russian), in Proc. 14th Int. Conf. "Supercomputing and mathematical modeling", Sarov, Russia, 2012, pp. 282-292. ]
18. Numerical study of separation phenomena in the dam-break flow interacting with a triangular obstacle / Khrabry A. [et al.] // Proc. 16th Int. Conf. on Fluid Flow Technologies (CMFF'15), Budapest, Hungary, September 1-4 2015. CD-ROM, CMFF15-144.pdf. 8 p.
19. Голуб Дж., Ван Лоун Ч. Матричные Вычисления. М.: Мир, 1999. 553 с. [G. Golub, С. Van Loan, Matrix Computations, (in Russian). Moscow: Mir publishing, 1999, 553 p. ]
20. Van Der Vorst H. A. Bi-CGSTAB: A fast and smoothly converging variant of Bi CG for the solution of nonsymmetric linear systems // SIAM Journal on Scientific and Statistical Computing. 1992. Vol. 13. P. 631-644.
ОБ АВТОРАХ
ХРАБРЫЙ Александр Иосифович, вед. программист каф. Гидроаэродинамики, канд. физ.-мат. наук по механике жидкости и газа (СПбПУ, 2015).
ЗАЙЦЕВ Дмитрий Кириллович, доцент каф. Гидроаэродинамики, канд. физ.-мат. наук по механике жидкости и газа (СПбПУ, 1986).
СМИРНОВ Евгений Михайлович, зав. каф. Гидроаэродинамики, проф., д-р физ.-мат. наук по механике жидкости и газа (СПбПУ, 1988).
METADATA
Title: Development and application of a specialized parallel code for numerical simulation of unsteady turbulent free-surface flows.
Authors: A. I. Khrabry1, D. K. Zaytsev1, E. M. Smirnov1
Affiliation:
Saint-Petersburg Polytechnic University (SPbPU), Russia.
Email: 1 aero@phmf.spbstu.ru.
Language: Russian.
Source: Vestnik UGATU (scientific journal of Ufa State Aviation Technical University), vol. 20, no. 3 (73), pp. 153-163, 2016. ISSN 2225-2789 (Online), ISSN 1992-6502 (Print).
Abstract: Computational method for free-surface flow modeling is presented with an emphasis on original developments aimed at improvement of solution accuracy and robustness of the VOF (Volume-Of-Fluid) method in cases of extra computational complexity. The computational techniques developed have been implemented in a specialized 3D unstructured-grid parallel code. Parallelization of computations is based on the "domain decomposition" technique. Examples of the code application to modeling of 2D and 3D unsteady turbulent flows interacting with obstacles of different shapes are presented. Computational results are compared with experimental data.
Key words: numerical simulation; parallel code; domain decomposition; free-surface flow; flow-obstacle interaction; turbulent flows; VOF method; viscous separation.
About authors:
KHRABRY, Alexander Iosifovich, Programmer, Dept. of Hydro-and Aerodynamics. PhD in Phys. & Math. (SPbPU, 2015).
ZAITSEV, Dmitri Kirillovich, Assistant prof., Dept. of Hydro- and Aerodynamics. PhD in Phys. & Math. (SPbPU, 1978).
SMIRNOV, Evgueni Mikhailovich, The head of Dept. of Hydro-and Aerodynamics, prof., Dr. of Science in Phys. & Math. (SPbPU, 1988).
