<!DOCTYPE article PUBLIC "-//NLM//DTD JATS (Z39.96) Journal Archiving and Interchange DTD v1.0 20120330//EN" "JATS-archivearticle1.dtd">
<article xmlns:xlink="http://www.w3.org/1999/xlink">
  <front>
    <journal-meta>
      <journal-title-group>
        <journal-title>EUCASS</journal-title>
      </journal-title-group>
    </journal-meta>
    <article-meta>
      <title-group>
        <article-title>Санкт-Петербургский политехнический университет Петра Великого</article-title>
      </title-group>
      <pub-date>
        <year>2005</year>
      </pub-date>
      <volume>4</volume>
      <fpage>4</fpage>
      <lpage>7</lpage>
      <abstract>
        <p>Излагается вычислительная методика для моделирования течений жидкости со свободной поверхностью, с акцентом на оригинальные элементы, повышающие точность решения и обеспечивающие корректную работу метода VOF (Volume-Of-Fluid) в особо сложных с вычислительной точки зрения случаях. Методика реализована в специализированном трехмерном коде, оперирующем с неструктурированными расчетными сетками. Параллелизация вычислений осуществлена по технологии “domain decomposition”. Даются примеры приложения разработанного кода к решению двумерных и трехмерных задач о нестационарном натекании турбулентного потока воды на различные препятствия. Результаты расчетов сопоставляются с экспериментальными данными. Ключевые слова: численное моделирование, параллельный код, domain decomposition, течение со свободной поверхностью, натекание на препятствие, турбулентные течения, метод VOF, вязкий отрыв.</p>
      </abstract>
    </article-meta>
  </front>
  <body>
    <sec id="sec-1">
      <title>-</title>
      <p>Численное моделирование трехмерных нестационарных
течений со свободной поверхностью:
итоги разработки и примеры приложения
специализированного параллельного кода
щее время является метод Volume-of-Fluid (VOF) [5], позволяющий проводить расчеты течений с
сильной деформацией свободной поверхности (включая опрокидывание волны), обеспечивая при
этом высокую детализацию структуры эволюционирующего потока. Сопоставление методик,
основанных на приближении «мелкой воды» и на полных уравнениях гидродинамики,
проводится, например, в работе [2], где рассматривается натекания потока на трапециевидное препятствие.
Результаты данной работы, полученные на основе полных уравнений, свидетельствуют, в
частности, о том, что с наветренной стороны препятствия может формироваться отрывная зона,
вызванная ростом давления вдоль основного потока при его натекании на препятствие. Модель «мелкой
воды» в принципе не позволяет воспроизвести этот эффект.</p>
      <p>Применение численного моделирования для детального анализа столь сложного течения как
набегание потока на препятствие предполагает проверку влияния схемных факторов на
получаемое решение. Причем это относится как к задаче расчета формы свободной поверхности,
так и к воспроизведению эффектов пристеночного трения. Для рассматриваемых течений
характерны высокие числа Рейнольдса, и обычно для определения величины трения на стенках
применяется метод пристенных функций (вкупе с «высокорейнольдсовыми» расчетными
сетками), а нередко расчеты проводятся и вовсе с заданием условия проскальзывания на стенках.
Известно, однако, что пристенные функции могут плохо работать в условиях отрыва
пограничного слоя, вызванного «неблагоприятным» градиентом давления.</p>
      <p>Получение решения, свободного от влияния схемных факторов, требует весьма густых
расчетных сеток в целом по расчетной области с дополнительным сильным сгущением к твердым
границам (так называемых «низкорейнольдсовых» сеток). Применительно к нестационарным
трехмерным расчетам, это требование сопряжено с особо высокой затратностью
вычислительных ресурсов, что до недавнего времени делало такого рода расчеты невозможными. Сегодня,
благодаря распространению мощных многопроцессорных систем, проведение аккуратных
расчетов трехмерных нестационарных турбулентных течений со свободной поверхностью
становится все более реализуемым. Соответственно, возникает потребность в разработке и
апробации эффективных вычислительных средств, позволяющих достичь высокого качества
пространственного и временного разрешения течения.</p>
      <p>Настоящая работа направлена на разработку основанного на методе VOF численного
алгоритма для расчета нестационарных трехмерных течений жидкости со свободной поверхностью,
его программную реализацию в параллельном коде и применение к решению модельных задач
практической направленности о нестационарном натекании потока воды на различные
препятствия. Излагается вычислительная методика с акцентом на оригинальные элементы,
повышающие точность решения. Дается описание вычислительных приемов, обеспечивающих
корректную работу метода в особо сложных с вычислительной точки зрения случаях.
2. Численная методика и программная реализация</p>
      <p>Реализация и тестирование численных схем для решения системы уравнений метода VOF
