Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Численное моделирование трехмерных нестационарных течений со свободной поверхностью: итоги разработки и примеры приложения специализированного параллельного кода А.И. Храбрый, Д.К. Зайцев, Е.М. Смирнов Санкт-Петербургский политехнический университет Петра Великого Излагается вычислительная методика для моделирования течений жидкости со сво- бодной поверхностью, с акцентом на оригинальные элементы, повышающие точ- ность решения и обеспечивающие корректную работу метода VOF (Volume-Of-Fluid) в особо сложных с вычислительной точки зрения случаях. Методика реализована в специализированном трехмерном коде, оперирующем с неструктурированными рас- четными сетками. Параллелизация вычислений осуществлена по технологии “domain decomposition”. Даются примеры приложения разработанного кода к решению дву- мерных и трехмерных задач о нестационарном натекании турбулентного потока воды на различные препятствия. Результаты расчетов сопоставляются с эксперименталь- ными данными. Ключевые слова: численное моделирование, параллельный код, domain decomposition, течение со свободной поверхностью, натекание на препятствие, тур- булентные течения, метод VOF, вязкий отрыв. 1. Введение Многие из интересных с практической точки зрения течений со свободной поверхностью являются нестационарными, сопровождаются натеканием потока на те или иные препятствия, что может приводить к сильной деформации свободной поверхности. Задачи, стоящие перед исследователями, включают как изучение детальной структуры потока, так и определение ве- личин действующих на препятствия нагрузок. К настоящему времени доступны экспериментальные данные для ряда модельных задач о натекании потока на препятствия различной формы (треугольник – [1], трапеция [2], вертикаль- ная стенка [3] и др.). Как правило, в экспериментальных работах приводятся фотографии сво- бодной поверхности в различные моменты времени. В последнее время численное моделирование становится применимым ко все более широ- кому классу течений, заменяя собой экспериментальные методы исследования. Те или иные численные модели характеризуются различными областями применимости и затратами вычис- лительных ресурсов. Эксперимент при этом во все большей степени выступает в роли источни- ка данных для валидации математических моделей и программных кодов. До недавнего времени численное моделирование нестационарного набегания потока воды на препятствия, как правило, выполнялось в рамках приближения «мелкой воды» (см., к при- меру, [1, 4]), позволяющего удовлетворительно предсказывать эволюцию подъема жидкости перед препятствием. Однако такие расчеты не способны воспроизвести многие детали взаимо- действия потока с препятствием, например, опрокидывание отраженной от препятствия волны, с невысокой точностью предсказывают распределение давления на стенках препятствия (особенно в случае «крутых» препятствий). Кроме того, как отмечено в некоторых из вышеупомянутых ра- бот, при использовании приближения «мелкой воды» недостаточно точно предсказывается ско- рость фронта основной волны и отраженных волн. Альтернативой использованию приближения «мелкой воды» является решение полных (трехмерных или двумерных) уравнений Навье-Стокса, осредненных по Рейнольдсу и замкнутых с привлечением подходящей модели турбулентности. При решении уравнений Навье-Стокса наиболее популярным методом для отслеживания положения свободной поверхности в настоя- 361 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt щее время является метод Volume-of-Fluid (VOF) [5], позволяющий проводить расчеты течений с сильной деформацией свободной поверхности (включая опрокидывание волны), обеспечивая при этом высокую детализацию структуры эволюционирующего потока. Сопоставление методик, основанных на приближении «мелкой воды» и на полных уравнениях гидродинамики, проводит- ся, например, в работе [2], где рассматривается натекания потока на трапециевидное препятствие. Результаты данной работы, полученные на основе полных уравнений, свидетельствуют, в частно- сти, о том, что с наветренной стороны препятствия может формироваться отрывная зона, вызван- ная ростом давления вдоль основного потока при его натекании на препятствие. Модель «мелкой воды» в принципе не позволяет воспроизвести этот эффект. Применение численного моделирования для детального анализа столь сложного течения как набегание потока на препятствие предполагает проверку влияния схемных факторов на получаемое решение. Причем это относится как к задаче расчета формы свободной поверхности, так и к воспроизведению эффектов пристеночного трения. Для рассматриваемых течений характерны высокие числа Рейнольдса, и обычно для определения величины трения на стенках применяется метод пристенных функций (вкупе с «высокорейнольдсовыми» расчетными сетками), а нередко расчеты проводятся и вовсе с заданием условия проскальзывания на стенках. Известно, однако, что пристенные функции могут плохо работать в условиях отрыва пограничного слоя, вызванного «неблагоприятным» градиентом давления. Получение решения, свободного от влияния схемных факторов, требует весьма густых рас- четных сеток в целом по расчетной области с дополнительным сильным сгущением к твердым границам (так называемых «низкорейнольдсовых» сеток). Применительно к нестационарным трехмерным расчетам, это требование сопряжено с особо высокой затратностью вычислитель- ных ресурсов, что до недавнего времени делало такого рода расчеты невозможными. Сегодня, благодаря распространению мощных многопроцессорных систем, проведение аккуратных рас- четов трехмерных нестационарных турбулентных течений со свободной поверхностью стано- вится все более реализуемым. Соответственно, возникает потребность в разработке и апроба- ции эффективных вычислительных средств, позволяющих достичь высокого качества про- странственного и временного разрешения течения. Настоящая работа направлена на разработку основанного на методе VOF численного алго- ритма для расчета нестационарных трехмерных течений жидкости со свободной поверхностью, его программную реализацию в параллельном коде и применение к решению модельных задач практической направленности о нестационарном натекании потока воды на различные препят- ствия. Излагается вычислительная методика с акцентом на оригинальные элементы, повышаю- щие точность решения. Дается описание вычислительных приемов, обеспечивающих коррект- ную работу метода в особо сложных с вычислительной точки зрения случаях. 2. Численная методика и программная реализация Реализация и тестирование численных схем для решения системы уравнений метода VOF проводились в программной среде гидродинамического кода внутреннего пользования Flag-S, разрабатываемого сотрудникам кафедры гидроаэродинамики СПбПУ и исходно не предназна- чавшегося для расчета течений со свободной поверхностью. В итоге была создана версия кода для расчета течений данного класса, названная Flag-FS (Flag – Free Surface). Новая версия уна- следовала от кода Flag-S следующие основные свойства и элементы. Код оперирует неструкту- рированными расчетными сетками, состоящими из полиэдральных ячеек, что позволяет прово- дить расчеты течений в областях сложной геометрии. Дискретизация решаемых уравнений вы- полняется по методу конечных объемов (МКО), связывающему изменение той или иной физи- ческой величины в центре ячейки расчетной сетки с потоками этой величины через грани ячей- ки. Вычисление конвективных и диффузионных потоков основывается на схемах второго по- рядка, оперирующих значениями самой величины и ее градиента в центрах смежных по грани ячеек. «Перевязывание» скорости и давления в итерационном процессе, реализуемом на каж- дом шаге по времени, осуществляется по методу SIMPLEC [6]. Для предотвращения простран- ственных осцилляций в поле течения применяется коррекция Рхи-Чоу [7]. Реализованы не- сколько моделей турбулентности, включая используемую в представленных ниже расчетах SST модель Ментера [8] с обобщенными пристенными функциями [9]. 362 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 2.1 Определяющие уравнения и их дискретные аналоги В методе VOF для определения положения свободной поверхности вводится специальная маркер-функция C, представляющая собой объемную долю жидкости (газу соответствует зна- чение C=0, жидкости – C=1). При расчетах по данному методу движущаяся по сетке межфазная граница становится «размытой», приобретая вид переходного слоя толщиной в несколько рас- четных ячеек, в котором значений маркер-функции меняются от 0 до 1. Положение свободной поверхности обычно определяется изоповерхностью С=0.5. В так называемой одножидкостной формулировке метода VOF уравнения гидродинамики решаются «насквозь», как для единой среды с переменными материальными свойствами, вы- ражаемыми через величину C по следующим соотношениям:   Cl  1  С  g (1)   Cl  1  С  g (2) Здесь  – плотность среды,  – динамическая вязкость среды, индексы «l» и «g» относятся к жидкости и газу соответственно. Поскольку объемная доля жидкости сохраняется вдоль траекторий, ее динамика может  быть описана уравнением конвективного переноса ( v – скорость потока): dC C    v  C  0 dt t Однако, это уравнение не является консервативным и не подходит для использования в рамках концепции МКО. Для реализации в коде Flag-FS уравнение было переписано в консервативной форме (3), одновременно с записью уравнения неразрывности в форме (4), предполагающей несжимаемость жидкости и газа: C     Cv   0 (3) t  v  0 (4) Дискретные конечно-объемные аналоги данных уравнений имеют следующий вид: С p t V C F  0 f f f (5) F  0 f f (6) Здесь индекс P относится к центру ячейки, V – ее объем, суммирование проводится по всем   граням ячейки, индекс f относится к центру грани, F f  S f n f  v f – объемный расход через  грань, Sf – площадь грани, n f – внешняя нормаль к ней. Отметим, что форма (5) гарантирует сохранение общего количества маркер-функции в расчетной области в целом, таким образом, обеспечивая сохранение объемов, занимаемых жидкостью и газом. По причинам, детально изложенным в работе [10], уравнение баланса импульса было тож- дественно преобразовано из традиционной консервативной формы к форме (7) с оригинальным представлением конвективных членов.  v         v v   v   v   p    τ  g (7) t      Здесь p – давление, g – ускорение силы тяжести, τ    t  v  v  – тензор вязких напряжений, t – турбулентная вязкость (дается моделью турбулентности). Форма (7) может быть обобщена на случай уравнений переноса (8) любой физической ве- личины  – в ее роли может выступать компонента скорости или параметр турбулентности: 363 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt        v      v   RHS , (8) t где RHS – это сумма источниковых и диффузионных членов. Дискретный аналог уравнения (8) имеет следующий вид:  p p t V    F     F  RHS f f f f p f f f (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 реализо- вана оригинальная методика, вводящая искусственную диффузию маркер-функции, которая действует только вблизи твердой стенки, по направлению к ней: С     C      C v   max    , 0  .  (10) t  nb  nb    Здесь nb  d – орт направления, перпендикулярного ближайшей твердой границе,  – ко- эффициент искусственной диффузии, определяемый по формуле: 364 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt  0 , d  d cell  d  d ,  max cell , d  d cell   d cell где d – расстояние до ближайшей стенки, dcell и max – пользовательские параметры. Параметр dcell определяет расстояние от стенки, на котором действует искусственная диффузия, max зада- ет максимальное значение коэффициента искусственной диффузии, достигаемое на стенке. Дискретный аналог уравнения (10) имеет следующий вид: С p    C F  max   C  d n  d S , 0   V (11) t f f f f f f f  f  2.2 Алгоритм вычислений и вопросы распараллеливания Совместное решение определяющих уравнений (3), (4), (8) в коде Flag-FS осуществляется посредством организации глобального итерационного процесса, на каждой итерации которого решаются линейные уравнения относительной приращений скорости, давления и параметров турбулентности для нахождения уточненных значений данных величин. «Перевязка» прираще- ний скорости и давления осуществляется по методу SIMPLEC. Поскольку для корректного решения уравнения переноса маркер-функции (3) требуется бездивергентное поле скорости, обновление поля маркер-функции проводится не на каждой глобальной итерации по решению уравнений гидродинамики (4) и (8), а с некоторым шагом, достаточным для получения поля скорости, близкого к дивергентному (как правило, шаг со- ставлял десяток итераций, а общее количество глобальных итераций доходило до 100). Исходя из того, что ограничения на величину шага по времени для уравнений гидродинамики, как пра- вило, мягче, чем для уравнения (3), была реализована методика, позволяющая, в целях эконо- мии вычислительных ресурсов, выполнять дробные шаги по времени для уравнения (3) в пре- делах одного шага по времени для решения уравнений гидродинамики [10]. На каждом из этих дробных шагов дискретный аналог уравнения переноса маркер-функции – уравнение (11) – ре- шается по методу установления (как правило, требовалось 3-5 итераций). Для решения систем линейных уравнений в коде Flag-FS используются итерационные ме- тоды на основе подпространств Крылова. Уравнение Пуассона для приращения давления в ме- тоде SIMPLEC решается с помощью метода сопряженных градиентов (CG) с неполным разло- жением Холецкого в качестве предобуславливателя (cм. например, [19]). Для решения осталь- ных уравнений используется вариант метода би-сопряженных градиентов (BiCGStab, [20]) с предобуславливателем SGS (Symmetric Gauss-Seidel). Распараллеливание вычислений производится по технологии “domain decomposition”, в рамках которой расчетная область разбивается на блоки примерно равного размера, и для вы- полнения вычислений внутри каждого блока создается отдельный процесс, оперирующий сво- им набором данных в оперативной памяти и выполняющийся на отдельном процессорном ядре. Для обеспечения связности численного решения производятся обмены данными между блока- ми (стыковка) с удовлетворением требованию «прозрачности» границ, т.е. отсутствия влияния разбиения расчетной области на получаемое решение. В коде Flag-S разбиение расчетной стеки проводится на соприкасающиеся (не накладыва- ющиеся) блоки. Связанность боков обеспечивается созданием в каждом из них ряда «загранич- ных» ячеек, повторяющих приграничные ячейки смежного блока (см. рис. 1). Для прозрачности границ требуется, чтобы потоки величин через грани ячеек на межблочных границах были одинаковыми для смежных ячеек, находящихся в разных блоках. Поскольку в рамках реализо- ванных в коде Flag-S численных схем второго порядка для вычисления потока величины через смежную грань двух ячеек требуются лишь значения самой величины и ее градиента в центрах этих ячеек, прозрачность границ обеспечивается копированием приграничных значений вели- чин и их градиентов. Учитывая, что при вычислении градиента величины в ячейке сетки ис- пользуются значения этой величины из соседних ячеек, стыковка блоков осуществляется в два 365 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt этапа: сначала блоки обмениваются значениями величин, по которым в каждом из блоков вы- числяются их градиенты, а затем передаются вычисленные градиенты. Обмены данных между блоками производятся по технологии MPI. Для уменьшения числа обращений к сравнительно медленным функциям MPI, все данные по каждой стыковке пере- даются единым массивом; при этом для пересылки используются неблокирующие процедуры MPI_ISEND и MPI_IRECV, так что каждая стыковка обрабатывается по мере готовности дан- ных в стыкуемых блоках, независимо от других стыковок. Блок 1 Блок 2 Блок 1 «Заграничные» Блок 2 ячейки передачи данных Рис. 1. Передача данных при стыковке блоков 3. Расчеты натекания потока на треугольное препятствие В качестве одного из приложений разработанного кода выступала тестовая задача о нате- кании потока, вызванного обрушением дамбы, на треугольное препятствие. Постановка задачи соответствовала эксперименту [1]. Схема экспериментальной установки приведена на рис. 2. Перед началом эксперимента вода удерживалась вертикальной перегородкой («дамбой») в ле- вой части длинного канала прямоугольного сечения. Затем перегородка резко убиралась, вода начинала распространяться вправо, натекая на препятствие и перетекая через него. Боковые стенки канала были прозрачными; положения свободной поверхности в различные моменты времени фиксировались при помощи фотокамеры, расположенной сбоку. Число Рейнольдса, определенное по начальной высоте столба жидкости H и характерной скорости потока (2gH)1/2, составляло 1.6·105. Представляемые ниже расчеты, проведенные в двумерной и трехмерной постановках, были направлены на исследование влияния вязких эффектов вблизи дна и боковых стенок канала на эволюционирующую картину течения. В этих расчетах пользовательский параметр dcell, входя- щий в модель искусственной диффузии маркер-функции у стенок (10), полагался равным поло- вине характерного размера ячеек вдоль стенки, а коэффициент искусственной диффузии max = ·(gH3)1/2, где значение коэффициента =10-4 подобрано эмпирически. 3.1 Результаты расчетов в двумерной постановке Двумерные расчеты были проведены в двух вариантах: с постановкой условий прилипания на стенках (и, соответственно, с учетом пристеночных вязких эффектов) и с использованием условий проскальзывания. В первом случае сетка имела сгущение к нижней границе, обеспечи- вающее значения нормализованного расстояния от стенки до центров пристенной ячейки y1  y1  w  , не превосходящие единицы, и состояла из 25 000 прямоугольных ячеек (с не- большой скошенностью в области над препятствием). Во втором случае использовалась сетка из 16 000 ячеек, без сгущения. 366 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt торцевая стенка боковые стенки Рис. 2. Схема экспериментальной установки работы [1]. Все размеры даны в сантиметрах Форма свободной поверхности, полученная в расчетах с разными граничными условиями, иллюстрируется на рис. 3. Результаты расчетов наложены на экспериментальные фотографии из работы [1] для двух моментов времени, отсчитываемых от устранения удерживавшей воду перегородки. Как видно из рисунка, форма свободной поверхности в расчете с условиями про- скальзывания (рис. 3 а, в) весьма сильно отличается от наблюдаемой в эксперименте. Учет эф- фектов придонного трения приводит к существенному улучшению качества предсказания тече- ния. Основное отличие между двумя решениями состоит в появлении отрывных зон в расчете с учетом придонного трения. Рис. 3 б, г показывает, что сначала одна (t = 3.0 c), а затем две (t = 3.7 c) крупные отрывные зоны возникают перед препятствием и на его наветренной сто- роне. Эти зоны приводят к возникновению «холмов» на свободной поверхности, видимых и на экспериментальных фотографиях. а) б) в) г) Рис. 3. Влияние придонного трения на расчетную форму свободной поверхности: (а,в) – расчет с услови- ем проскальзывания, (б,г) – расчет с условием прилипания (также показаны отрывные зоны). Данные расчетов нанесены поверх экспериментальных фотографий [1] для (а,б) t=3.0 с и (в,г) t=3.7 с 3.2 Трехмерные расчеты Трехмерные расчеты проводились с целью анализа влияния вязких эффектов у боковых стенок канала на структуру течения и форму свободной поверхности. Рассматриваемая конфи- гурация обладает симметрией относительно средней вертикальной плоскости. Это позволило использовать расчетную область, соответствующую половине экспериментальной. Расчетная сетка состояла из одного миллиона шестигранных ячеек и была сгущена к нижней и к боковой стенкам; степень сгущения к обеим границам была такой же, как в двумерном расчете. В расчете были задействованы 62 вычислительных ядра кластера с процессорами типа AMD Opteron 6300 (ускорение примерно в 45 раз по сравнению с однопроцессорным вариан- том); общее время счета – около недели (2000 шагов по времени). 367 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt На рис. 4 представлена расчетная форма свободной поверхности в те же моменты времени, что и на рис. 3, а на рис. 5 дана векторная картина скоростного поля в момент времени 3.0 с. Видно, что, как и в двумерном расчете, возникают сначала одна, а потом две отрывные зоны, приводящие к появлению «холмов» на свободной поверхности. Примерно на половине полу- ширины канала, со стороны плоскости симметрии, высота подъема свободной поверхности на «холмах» практически постоянна поперек канала. Однако по мере приближения к боковой стенке свободная поверхность существенно отклоняется от двумерной формы, а «холмы» прак- тически вырождаются. а) б) плоскость плоскость симметрии симметрии Рис. 4. Свободная поверхность по результатам трехмерных расчетов для моментов времени (а) t=3.0 с и (б) t=3.7 с. Также показаны вихревые линии (пунктир), идущие из центров отрывных зон в плоскости симметрии Трехмерные эффекты еще в большей степени проявляются в конфигурации отрывных зон. Исходящие из центров отрывных зон в плоскости симметрии вихревые линии (обозначены пунктиром на рис. 4) на большей части своей длины отклоняются от прямой линии, которая имела бы место в двумерном течении, а вблизи боковой стенки канала «поворачивают» наверх и выходят на поверхность жидкости. Свидетельство выхода «оси» отрывной зоны на свобод- ную поверхность можно видеть и в картине течения на свободной поверхности (рис. 5), где от- четливо видна область рециркуляционного движения вблизи боковой стенки (над наветренной стороной препятствия). Рис. 5. Расчетное поле скорости на свободной поверхности и в плоскости симметрии, t = 3.0 с На рис. 6 расчетная форма свободной поверхности (вид сбоку) сопоставляется с экспе- риментальными фотографиями [1]. Здесь свободная поверхность представляется не линией (как в двумерном случае), а имеет вид ленты, изменения в видимой толщине которой отражают за- висимость высоты подъема жидкости от трансверсальной координаты. Примечательно, что 368 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt сходственные «ленты» видны и на экспериментальных фотографиях (темные области вдоль свободной поверхности). Сопоставление рис. 3 и 6 позволяет заключить, что результаты трех- мерного расчета намного лучше согласуются с экспериментальными данными. а) б) Рис. 6. Результаты трехмерных расчетов течения, соответствующего эксперименту [1] (мгновенные формы свободной поверхности) в сопоставлении с экспериментальными фотографиями для моментов времени 3.0 с (а) и 3.7 с (б) 4. Расчеты натекания потока на множественные (однорядные и дву- рядные) препятствия Другой пример приложения разработанного кода – это модельная задача о натекании пото- ка воды, возникшего в результате разрушения дамбы, на препятствия в форме параллелепипе- дов, выстроенных поперек потока в один и два бесконечных ряда (во втором случае – в шах- матном порядке, см. рис. 7). Данная модель имитирует постройки в прибрежной зоне. Основа- ния препятствий имели размеры 16х16 см, расстояние между рядами препятствий равнялось 16 см, шаг между препятствиями одного ряда – 32 см. Число Рейнольдса Re=l (2gH3)1/2 /l со- ставляло 1.8·106. В расчетах с одним рядом отсутствовал передний ряд (кубических) препят- ствий. 55 90 120 116 16 Рис. 7. Постановка задачи о натекании воды на ряд препятствий. Все размеры даны в сантиметрах С учетом присущей задаче симметрии, расчетная область в поперечном направлении рас- полагалась от середины низкого препятствия до середины соседнего с ним высокого препят- 369 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt ствия (на продольных вертикальных границах области ставились условия симметрии). Для обе- их конфигураций расчетные сетки состояли из примерно полумиллиона шестигранных ячеек и имели одинаковое сгущение ко дну, обеспечивающее значения y+, не превосходящие единицы. Результаты расчетов для вариантов с одним и двумя рядами препятствий приведены на рис. 8 (для удобства восприятия визуализированных картин течения расчетные поля размножены с учетом симметрии, присущей задаче). Рис. 8. Положения свободной поверхности и поля давления на лицевой стороне высокого препятствия в моменты времени 0.8 c (верхний ряд), 1.2 с (средний ряд) и 2.0 с (нижний ряд): результаты расчетов в отсутствие первого ряда препятствий (слева) и при его наличии (справа) 370 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Как видно из рисунка, характер течений, рассчитанных для вариантов с одним и двумя ря- дами препятствий, существенно различен. В отсутствие первого ряда поток ударяется в ниж- нюю часть высоких препятствий, часть жидкости движется вверх вдоль стенок, а затем опроки- дывается назад, навстречу потоку. При наличии первого ряда из кубических препятствий поток сначала ударяется о них. Часть потока «перелетает» через низкие препятсвия и ударяется о вто- рой ряд препятствий, достигая середины их высоты. Различия в картине течений заметно про- являются в поле давления на наветренной стороне высоких препятствий. В отсутствие первого ряда область максимальных значения давления располагается вблизи дна. При наличии первого ряда зона высокого давления несколько уменьшается, однако перемещается выше, создавая, таким образом, больший опрокидывающий момент (12.1 Н·м против 4.8 Н·м при t=0.8 c). Оче- видно, последнее представляет большую опасность с точки зрения потенциального разрушения прибрежных зданий под действием мощной набегающей волны. 5. Выводы На базе кода внутреннего пользования Flag-S создана специализированная версия, позво- ляющая на основе метода VOF проводить параллельные расчеты трехмерных и двумерных не- стационарных течений жидкости со свободной поверхностью. Разработан (и реализован) ряд оригинальных составляющих метода VOF для обеспечения его корректной работы в сложных с вычислительной точки зрения случаях. Произведен поиск оптимальных численных схем в плане чувствительности результатов расчетов к густоте используемой расчетной сетки и шагам по времени. Разработанный программный код применен к решению тестовых и модельных за- дач о нестационарном натекании потока воды препятствия различной формы, в частности, на одиночное треугольное препятствие и на множественные (однорядные и двурядные) препят- ствия в форме параллелепипедов. На примере задачи с треугольным препятствием показано, что пристенные вязкие эффекты могут оказывать весьма сильное влияние на картину течения через формирование крупных от- рывных зон, приводящих к возникновению «холмов» на свободной поверхности. Установлено, что вязкие эффекты вблизи боковых стенок канала делают картину течения существенно трех- мерной. Результаты трехмерного расчета хорошо согласуются с экспериментальными данными. Для натекания потока на множественные препятствия получен интересный для практики вывод о том, что наличие первого ряда препятствий может привести к увеличению опрокиды- вающего момента, действующего на препятствия второго ряда. Литература 1. Soares-Frazao S. Experiments of dam-break wave over a triangular bottom sill // Journal of Hy- draulic 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. 5. Hirt C.W., Nichols B.D. Volume of fluid (VOF) method for the dynamics of free boundaries // Jour- nal of computational physics. 1981. Vol. 39. P. 201-226. 6. Vandoormaal J.P., Raithby G.D. Enhancements of the SIMPLE method for predicting incompress- ible 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 trail- ing edge separation // AIAA Journal. 1983. Vol. 21. P. 1525-1532. 371 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 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- Turbu- lence Models Aimed at Grid Sensitivity Reduction // Proceedings of the European Conference for Aerospace Sciences (EUCASS 2005), 4–7 July 2005, Moscow. CD-ROM, ID 2.09.06. 7 p. 10. Храбрый А.И., Зайцев Д.К., Смирнов Е.М. Численное моделирование течений со свободной поверхностью на основе метода VOF // Труды Крыловского государственного научного центра. 2013. Вып. 78. С. 53–64. 11. Gaskell P.H., Lau A.K.C . Curvature-compensated convective-transport: SMART, A new bound- edness-preserving transport algorithm // Int. J. Numer. Methods Fluids. 1988.Vol. 8. Issue 6. P. 617-641. 12. Muzaferija S., Peric M., Sames P., Schelin T. A two-fluid Navier-Stokes solver to simulate water entry // Proceedings of the 22 Symposium 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 in- terface 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 Sev- eral High-Resolution Finite Volume Schemes: Test Computations / In: Computational Fluid Dy- namics 2010. Ed.: A.Kuzmin. New-York: Springer, 2011. 954 p. P. 535-540. 16. Khrabry A.I., Smirnov E.M., Zaytsev D.K. Free surface flow computations using the M-CICSAM scheme added with a sharpening procedure // Proceedings of the 6th European Congress on Com- putational Methods in Applied Sciences and Engineering (ECCOMAS 2012), Vienna, Austria, 10-14 September 2012. CD-ROM. ISBN: 978-3-9502481-9-7. 2 p. 17. Зайцев Д.К., Смирнов Е.М., Храбрый А.И. Расчет течений со свободными поверхностями: влияние схемных факторов и модели турбулентности // Труды 14-й междунар. конф. «Су- первычисления и математическое моделирование», Саров, 1-5 октября 2012. С. 282-292. 18. Khrabry A., Smirnov E., Zaytsev D., Goryachev V. Numerical study of separation phenomena in the dam-break flow interacting with a triangular obstacle // Proceedings of the 16th International Conference on Fluid Flow Technologies (CMFF’15), Budapest, Hungary, September 1-4 2015. CD-ROM, CMFF15-144.pdf. 8 p. 19. Голуб Дж., Ван Лоун Ч. Матричные Вычисления. М.: Мир, 1999. 553 с. 20. Van Der Vorst H.A. Bi-CGSTAB: A fast and smoothly converging variant of Bi CG for the solu- tion of nonsymmetric linear systems // SIAM Journal on Scientific and Statistical Computing. 1992. Vol. 13. P. 631-644. 372 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Numerical simulation of 3D unsteady free-surface flows: development and application of a specialized parallel code A.I. Khrabry, D.K. Zaytsev, E.M. Smirnov Saint-Petersburg Polytechnic University 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 computa- tional 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. Keywords: numerical simulation, parallel code, domain decomposition, free-surface flow, flow-obstacle interaction, turbulent flows, VOF method, viscous separation. References 1. Soares-Frazao S. Experiments of dam-break wave over a triangular bottom sill // Journal of Hy- draulic 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. Fedotova Z.I., Chubarov L.B. Chislennoe modelirovanie nakata tsunami [Numerical modeling of tsunami rolling] // Trudy Mezhdunarodnoy konferentsii RDAMM-2001 [Proceedings of Interna- tional Conference RDAMM–2001], 2001. Vol. 6. Part 2. Special Issue. P. 380-396. 5. Hirt C.W., Nichols B.D. Volume of fluid (VOF) method for the dynamics of free boundaries // Jour- nal of computational physics. 1981. Vol. 39. P. 201-226. 6. Vandoormaal J.P., Raithby G.D. Enhancements of the SIMPLE method for predicting incompress- ible 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 trail- ing 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- Turbu- lence Models Aimed at Grid Sensitivity Reduction // Proceedings of the European Conference for Aerospace Sciences (EUCASS 2005), 4–7 July 2005, Moscow. CD-ROM, ID 2.09.06. 7 p. 10. Khrabry A.I., Smirnov E.M., Zaytsev D.K. Chislennoe modelirovanie techeniy so svobodnoy pov- erkhnost'yu na osnove metoda VOF [Numerical simulation of free surface flows on the base of the VOF method] // Trudy Krylovskogo gosudarstvennogo nauchnogo tsentra [Proceedings of Krylov state research centre]. 2013. Vol. 78. P. 53–64. 11. Gaskell P.H., Lau A.K.C . Curvature-compensated convective-transport: SMART, A new bound- edness-preserving transport algorithm // Int. J. Numer. Methods Fluids. 1988.Vol. 8. Issue 6. P. 617-641. 373 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 12. Muzaferija S., Peric M., Sames P., Schelin T. A two-fluid Navier-Stokes solver to simulate water entry // Proceedings of the 22 Symposium 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 in- terface 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 Sev- eral High-Resolution Finite Volume Schemes: Test Computations / In: Computational Fluid Dy- namics 2010. Ed.: A.Kuzmin. New-York: Springer, 2011. 954 p. P. 535-540. 16. Khrabry A.I., Smirnov E.M., Zaytsev D.K. Free surface flow computations using the M-CICSAM scheme added with a sharpening procedure // Proceedings of the 6th European Congress on Com- putational Methods in Applied Sciences and Engineering (ECCOMAS 2012), Vienna, Austria, 10-14 September 2012. CD-ROM. ISBN: 978-3-9502481-9-7. 2 p. 17. Zaytsev D.K., Smirnov E.M., Khrabry A.I. Numerical simulation of free surface flows: influence of scheme factors and turbulence model // Proc. 14th International Conference “Supercomputing and mathematical modeling”, Sarov, Russia. 2012. P. 282-292. 18. Khrabry A., Smirnov E., Zaytsev D., Goryachev V. Numerical study of separation phenomena in the dam-break flow interacting with a triangular obstacle // Proceedings of the 16th International Conference on Fluid Flow Technologies (CMFF’15), Budapest, Hungary, September 1-4 2015. CD-ROM, CMFF15-144.pdf. 8 p. 19. Golub G., Van Loan С. Matrichnye Vychisleniya [Matrix Computations]. M.: Mir [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 solu- tion of nonsymmetric linear systems // SIAM Journal on Scientific and Statistical Computing. 1992. Vol. 13. P. 631-644. 374