проводились в программной среде гидродинамического кода внутреннего пользования Flag-S,
разрабатываемого сотрудникам кафедры гидроаэродинамики СПбПУ и исходно не
предназначавшегося для расчета течений со свободной поверхностью. В итоге была создана версия кода
для расчета течений данного класса, названная Flag-FS (Flag – Free Surface). Новая версия
унаследовала от кода Flag-S следующие основные свойства и элементы. Код оперирует
неструктурированными расчетными сетками, состоящими из полиэдральных ячеек, что позволяет
проводить расчеты течений в областях сложной геометрии. Дискретизация решаемых уравнений
выполняется по методу конечных объемов (МКО), связывающему изменение той или иной
физической величины в центре ячейки расчетной сетки с потоками этой величины через грани
ячейки. Вычисление конвективных и диффузионных потоков основывается на схемах второго
порядка, оперирующих значениями самой величины и ее градиента в центрах смежных по грани
ячеек. «Перевязывание» скорости и давления в итерационном процессе, реализуемом на
каждом шаге по времени, осуществляется по методу SIMPLEC [6]. Для предотвращения
пространственных осцилляций в поле течения применяется коррекция Рхи-Чоу [7]. Реализованы
несколько моделей турбулентности, включая используемую в представленных ниже расчетах SST
модель Ментера [8] с обобщенными пристенными функциями [9].
2.1 Определяющие уравнения и их дискретные аналоги</p>
      <p>В методе VOF для определения положения свободной поверхности вводится специальная
маркер-функция C, представляющая собой объемную долю жидкости (газу соответствует
значение C=0, жидкости – C=1). При расчетах по данному методу движущаяся по сетке межфазная
граница становится «размытой», приобретая вид переходного слоя толщиной в несколько
расчетных ячеек, в котором значений маркер-функции меняются от 0 до 1. Положение свободной
поверхности обычно определяется изоповерхностью С=0.5.</p>
      <p>В так называемой одножидкостной формулировке метода VOF уравнения гидродинамики
решаются «насквозь», как для единой среды с переменными материальными свойствами,
выражаемыми через величину C по следующим соотношениям:</p>
      <p>Здесь  – плотность среды,  – динамическая вязкость среды, индексы «l» и «g» относятся
к жидкости и газу соответственно.</p>
      <p>Поскольку объемная доля жидкости сохраняется вдоль траекторий, ее динамика может

быть описана уравнением конвективного переноса ( v – скорость потока):
Однако, это уравнение не является консервативным и не подходит для использования в рамках
концепции МКО. Для реализации в коде Flag-FS уравнение было переписано в консервативной
форме (3), одновременно с записью уравнения неразрывности в форме (4), предполагающей
несжимаемость жидкости и газа:
  Cl  1  Сg
  Cl  1  Сg
dC
dt</p>
      <p>t

C </p>
      <p> v C  0
Дискретные конечно-объемные аналоги данных уравнений имеют следующий вид:
Здесь индекс P относится к центру ячейки, V – ее объем, суммирование проводится по всем
 
граням ячейки, индекс f относится к центру грани, Ff  S f n f  v f – объемный расход через

грань, Sf – площадь грани, n f – внешняя нормаль к ней. Отметим, что форма (5) гарантирует
сохранение общего количества маркер-функции в расчетной области в целом, таким образом,
обеспечивая сохранение объемов, занимаемых жидкостью и газом.</p>
      <p>
        По причинам, детально изложенным в работе [
        <xref ref-type="bibr" rid="ref5">10</xref>
        ], уравнение баланса импульса было
тождественно преобразовано из традиционной консервативной формы к форме (7) с оригинальным
представлением конвективных членов.
      </p>
      <p>
 v    vv v   v  p    τ  g (7)
Здесь p – давление, g – ускорение силы тяжести, τ    t v  v  – тензор вязких

напряжений, t – турбулентная вязкость (дается моделью турбулентности).</p>
      <p>Форма (7) может быть обобщена на случай уравнений переноса (8) любой физической
величины  – в ее роли может выступать компонента скорости или параметр турбулентности:
где RHS – это сумма источниковых и диффузионных членов.</p>
      <p>Дискретный аналог уравнения (8) имеет следующий вид:
 p V   f  f Ff  p  f Ff  RHS
p t</p>
      <p>f f
Данный способ дискретизации конвективных слагаемых позволяет избежать нефизичных
осцилляций в решении, вызываемых большим перепадом значений плотности вблизи
межфазной границы.</p>
      <p>
        Фигурирующие в дискретных аналогах определяющих уравнений значения величин на
гранях ячеек находятся с помощью процедуры реконструкции, которая задействует значения из
центров ближайших ячеек. Для компонент скорости и параметров турбулентности в коде
Flag-FS используются противопоточные схемы второго порядка. Однако такие схемы не
подходят для вычисления значений маркер-функции Cf на гранях ячеек, так как вносимые ими
эффекты численной диффузии приводят к размытию межфазной границы на множество ячеек
расчетной сетки. В литературе предложен ряд специализированных, т.н. «сжимающих»
численных схем для нахождения значений Cf. Большинство таких схем удовлетворяют т.наз. CBC
критерию [11], который гарантирует локальную ограниченность получаемого решения в случае
одномерной постановки и равномерной расчетной сетки. Однако свойства схем в многомерном
случае и при расчетах на неравномерных/скошенных сетках каждый раз требуют отдельного
изучения. В частности, результаты детального сопоставления возможностей двух популярных
схем HRIC [12] и CICSAM [13], а также предложенной позднее схемы M-CICSAM [
        <xref ref-type="bibr" rid="ref6">14</xref>
        ],
представлены в работе [
        <xref ref-type="bibr" rid="ref7">15</xref>
        ]. В проведенных тестах схема M-CICSAM продемонстрировала явное
превосходство над другими двумя схемами: меньшую зависимость качества результатов расчетов от
густоты и скошенности расчетной сетки, величины шага по времени. С учетом результатов
работы [
        <xref ref-type="bibr" rid="ref7">15</xref>
        ], для реализации в коде Flag-FS была выбрана схема M-CICSAM в сочетании со схемой
Кранка-Николсон для аппроксимации по времени. Такая комбинация хорошо работает в
широком диапазоне чисел Куранта, практически не искажая межфазную границу и не допуская
сколько-нибудь заметных «вылетов» значений маркер-функции за пределы диапазона [0, 1]. Во
избежание получения нефизичных (отрицательных) значений плотности и вязкости, в
выражениях (1) и (2) упомянутые «вылеты» (которые, как правило, не превосходили 0.001) обрезались.
      </p>
      <p>
        Кроме того, для устранения небольшого остаточного размытия межфазной границы
(которое, как показано в работе [16], в отдельных случаях может приводить к существенному
искажению поля течения), в коде Flag-FS была реализована оригинальная методика «обострения»
фронта [16]. Отметим, что схема Кранка-Николсон использовалась для аппроксимации и
остальных уравнений, что позволило обеспечить достаточно высокую точность решения при
числах Куранта, достигающих 0.8. Именно с таким значением числа Куранта проводились
представленные ниже расчеты (шаг по времени подирался автоматически). Вопросы
чувствительности решения к густоте расчетной сетки исследовались в работе [
        <xref ref-type="bibr" rid="ref8">17</xref>
        ] на тестовой задаче
натекания потока на квадратное препятствие.
      </p>
      <p>
        Методические расчеты показали также, что при использовании «низкорейнольдсовых»
сеток для решения задач о растекании жидкости вдоль сухой твердой стенки, вблизи последней
формируется тонкая нефизичная газовая прослойка, приводящая к существенному искажению
величины трения на стенке [
        <xref ref-type="bibr" rid="ref9">18</xref>
        ]. Для устранения описанного артефакта в коде Flag-FS
реализована оригинальная методика, вводящая искусственную диффузию маркер-функции, которая
действует только вблизи твердой стенки, по направлению к ней:
С
t
   C v  max nb   nCb , 0 .
      </p>
      <p>

Здесь nb  d – орт направления, перпендикулярного ближайшей твердой границе,  –
коэффициент искусственной диффузии, определяемый по формуле:
(9)
(10)
dcell  d
dcell
,d  dcell
,d  dcell</p>
      <p>,
где d – расстояние до ближайшей стенки, dcell и max – пользовательские параметры. Параметр
dcell определяет расстояние от стенки, на котором действует искусственная диффузия, max
задает максимальное значение коэффициента искусственной диффузии, достигаемое на стенке.
Дискретный аналог уравнения (10) имеет следующий вид:</p>
      <p> 
Сp V  C f Ff  max  f C f  d n f  d S f , 0
t f  f 
(11)
2.2 Алгоритм вычислений и вопросы распараллеливания</p>
      <p>Совместное решение определяющих уравнений (3), (4), (8) в коде Flag-FS осуществляется
посредством организации глобального итерационного процесса, на каждой итерации которого
решаются линейные уравнения относительной приращений скорости, давления и параметров
турбулентности для нахождения уточненных значений данных величин. «Перевязка»
приращений скорости и давления осуществляется по методу SIMPLEC.</p>
      <p>
        Поскольку для корректного решения уравнения переноса маркер-функции (3) требуется
бездивергентное поле скорости, обновление поля маркер-функции проводится не на каждой
глобальной итерации по решению уравнений гидродинамики (4) и (8), а с некоторым шагом,
достаточным для получения поля скорости, близкого к дивергентному (как правило, шаг
составлял десяток итераций, а общее количество глобальных итераций доходило до 100). Исходя
из того, что ограничения на величину шага по времени для уравнений гидродинамики, как
правило, мягче, чем для уравнения (3), была реализована методика, позволяющая, в целях
экономии вычислительных ресурсов, выполнять дробные шаги по времени для уравнения (3) в
пределах одного шага по времени для решения уравнений гидродинамики [
        <xref ref-type="bibr" rid="ref5">10</xref>
        ]. На каждом из этих
дробных шагов дискретный аналог уравнения переноса маркер-функции – уравнение (11) –
решается по методу установления (как правило, требовалось 3-5 итераций).
      </p>
      <p>
        Для решения систем линейных уравнений в коде Flag-FS используются итерационные
методы на основе подпространств Крылова. Уравнение Пуассона для приращения давления в
методе SIMPLEC решается с помощью метода сопряженных градиентов (CG) с неполным
разложением Холецкого в качестве предобуславливателя (cм. например, [
        <xref ref-type="bibr" rid="ref10">19</xref>
        ]). Для решения
остальных уравнений используется вариант метода би-сопряженных градиентов (BiCGStab, [
        <xref ref-type="bibr" rid="ref11">20</xref>
        ]) с
предобуславливателем SGS (Symmetric Gauss-Seidel).
      </p>
      <p>Распараллеливание вычислений производится по технологии “domain decomposition”, в
рамках которой расчетная область разбивается на блоки примерно равного размера, и для
выполнения вычислений внутри каждого блока создается отдельный процесс, оперирующий
своим набором данных в оперативной памяти и выполняющийся на отдельном процессорном ядре.
Для обеспечения связности численного решения производятся обмены данными между
блоками (стыковка) с удовлетворением требованию «прозрачности» границ, т.е. отсутствия влияния
разбиения расчетной области на получаемое решение.</p>
      <p>В коде Flag-S разбиение расчетной стеки проводится на соприкасающиеся (не
накладывающиеся) блоки. Связанность боков обеспечивается созданием в каждом из них ряда
«заграничных» ячеек, повторяющих приграничные ячейки смежного блока (см. рис. 1). Для прозрачности
границ требуется, чтобы потоки величин через грани ячеек на межблочных границах были
одинаковыми для смежных ячеек, находящихся в разных блоках. Поскольку в рамках
реализованных в коде Flag-S численных схем второго порядка для вычисления потока величины через
смежную грань двух ячеек требуются лишь значения самой величины и ее градиента в центрах
этих ячеек, прозрачность границ обеспечивается копированием приграничных значений
величин и их градиентов. Учитывая, что при вычислении градиента величины в ячейке сетки
используются значения этой величины из соседних ячеек, стыковка блоков осуществляется в два
этапа: сначала блоки обмениваются значениями величин, по которым в каждом из блоков
вычисляются их градиенты, а затем передаются вычисленные градиенты.</p>
      <p>Обмены данных между блоками производятся по технологии MPI. Для уменьшения числа
обращений к сравнительно медленным функциям MPI, все данные по каждой стыковке
передаются единым массивом; при этом для пересылки используются неблокирующие процедуры
MPI_ISEND и MPI_IRECV, так что каждая стыковка обрабатывается по мере готовности
данных в стыкуемых блоках, независимо от других стыковок.</p>
      <p>Блок 1
Блок 2
Блок 1
«Заграничные»
ячейки</p>
      <p>Блок 2
передачи
данных
Рис. 1. Передача данных при стыковке блоков
3. Расчеты натекания потока на треугольное препятствие
В качестве одного из приложений разработанного кода выступала тестовая задача о
натекании потока, вызванного обрушением дамбы, на треугольное препятствие. Постановка задачи
соответствовала эксперименту [1]. Схема экспериментальной установки приведена на рис. 2.
Перед началом эксперимента вода удерживалась вертикальной перегородкой («дамбой») в
левой части длинного канала прямоугольного сечения. Затем перегородка резко убиралась, вода
начинала распространяться вправо, натекая на препятствие и перетекая через него. Боковые
стенки канала были прозрачными; положения свободной поверхности в различные моменты
времени фиксировались при помощи фотокамеры, расположенной сбоку. Число Рейнольдса,
определенное по начальной высоте столба жидкости H и характерной скорости потока (2gH)1/2,
составляло 1.6·105.</p>
      <p>Представляемые ниже расчеты, проведенные в двумерной и трехмерной постановках, были
направлены на исследование влияния вязких эффектов вблизи дна и боковых стенок канала на
эволюционирующую картину течения. В этих расчетах пользовательский параметр dcell,
входящий в модель искусственной диффузии маркер-функции у стенок (10), полагался равным
половине характерного размера ячеек вдоль стенки, а коэффициент искусственной диффузии
max = ·(gH3)1/2, где значение коэффициента =10-4 подобрано эмпирически.
3.1 Результаты расчетов в двумерной постановке</p>
      <p>Двумерные расчеты были проведены в двух вариантах: с постановкой условий прилипания
на стенках (и, соответственно, с учетом пристеночных вязких эффектов) и с использованием
условий проскальзывания. В первом случае сетка имела сгущение к нижней границе,
обеспечивающее значения нормализованного расстояния от стенки до центров пристенной ячейки
y1  y1 w  , не превосходящие единицы, и состояла из 25 000 прямоугольных ячеек (с
небольшой скошенностью в области над препятствием). Во втором случае использовалась сетка
из 16 000 ячеек, без сгущения.
торцевая стенка</p>
      <p>боковые стенки
Рис. 2. Схема экспериментальной установки работы [1]. Все размеры даны в сантиметрах
Форма свободной поверхности, полученная в расчетах с разными граничными условиями,
иллюстрируется на рис. 3. Результаты расчетов наложены на экспериментальные фотографии
из работы [1] для двух моментов времени, отсчитываемых от устранения удерживавшей воду
перегородки. Как видно из рисунка, форма свободной поверхности в расчете с условиями
проскальзывания (рис. 3 а, в) весьма сильно отличается от наблюдаемой в эксперименте. Учет
эффектов придонного трения приводит к существенному улучшению качества предсказания
течения. Основное отличие между двумя решениями состоит в появлении отрывных зон в расчете с
учетом придонного трения. Рис. 3 б, г показывает, что сначала одна (t = 3.0 c), а затем две
(t = 3.7 c) крупные отрывные зоны возникают перед препятствием и на его наветренной
стороне. Эти зоны приводят к возникновению «холмов» на свободной поверхности, видимых и на
экспериментальных фотографиях.
Рис. 3. Влияние придонного трения на расчетную форму свободной поверхности: (а,в) – расчет с
условием проскальзывания, (б,г) – расчет с условием прилипания (также показаны отрывные зоны).
Данные расчетов нанесены поверх экспериментальных фотографий [1] для (а,б) t=3.0 с и (в,г) t=3.7 с
3.2 Трехмерные расчеты</p>
      <p>Трехмерные расчеты проводились с целью анализа влияния вязких эффектов у боковых
стенок канала на структуру течения и форму свободной поверхности. Рассматриваемая
конфигурация обладает симметрией относительно средней вертикальной плоскости. Это позволило
использовать расчетную область, соответствующую половине экспериментальной. Расчетная
сетка состояла из одного миллиона шестигранных ячеек и была сгущена к нижней и к боковой
стенкам; степень сгущения к обеим границам была такой же, как в двумерном расчете.</p>
      <p>В расчете были задействованы 62 вычислительных ядра кластера с процессорами типа
AMD Opteron 6300 (ускорение примерно в 45 раз по сравнению с однопроцессорным
вариантом); общее время счета – около недели (2000 шагов по времени).
На рис. 4 представлена расчетная форма свободной поверхности в те же моменты времени,
что и на рис. 3, а на рис. 5 дана векторная картина скоростного поля в момент времени 3.0 с.
Видно, что, как и в двумерном расчете, возникают сначала одна, а потом две отрывные зоны,
приводящие к появлению «холмов» на свободной поверхности. Примерно на половине
полуширины канала, со стороны плоскости симметрии, высота подъема свободной поверхности на
«холмах» практически постоянна поперек канала. Однако по мере приближения к боковой
стенке свободная поверхность существенно отклоняется от двумерной формы, а «холмы»
практически вырождаются.
плоскость
симметрии
Рис. 4. Свободная поверхность по результатам трехмерных расчетов для моментов времени
(а) t=3.0 с и (б) t=3.7 с. Также показаны вихревые линии (пунктир),</p>
      <p>идущие из центров отрывных зон в плоскости симметрии
Трехмерные эффекты еще в большей степени проявляются в конфигурации отрывных зон.
Исходящие из центров отрывных зон в плоскости симметрии вихревые линии (обозначены
пунктиром на рис. 4) на большей части своей длины отклоняются от прямой линии, которая
имела бы место в двумерном течении, а вблизи боковой стенки канала «поворачивают» наверх
и выходят на поверхность жидкости. Свидетельство выхода «оси» отрывной зоны на
свободную поверхность можно видеть и в картине течения на свободной поверхности (рис. 5), где
отчетливо видна область рециркуляционного движения вблизи боковой стенки (над наветренной
стороной препятствия).</p>
      <p>Рис. 5. Расчетное поле скорости на свободной поверхности и в плоскости симметрии, t = 3.0 с
На рис. 6 расчетная форма свободной поверхности (вид сбоку) сопоставляется с
экспериментальными фотографиями [1]. Здесь свободная поверхность представляется не линией (как
в двумерном случае), а имеет вид ленты, изменения в видимой толщине которой отражают
зависимость высоты подъема жидкости от трансверсальной координаты. Примечательно, что
сходственные «ленты» видны и на экспериментальных фотографиях (темные области вдоль
свободной поверхности). Сопоставление рис. 3 и 6 позволяет заключить, что результаты
трехмерного расчета намного лучше согласуются с экспериментальными данными.
Рис. 6. Результаты трехмерных расчетов течения, соответствующего эксперименту [1]
(мгновенные формы свободной поверхности) в сопоставлении с экспериментальными фотографиями
для моментов времени 3.0 с (а) и 3.7 с (б)
4. Расчеты натекания потока на множественные (однорядные и
двурядные) препятствия</p>
      <p>Другой пример приложения разработанного кода – это модельная задача о натекании
потока воды, возникшего в результате разрушения дамбы, на препятствия в форме
параллелепипедов, выстроенных поперек потока в один и два бесконечных ряда (во втором случае – в
шахматном порядке, см. рис. 7). Данная модель имитирует постройки в прибрежной зоне.
Основания препятствий имели размеры 16х16 см, расстояние между рядами препятствий равнялось
16 см, шаг между препятствиями одного ряда – 32 см. Число Рейнольдса Re=l (2gH3)1/2 /l
составляло 1.8·106. В расчетах с одним рядом отсутствовал передний ряд (кубических)
препятствий.</p>
      <p>Рис. 7. Постановка задачи о натекании воды на ряд препятствий. Все размеры даны в сантиметрах
С учетом присущей задаче симметрии, расчетная область в поперечном направлении
располагалась от середины низкого препятствия до середины соседнего с ним высокого
препят120</p>
      <p>55
ствия (на продольных вертикальных границах области ставились условия симметрии). Для
обеих конфигураций расчетные сетки состояли из примерно полумиллиона шестигранных ячеек и
имели одинаковое сгущение ко дну, обеспечивающее значения y+, не превосходящие единицы.
Результаты расчетов для вариантов с одним и двумя рядами препятствий приведены на рис. 8
(для удобства восприятия визуализированных картин течения расчетные поля размножены с
учетом симметрии, присущей задаче).</p>
      <p>Рис. 8. Положения свободной поверхности и поля давления на лицевой стороне высокого
препятствия в моменты времени 0.8 c (верхний ряд), 1.2 с (средний ряд) и 2.0 с (нижний ряд):
результаты расчетов в отсутствие первого ряда препятствий (слева) и при его наличии (справа)
Как видно из рисунка, характер течений, рассчитанных для вариантов с одним и двумя
рядами препятствий, существенно различен. В отсутствие первого ряда поток ударяется в
нижнюю часть высоких препятствий, часть жидкости движется вверх вдоль стенок, а затем
опрокидывается назад, навстречу потоку. При наличии первого ряда из кубических препятствий поток
сначала ударяется о них. Часть потока «перелетает» через низкие препятсвия и ударяется о
второй ряд препятствий, достигая середины их высоты. Различия в картине течений заметно
проявляются в поле давления на наветренной стороне высоких препятствий. В отсутствие первого
ряда область максимальных значения давления располагается вблизи дна. При наличии первого
ряда зона высокого давления несколько уменьшается, однако перемещается выше, создавая,
таким образом, больший опрокидывающий момент (12.1 Н·м против 4.8 Н·м при t=0.8 c).
Очевидно, последнее представляет большую опасность с точки зрения потенциального разрушения
прибрежных зданий под действием мощной набегающей волны.
5. Выводы</p>
      <p>На базе кода внутреннего пользования Flag-S создана специализированная версия,
позволяющая на основе метода VOF проводить параллельные расчеты трехмерных и двумерных
нестационарных течений жидкости со свободной поверхностью. Разработан (и реализован) ряд
оригинальных составляющих метода VOF для обеспечения его корректной работы в сложных с
вычислительной точки зрения случаях. Произведен поиск оптимальных численных схем в
плане чувствительности результатов расчетов к густоте используемой расчетной сетки и шагам
по времени. Разработанный программный код применен к решению тестовых и модельных
задач о нестационарном натекании потока воды препятствия различной формы, в частности, на
одиночное треугольное препятствие и на множественные (однорядные и двурядные)
препятствия в форме параллелепипедов.</p>
      <p>На примере задачи с треугольным препятствием показано, что пристенные вязкие эффекты
могут оказывать весьма сильное влияние на картину течения через формирование крупных
отрывных зон, приводящих к возникновению «холмов» на свободной поверхности. Установлено,
что вязкие эффекты вблизи боковых стенок канала делают картину течения существенно
трехмерной. Результаты трехмерного расчета хорошо согласуются с экспериментальными данными.</p>
      <p>Для натекания потока на множественные препятствия получен интересный для практики
вывод о том, что наличие первого ряда препятствий может привести к увеличению
опрокидывающего момента, действующего на препятствия второго ряда.
Литература
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.
11. Gaskell P.H., Lau A.K.C . Curvature-compensated convective-transport: SMART, A new
boundedness-preserving transport algorithm // Int. J. Numer. Methods Fluids. 1988.Vol. 8. Issue 6.</p>
      <p>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.
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
Computational Methods in Applied Sciences and Engineering (ECCOMAS 2012), Vienna, Austria,
10-14 September 2012. CD-ROM. ISBN: 978-3-9502481-9-7. 2 p.
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.</p>
      <p>CD-ROM, CMFF15-144.pdf. 8 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.</p>
      <p>Numerical simulation of 3D unsteady free-surface flows:
development and application of a specialized parallel code</p>
    </sec>
    <sec id="sec-2">
      <title>A.I. Khrabry, D.K. Zaytsev, E.M. Smirnov</title>
    </sec>
    <sec id="sec-3">
      <title>Saint-Petersburg Polytechnic University</title>
      <p>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.</p>
      <p>Keywords: numerical simulation, parallel code, domain decomposition, free-surface flow,
flow-obstacle interaction, turbulent flows, VOF method, viscous separation.
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.</p>
      <p>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.</p>
      <p>Hu C., Sueyoshi M. Numerical Simulation and Experiment on Dam Break Problem // Journal of
Marine Science and Application. 2010. Vol. 9. P. 109-114.
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.
9. Smirnov, E.M., Zaitsev D.K. Modification of Wall Boundary Conditions for Low-Re k-
Turbulence 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.
11. Gaskell P.H., Lau A.K.C . Curvature-compensated convective-transport: SMART, A new
boundedness-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.
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
Computational Methods in Applied Sciences and Engineering (ECCOMAS 2012), Vienna, Austria,
10-14 September 2012. CD-ROM. ISBN: 978-3-9502481-9-7. 2 p.</p>
    </sec>
  </body>
  <back>
    <ref-list>
      <ref id="ref1">
        <mixed-citation>
          4. Fedotova
          <string-name>
            <surname>Z.I.</surname>
          </string-name>
          , Chubarov L.B.
          <article-title>Chislennoe modelirovanie nakata tsunami [Numerical modeling</article-title>
          of tsunami rolling] // Trudy Mezhdunarodnoy konferentsii RDAMM-2001
          <source>[Proceedings of International Conference RDAMM-2001]</source>
          ,
          <year>2001</year>
          . Vol.
          <volume>6</volume>
          . Part 2.
          <string-name>
            <given-names>Special</given-names>
            <surname>Issue</surname>
          </string-name>
          . P.
          <volume>380</volume>
          -
          <fpage>396</fpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref2">
        <mixed-citation>
          <string-name>
            <given-names>Hirt C.W.</given-names>
            ,
            <surname>Nichols</surname>
          </string-name>
          <string-name>
            <surname>B.D.</surname>
          </string-name>
          <article-title>Volume of fluid (VOF) method for the dynamics of free boundaries //</article-title>
          <source>Journal of computational physics. 1981</source>
          . Vol.
          <volume>39</volume>
          . P.
          <volume>201</volume>
          -
          <fpage>226</fpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref3">
        <mixed-citation>
          <string-name>
            <given-names>Vandoormaal J.P.</given-names>
            ,
            <surname>Raithby</surname>
          </string-name>
          <string-name>
            <surname>G.D.</surname>
          </string-name>
          <article-title>Enhancements of the SIMPLE method for predicting incompressible</article-title>
          fluid flows // Numerical Heat Transfer.
          <year>1984</year>
          . Vol.
          <volume>7</volume>
          , issue 2. P.
          <volume>147</volume>
          -
          <fpage>163</fpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref4">
        <mixed-citation>
          <string-name>
            <surname>Menter</surname>
            ,
            <given-names>F.R.</given-names>
          </string-name>
          <article-title>Two equation eddy-viscosity turbulence models for engineering applications</article-title>
          // AIAA Journal.
          <year>1994</year>
          . Vol.
          <volume>32</volume>
          . P.
          <volume>1598</volume>
          -
          <fpage>1605</fpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref5">
        <mixed-citation>
          10.
          <string-name>
            <surname>Khrabry</surname>
            <given-names>A.I.</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Smirnov</surname>
            <given-names>E.M.</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Zaytsev</surname>
            <given-names>D.K.</given-names>
          </string-name>
          <article-title>Chislennoe modelirovanie techeniy so svobodnoy poverkhnost'yu na osnove metoda VOF [Numerical simulation of free surface flows on the base of the VOF method] // Trudy Krylovskogo gosudarstvennogo nauchnogo tsentra [</article-title>
          <source>Proceedings of Krylov state research centre]</source>
          .
          <year>2013</year>
          . Vol.
          <volume>78</volume>
          . P.
          <volume>53</volume>
          -
          <fpage>64</fpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref6">
        <mixed-citation>
          14.
          <string-name>
            <surname>Waclawczyk</surname>
            <given-names>T.</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Koronowicz</surname>
            <given-names>T.</given-names>
          </string-name>
          <article-title>Remarks on prediction of wave drag using VOF method with interface capturing approach</article-title>
          // Archives of Civil and Mechanical Engineering.
          <year>2008</year>
          . Vol.
          <volume>8</volume>
          . P. 5-
          <fpage>14</fpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref7">
        <mixed-citation>
          15.
          <string-name>
            <surname>Khrabry</surname>
            <given-names>A.I.</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Smirnov</surname>
            <given-names>E.M.</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Zaytsev</surname>
            <given-names>D.K.</given-names>
          </string-name>
          <article-title>Solving the Convective Transport Equation with Several High-Resolution Finite Volume Schemes: Test Computations</article-title>
          / In: Computational Fluid Dynamics 2010. Ed.:
          <string-name>
            <given-names>A.</given-names>
            <surname>Kuzmin</surname>
          </string-name>
          . New-York: Springer,
          <year>2011</year>
          . 954 p. P.
          <volume>535</volume>
          -
          <fpage>540</fpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref8">
        <mixed-citation>
          17.
          <string-name>
            <surname>Zaytsev</surname>
            <given-names>D.K.</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Smirnov</surname>
            <given-names>E.M.</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Khrabry</surname>
            <given-names>A.I.</given-names>
          </string-name>
          <article-title>Numerical simulation of free surface flows: influence of scheme factors and turbulence model //</article-title>
          <source>Proc. 14th International Conference “Supercomputing and mathematical modeling”</source>
          , Sarov, Russia.
          <year>2012</year>
          . P.
          <volume>282</volume>
          -
          <fpage>292</fpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref9">
        <mixed-citation>
          18.
          <string-name>
            <surname>Khrabry</surname>
            <given-names>A.</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Smirnov</surname>
            <given-names>E.</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Zaytsev</surname>
            <given-names>D.</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Goryachev</surname>
            <given-names>V</given-names>
          </string-name>
          .
          <article-title>Numerical study of separation phenomena in the dam-break flow interacting with a triangular obstacle //</article-title>
          <source>Proceedings of the 16th International Conference on Fluid Flow Technologies (CMFF'15)</source>
          , Budapest, Hungary, September 1-4
          <year>2015</year>
          .
          <article-title>CD-ROM</article-title>
          ,
          <fpage>CMFF15</fpage>
          -
          <lpage>144</lpage>
          .pdf. 8 p.
        </mixed-citation>
      </ref>
      <ref id="ref10">
        <mixed-citation>
          19. Golub G.,
          <string-name>
            <given-names>Van Loan С. Matrichnye</given-names>
            <surname>Vychisleniya</surname>
          </string-name>
          [Matrix Computations]. M.: Mir [Moscow, Mir publishing],
          <year>1999</year>
          . 553 p.
        </mixed-citation>
      </ref>
      <ref id="ref11">
        <mixed-citation>
          20.
          <string-name>
            <surname>Van Der Vorst H.A. Bi-CGSTAB</surname>
          </string-name>
          :
          <article-title>A fast and smoothly converging variant of Bi CG for the solution of nonsymmetric linear systems //</article-title>
          <source>SIAM Journal on Scientific and Statistical Computing</source>
          .
          <year>1992</year>
          . Vol.
          <volume>13</volume>
          . P.
          <volume>631</volume>
          -
          <fpage>644</fpage>
          .
        </mixed-citation>
      </ref>
    </ref-list>
  </back>
</article>