CC BY
30
УДК 532.5
А. С. Козелков1'2, Д. Ю. Стрелец3, В. Р. Ефремов4, Ю. Г. Нечепурепко4, В. В. Курулин1, Е. С. Тятюшкина1, А. В. Корнев3
1 Российский федеральный ядерный центр - Всероссийский научно-исследовательский институт экспериментальной физики 2Нижегородский государственный технический университет им. Р. Е. Алексеева Московский авиационный институт 4Конструкторское бюро приборостроения им. академика А. Г. Шипунова
Исследование свойств схем дискретизации уравнения переноса объемной доли при расчете многофазных течений методом VOF
В работе представлены результаты исследования свойств схем дискретизации уравнения переноса объемной доли при расчете многофазных течений методом VOF (Volume of Fluid) на произвольных неструктурированных сетках. Представлено описание способов сохранения формы поверхности и проанализированы свойства наиболее известных специальных схем дискретизации уравнения переноса объемной доли, а также влияние числа Куранта на перенос фронта. В результате исследования определено предельное значение числа Куранта для различных схем и типов расчетных сеток, необходимое для сохранения формы скаляра после его переноса.
Ключевые слова: уравнение Навье-Стокса, численное моделирование, многофазные течения, свободная поверхность, неструктурированные сетки, перенос фронта, пакет программ ЛОГОС.
A.S. Kozelkov1'2, D. Yu. Strelets3, V. R. Efremov4, Yu. G. Nechepurenko4, V. V. Kurulin1, E. S. Tyatyushkina1, A. V. Kornev3
1 FSUE «RFNC-VNIIEF» 2 Nizhny Novgorod State Technical University n.a. R.E. Alekseev 3 Moscow Aviation Institute (National Research University) 4 JSC «KBP named after Academician A. Shipunov»
Investigation of the properties of discretization schemes for the volume fraction transfer equation in the calculation of multiphase flows by the VOF method
The paper summarizes the results of studying the properties of discretization schemes for the volume fraction transport equation in VOF (Volume of Fluid) multiphase flow simulations on arbitrary unstructured meshes. Surface shape preservation techniques are presented, the properties of the most popular ad hoc discretization schemes for the volume fraction transport equation are analyzed, and the effect of the Courant number on front transport is examined. The obtained results provide the limiting value of the Courant number using different schemes and types of a computational mesh required to preserve the shape of the scalar field after its transport.
Key words: Navier-Stokes equations, numerical vodeling, multiphase flows, free surface, unstructured meshes, front transport, LOGOS.
© Козелков А. С., Стрелец Д. Ю., Ефремов В. Р., Нечепуренко Ю. Г., Курулин В. В.,
Тятюшкина Е. С., Корнев А. В., 2017 © Федеральное государственное автономное образовательное учреждение высшего образования «Московский физико-технический институт (государственный университет)», 2017
1. Введение
В настоящее время существует несколько методов моделирования многофазных течений со свободной поверхностью, которые отличаются способом ее расчета. Первый метод основан на «лагранжевом» подходе, в котором свободная поверхность отслеживается либо путем перемещения узлов сетки, либо частицами [1,2]. Второй метод основан на «эйлеровом» подходе, в котором для отслеживания свободной поверхности вводятся специальные маркеры, в роли которых могут выступать либо частицы [3], либо пространственные маркер-функции, подчиняющиеся уравнению конвективного переноса [4-6]. Наиболее применимым на практике является второй метод с использованием в качестве маркер-функции объемной доли жидкости (VOF - Volume-of-fluid) [6]. В данном методе система жидкость-газ рассматривается как единая односкоростная среда с переменными физическими свойствами. Метод легко обобщается на случай произвольных неструктурированных сеток и произвольного количества фаз [7,8]. Применение этого метода требует решения уравнения переноса объемных долей. Выбор корректной численной схемы для дискретизации данного уравнения напрямую влияет на итоговую форму свободной поверхности. Основное требование, предъявляемое к используемой схеме - это сохранение «резкости» фронта раздела между жидкостью и газом, а также сохранение формы свободной поверхности при ее параллельном переносе и вращении. В настоящей работе проанализированы свойства наиболее известных специальных схем дискретизации уравнения переноса объемной доли, а также влияние числа Куранта на перенос фронта. В результате исследования определено предельное значение числа Куранта для различных схем и типов расчетных сеток, необходимое для сохранения формы скаляра после его переноса.
2. Основные уравнения и численные схемы
Будем предполагать, что в многофазной системе все фазы несжимаемы. Односкоростная модель, в которой движение всех фаз системы имеет общую скорость, а все параметры являются среднеобъемными характеристиками многофазной смеси, опирается на систему уравнений Навье-Стокса, дополненную уравнением переноса объемной доли [7,8]:
V- и = 0,
< ^ + V ■ (puu) = —Vp + V ■ (тш + Tt) + рд + F, (i)
a(k)p(k) + V ■ a(k)p(k)u(k) = 0,
V
где u = 2=1 a(fc)p(fc)-u(fc))/p - вектор общей скорости многофазной смеси по совокупности фаз k\ п - количество фаз; р = ^a(fc)p(fc) - плотность многофазной смеси; тш и rt - молекулярная и турбулентная составляющие тензора касательных напряжений многофазной смеси; F - силы поверхностного натяжения. Дискретизация системы (1) может быть осуществлена любым известным способом [9-11]. Если при разработке численного метода ориентироваться на решение задач на произвольных неструктурированных сетках с ячейками произвольной формы, оптимальным выбором дискретизации является метод конечных объемов. Численное решение первых двух уравнений системы (1) для моделирования течений с несжимаемыми фазами целесообразно осуществлять с помощью известного алгоритма SIMPLE [9-11] в полностью неявной модификации, снимающей жесткие ограничения на шаг по времени [8]. Для моделирования границ раздела фаз после решения первых двух уравнений системы (1) решается уравнение переноса объемных долей (третье уравнение системы (1)) для (п — 1) объемной доли фаз. Его дискретизация методом конечных объемов осуществляется по схеме, полностью аналогичной той, что используется для уравнения сохранения импульса (второе уравнение системы (1)) в методе конечных объемов. Дискретизация нестационарного слагаемого осуществляется по схемам Эйлера
или Адамеа Бэшфорта [8]. Отличительной численной особенностью этого уравнения является аппроксимация конвективного слагаемого, которая должна обеспечивать сохранение минимально возможной толщины границы раздела сред, т.е. не допустить размытие контактной границы, и сохранение ограниченности значения объемной доли в интервале от 0 до 1. Для аппроксимации конвективного слагаемого уравнения переноса объемных долей используют так называемые «сжимающие» схемы [7,8]. Принцип построения сжимающих схем следующий. Рассмотрим внутреннюю грань являющуюся смежной для ячеек Р и N (рис. 1).
Рис. 1. Схема дискретизации конвективного слагаемого (п нормаль к грани)
Введем нормализованную переменную [7,12]:
* У -Уи
У =
. (2)
Ум - Уи
где Уи ~ значение в ячейке выше по потоку ячейки Р.
Значения нормализованной переменной в ячейке Р на храни { равны:
* Ур- Уи * У -Уи
Ур =-, У* =-. (3)
УМ - Уи Ум - Уи
В случае неструктурированной сетки выражение для нормализованной переменной
можно представить в следующем виде [13]:
У = 1 -
Ум - У 2( Й )рЛг
(4)
здесь йг - расстояние между центрами ячеек.
У
У = (1 - 7г) •УР + Ц • УМ,
(о)
где 7/ =(У* У*р)/(1 -У*Р).
Согласно критерию локальной охраниченности СВС [12] одномерное решение конвективного переноса будет ограниченным, если значение нормализованной переменной на гра-
У
У = У*м ^и У*м < 0 или У^ > 1, Ум <У*{ < 1, при 0 <У*М < 1. Для явных схем критерий СВС является более жестким и имеет вид
У = У*м, при У*м < 0 или У*м > 1, Ум < У < шш(1, ^^), при 0 < Ум < 1,
(6)
где CFL - число Куранта [10].
Одна из самых известных сжимающих схем, схема HRIC, работает в несколько этапов [14]. На первом этапе нормализованное значение величины р на грани f рассчитывается с использованием схемы, соответствующей верхней границе критерия СВС при CFL = 0.5:
^ при < 0 или > 1, <Pt* = <| при 0 < ¡p*N < 0.5, (8)
1 при 0.5 < ip*N < 1.
На следующем этапе проводится корректировка схемы с учетом расположения межфазной границы - такой этап необходим для подавления искусственного искривления межфазной границы. Получим значение
(^*)' = А/ ■ ^* + (1 — Xf) ■ , (9)
где \f = \J| cos Of|. Угол в определяется между вектором, соединяющим центры соседних ячеек N и Р, и вектором нормали к межфазной границе.
На последнем этапе производится учет числа Куранта. Схема корректируется для приведения ее к критерию СВС:
'(<£**)' при CFL < 0.3, р* = <| <Р%пР* CFL > 0.7, (10)
<P*N + ((&*)' — <Pn)Ж—Шг ПРИ 0.3 < CFL < 0.7.
Данная схема в меньшей степени зависит от числа Куранта.
Схема HYPER-C обеспечивает ограниченность решения таким образом, чтобы схема аппроксимации величины на грань между ячейками удовлетворяла верхней границе критерия СВС:
<p*N при <p*N < 0 или <p*N > 1, min(1 ) при 0 < <p*N < 1.
)HYPER—С = { 1 ____n ^ 1 (П)
Схема М-СЮБАМ описана в работе [15]. Данная схема обладает лучшими качествами схем НШС и СК'ЯЛ.М. в формулировке которой не содержится числа Куранта. За основу была взята схема НШС, соответствующая верхней границе критерия СВС при СРЬ = 0.5:
(<£*)м-С1СЯАМ = А ■ <Р** + (1 - А) ■ (<£*)ркомм, (12)
1
где fo = (| cos 9f\) 4
ip*N при < 0 или > 1,
)рКОММ = < . ,1 , * / ч П/ * (13)
1 1 шт(4 + , (<р**)н¥РЕК-с) при о < < 1.
Таким образом, схема М-СЮБАМ обеспечивает ограниченность решения, удовлетворяет критерию локальной ограниченности (критерию СВС).
В алгебраической форме система уравнений для определения характеристик к-й фазы имеет вид
(к)Лк) , ^ „(*)„,(*)_ к(к)
aw a
Р + £ a^ ог£ = bf. (14)
N=NB ( Р)
Коэффициенты матрицы неявного решения уравнения (14) имеют вид [8]
ар = ам ' т> (15)
а$ = ш1п(0, и(а 113{),
¿к) _ _ (к) + V
ам + г,
N=МВ(Р)
Ак) (п- 1)су (к) (к)Л 1 V (к),п-1
Ьр = -щ 3{(аумс -а>и'в) + ~ар , где и - значения объемной доли на грани, найденные по схеме МСЮБАМ и по
(к),п— 1
противопоточноп схеме соответственно; ар - значение объемной доли на предыдущий шаг по времени. Представленные слагаемые получены с помощью дискретизации методом отложенной коррекции [7]. Для численного решения итоговая система уравнений должна быть дополнена начальными и граничными условиями. На твердых стенках градиент давления и объемных долей равен нулю: ^ = 0 = 0; значение скорости равно нулю: и = 0,г/ = 0,ш = 0. На входе задано фиксированное значение скорости и объемных долей, градиент давления равен нулю: ^ = 0. На выходе фиксируется статическое давление; градиенты скорости и объемных долей равны нулю: ^ = 0 ^ = 0 ^¡^ = 0 = 03. Технология расчета
Существует несколько способов сохранения резкости фронта:
1) применение «сжимающих» численных схем, описанных выше;
2) применение дополнительных слагаемых в уравнении переноса объемной доли, которые работают как антидиффузия в зоне фронта.
Первый способ имеет свои преимущества: он не изменяет исходного уравнения, а значения объемной доли на грани вычисляются с помощью специальных схем (НШС, МСГСБАМ). Однако применение таких схем на практике приводит к существенному ограничению шага по времени. Так, в [17] показано, что приемлемые результаты при использовании схем НШС и МСГСБАМ получаются при значении числа Куранта СЕЬ < 0.75. Такие оценки ограничения шага по времени получены на структурированной расчетной сетке при заданной фиксированной скорости. Оценки делаются только на основании решения уравнения переноса, т.е. уравнения движения и неразрывности не решаются. Вероятно, использование неструктурированной расчетной сетки потребует еще больших ограничений шага по времени.
Второй способ известен из работ [17,18]. В этом способе для предотвращения увеличения ширины фронта предложено использовать дополнительное слагаемое, которое приведет третье уравнение системы (1) к виду (3) (номер фазы к опущен):
да
— + V- (аи)+ V- (а(1 - а)у) = 0. (16)
д
а
Данное слагаемое отлично от нуля в переходной области, где объемная доля изменяется в 0< а<1
полем скорости у, которое должно вычисляться таким образом, чтобы объемная доля пе-
а
скорость у на грани { таким образом:
••-^'-[У? (17)
где - объемный поток, С - коэффициент сжатия. Такой механизм позволяет значительно уменьшить ширину фронта, однако при больших значениях коэффициента сжатия возможно появление численных осцилляций решения.
Также существует и третий способ численной дискретизации уравнения переноса объемной доли, в котором применяются оба представленных подхода одновременно: для дискретизации оригинального уравнения используются сжимающие схемы, а для дополнительного сжатия фронта используют модифицированное уравнение (3).
В данной работе исследуется вопрос об оптимальном варианте применения описанных выше подходов для конечно-объемной реализации с применением полностью неявной дискретизации на неструктурированной расчетной сетке. Эффективность каждого из подходов оценивается по максимальному шаху по времени, при котором подход дает приемлемый результат на сетках различной структуры.
4. Численные эксперименты
Представленные ниже численные эксперименты проводились с использованием пакета программ ЛОГОС. Пакет программ ЛОГОС предназначен для решения сопряженных трехмерных задач конвективного тепломассопереноса, аэродинамики и гидродинамики на параллельных ЭВМ [19 21]. Пакет программ ЛОГОС успешно прошел верификацию и показал достаточно хорошие результаты на серии различных гидродинамических задач, включая расчеты турбулентных и нестационарных течений [22 24]. В настоящее время пакет программ ЛОГОС уже успешно используется для решения фундаментальных задач геофизики [25 27], а также для решения индустриальных задач [30,31]. Ускорение расчетов с помощью высокопараллельных вычислительных систем осуществляется многосеточными технологиями, детальное описание которых можно найти в [32,33].
4.1. Задача о переносе скаляра
Будем рассматривать классическую задачу о переносе скаляра постоянным полем скорости [7]. Данная задача наиболее удобна для оценки свойств представленных подходов, поскольку не связана с решением уравнений Навье Стокса, а отдельно рассматривает решение уравнения переноса объемной доли для скаляра.
Задача состоит в следующем. В начальный момент поле скалярной величины ограничено двумя квадратными рамками (рис. 2). Размер стороны внешней рамки 0.8 м., внутренней 0.4 м. В начальном положении центр рамки имеет координаты (0.8; 0.8) м. Расчетная область представляет собой квадрат со стороной 4 м.
а = 0
а = 1
Рис. 2. Схематичная постановка задачи
Были построены три расчетные сетки. Первая из них это равномерная сетка, которая состоит из правильных шестигранников, общее количество элементов составляет 40 тыс. ячеек (рис. За). Вторая сетка состоит из призм, общее количество элементов составляет 89 тыс. ячеек (рис. 36). Третья сетка состоит из произвольных многогранников с общим количеством элементов 33 тыс. ячеек (рис. Зв). Характерный размер ячейки для всех сеток одинаков и составляет 0.02 м.
а) б) в)
Рис. 3. Фрагмент расчетной сетки в области первоначального положения скалярной величины: (а структурированная шестигранная сетка, б неструктурированная сетка из треугольников, в полиэдральная сетка)
Скалярная величина перемещается фиксированным нолем скорости, направленной иод углом к линиям расчетной сетки, равной 8 м/с и 4 м/с но осям соответственно (рис. 2). Скаляр переносится на расстояние, равное 2.5 м, что составляет около трех размеров рамки и является достаточным, чтобы провести сравнительный анализ конечной формы скаляра. При таком плоекопараллельном перемещении конечная форма скаляра, заключенная в данные рамки, должна оставаться неизменной. Однако присущая всем численным схемам диффузия будет вносить искажение в первоначальную форму. По величине этого искажения можно судить о погрешности рассматриваемых подходов для каждой из применяемых схем.
Конечная форма скаляра зависит от выбранной дифференциальной схемы аппроксимации конвективного уравнения переноса, а также от значения числа Куранта. Задача первого этана исследования заключается в определении влияния значения числа Куранта и тина расчетной сетки на качество результатов: в конечном счете, будет определено предельное значение числа Куранта для различных сеток, необходимое для сохранения формы скаляра после его переноса.
Для аппроксимации конвективного уравнения переноса берутся следующие схемы: стандартная противопоточная схема I Г) [9,22] и схемы, относящиеся к классу «сжимающих», I II» ГС и МСГСБАМ. Аппроксимация нестационарного слагаемого осуществляется но неявной схеме Адамса Бэшфорта. В таблице 1 приведены параметры численных экспериментов указан номер эксперимента, схема и шаг но времени, соответствующий определенному числу Куранта.
Т а б л и ц а 1
Параметры численных экспериментов для структурированной сетки
Дифференциальная схема , 10-3 с. Число Куранта, СП
НИ 1С. МП СНАМ 0.418 0.25
НИ 1С. МП СНАМ 0.835 0.5
НИ 1С. МП СНАМ 1.25 0.75
НИ 1С. МП СНАМ 1.67 1
НИ 1С. МП СНАМ 2.084 1.25
НИ 1С. МП СНАМ 2.5 1.5
Результаты расчетов приведены на рис. 4. Видно, что с ростом числа Куранта качество решения, получаемого с помощью всех схем, ухудшается: сильным изменениям подвергается форма квадратной рамки, правый верхний и левый нижний углы вытягиваются но направлению потока, а правый нижний угол становится более округлым. Также но-
являются искажения решения в виде ненулевых объемных долей, выходящих за пределы внешней квадратной рамки. Такие изменения наблюдаются для всех случаев при значениях CFL > 1. Это свидетельствует о неустойчивости численного решения для данных чисел Куранта. При значениях CFL < 1 конечная форма скаляра несколько искажается: правые верхний и нижний углы подвергаются искажению, верхний угол слегка вытягивается, а нижний округляется. При значениях CFL 0.25 и CFL 0.5 конечная форма скаляра остается практи чески неизменной.
Можно отметить, что схемы HRIC и MCICSAM практически одинаково искажают конечную форму после переноса скаляра при всех значениях числа Куранта, при этом схема MCICSAM сильнее размывает конечную форму скаляра. Рассмотрим значения числа Куранта, при которых конечная форма скаляра подверглась искажению (рис. 4). Коэффициент сжатия будем варьировать в диапазоне 0.1 < С < 1с шагом 0.1. Так как искажение конечной формы скаляра при увеличении числа Куранта происходит практически одинаково для обеих схем HRIC и MCICSAM, будем представлять следующие результаты только для первой из них.
GFL 0.25, HRIC
GFL 0.75, HRIC
GFL 1.25, HRIC
GFL 0.25, MCICSAM
GFL 0.75, MCICSAM
GFL 1.25, MCICSAM
Рис. 4. Коночная форма скаляра поело переноса с использованием схем НШС и МСТСБАМ
Видно, что искажение формы скаляра и присутствие ненулевой объемной доли вне квадратной рамки растут с увеличением значения числа Куранта, поэтому рассматривать значения числа Куранта выше 1.5 уже не имеет смысла, что в итоге может привести к развалу решения. Предельным значениям на структурированной сетке для схем РШ 1С и МСГСБАМ можно считать значение СЕЬ 0.75.
Однако появление численных неустойчивостей при значении GFL > 1 можно предотвратить путем применения коэффициента сжатия, т.е. третьих) способа для решения уравнения конвективного переноса. Использование такого подхода, возможно, позволит повысить предельное значение числа Куранта. На рисунках ниже приведены результаты расчетов (рис. 5).
С 0.1, GFL 1
С 0.5, GFL 1
С 0.1, GFL 1.25
С 0.5, GFL 1.25
С 0.1, CFL 1.5
С 0.5, CFL 1.5
Рис. 5. Конечная форма скаляра после переноса с использованием схемы HRIC при различных CFL
Как видно из рис. 5, при С 0.1 для CFL 1,25 и CFL 1,5 все еще присутствуют ненулевые объемные доли, выходящие за пределы квадратной рамки. Однако с увеличением коэффициента сжатия, начиная со значения С 0.5, данные неустойчивости исчезают, и одновременно конечная форма скаляра начинает симметрично изменяться, вытягивая левый верхний и правый нижний углы поперек направлению потока, а также округляя противоположные углы. При значении С 0.5 форма квадратной рамки вытягивается
симметрично относительно диагонали квадрата и при этом достаточно сильно срезает противоположные углы. При CFL 1.5 вплоть до С 0.5 в решении присутствуют численные осцилляции, которые в дальнейшем исчезают.
Таким образом, предельными значениями числа Куранта и коэффициента сжатия в данной серии расчетов (показанной на рис. 5) является случай при CFL 1.5 и С 0.5. При этих значениях еще возможно проведение численного эксперимента с удовлетворительным качеством. Их повышение уже не гарантирует решения задачи.
Универсальность выбранного критерия проверим на неструктурированных сетках. Для этого проведем серию численных экспериментов при значениях CFL 1.5 и С 0.5. Расчеты на треугольной сетке проводятся с шагом но времени 0.00105 на полиэдральной сетке с шагом 0.00133, которым на структурированной сетке соответствует значение числа Куранта 1.5.
На рис. 6 представлены результаты расчетов. Видно, что возникающие осцилляции решения при таком значении числа Куранта, скорее всего, приведут к развалу решения. Выбранное предельное значение для числа Куранта и коэффициента сжатия на структурированных сетках не подходит для неструктурированных сеток. Для определения общего критерия для неструктурированных сеток понизим значения числа Куранта до 0.75 и проведем серию расчетов, параметры для которых приведены в таблицах 2 и 3 соответственно.
Треугольная сетка Полиэдральная сетка
Рис. 6. Коночная форма скаляра поело переноса на треугольной и полиэдральной сетках при СТО 1.5 и С 0.5
Т а б л и ц а 2
Параметры численных экспериментов для треугольной сетки
Ai,10-3 с Число Куранта, CFL
0.525 0.75
0.7 1
0.871 1.25
Т а б л и ц а 3
Параметры численных экспериментов для полиэдральной сетки
At,10-3 с Число Куранта, CFL
0.664 0.75
0.89 1
1.33 1.25
На рис. 7 видно, что при СРЬ 0.75 размытие границы скаляра на треугольной сетке практически отсутствует. С увеличением числа Куранта, как и следовало ожидать, растут ненулевые включения вне квадратной области скаляра, что, так же, как и в предыдущих случаях, может привести к развалу решения при использовании дополнительного уточнения решения, т.е. третьих) подхода.
СИ 0.75 СП 1 СП 1.25
Рис. 7. Коночная форма скаляра поело переноса на треугольной сетке для нескольких значений СРЬ
На рис. 8 мы видим результаты расчетов на полиэдральной сетке. Отметим, что на данной сетке при любом значении Куранта мы получили приемлемые результаты. По результатам представленных расчетов можно сделать вывод, что применяя только сжимающую схему без каких-либо дополнительных соотношений, предельным значением на неструктурированных сетках является значение числа Куранта 0.75. Однако, если рассматривать только полиэдральную сетку, для нее предельным в данном исследовании является значение числа Куранта, равное 1.5.
СП 0.75
СИ 1
СИ 1.25
Рис. 8. Конечная форма скаляра после переноса на полиэдральной сетке для нескольких значений СРЬ
Остается вопрос, можно ли повысить значение числа Куранта для треугольной сетки, одновременно сохраняя конечную форму скаляра после переноса. Попробуем уточнить решение, получаемое на треугольной сетке. Для этох'о проведем серию расчетов с использованием третьих) подхода решения конвективного уравнения переноса на треухххльной сетке для значений числа Куранта 0.75, 1, 1.25, а коэффициент сжатия будем изменять в диапазоне предельнохх) значения для структурированных сеток до 0.4 с шагом 0.1.
На рис. 9 видно, что применение коэффициента сжатия привело к некоторому улучшению конечнохх) результата, и нам удалось повысить предельное значение числа Куранта с 0.75 до 1 с применением сжатия при С 0.1. Но очевидно, что с увеличением коэффициента сжатия конечная форма скаляра искажается все сильнее. Для СРЬ 1.25 ни с одним значением коэффициента сжатия улучшить результат не удалось. Таким образом, можно сказать, что на треугольных сетках мы выдерживаем качественное решение при СРЬ 1 и С 0.1.
Так как мы определяем общий критерий для структурированных сеток, проведем аналогичные расчеты на полиэдральной сетке. Результаты приведены на рис. 10. Видно, что с увеличением коэффициента сжатия происходит развал решения, и предельным, как и
для треугольных сеток, является коэффициент сжатия С 0.1. Однако при СРЬ 1.25 и С 0.1 мы также получаем приемлемые результаты в отличие от треугольной сетки, где не удалось убрать численные осцилляции путем применения дополнительного сжатия.
С 0.1, С 0.2, С 0.1, С 0.2, С 0.1, С 0.2, СРЬ 0.75 СРЬ 0.75 СРЬ 1 СРЬ 1 СРЬ 1.25 СП 1.25
Рис. 9. Коночная форма скаляра на треугольной сотко при различных СРЬ
С 0.1, С 0.2, С 0.1, С 0.2, С 0.1, С 0.2,
СРЬ 0.75 СРЬ 0.75 СРЬ 1 СРЬ 1 СРЬ 1.25 СРЬ 1.25
Рис. 10. Коночная форма скаляра на полиэдральной сотко при СРЬ 0.75
По результатам исследования можно сказать, что для неструктурированных сеток общим предельным значением числа Куранта является СРЬ 0.75 (рис. 10). При использовании дополнительного сжатия на величину С 0.1, т.е. третьего подхода для решения конвективного уравнения переноса, удалось повысить значение числа Куранта до 1 (рис. 11, 12). Предельным значением для структурированных сеток является значение СРЬ 0.75 (рис. 13).
Треугольная сетка, НШС
Полиэдральная сетка, НШС
Треугольная сетка, МПСЯЛМ
Полиэдральная сетка, МПСЯЛМ
Рис. 11. Коночная форма скаляра поело переноса на неструктурированных сетках с использованием первого подхода для решения конвективного уравнения переноса
Однако так же, как и для неструктурированных сеток, применяя дополнительное сжатие на величину С 0.4, удалось убрать присутствующие в решении неустойчивости и повысить предельное значение числа Куранта до 1.5.
Описанные выше результаты были получены с использованием сжимающих схем. Проведем расчеты на структурированной сетке с использованием стандартной схемы I, Г) (рис. 14). Для улучшения результатов применим сжатие.
Видно, что стандартные схемы вообще являются неприемлемыми для такого типа задач. Однако использование дополнительного сжатия (второй подход) при моделировании может уточнить форму переносимой величины.
П
Треугольная сетка, НШС
Полиэдральная сетка, НШС
Треугольная сетка, МСК'НАМ
Полиэдральная сетка, МСК'НАМ
Рис. 12. Коночная форма скаляра поело переноса на неструктурированных сетках с использованием третьего подхода для решения конвективного уравнения переноса
Первый подход,
ншс
Первый подход, МСГСБАМ
Третий подход, НШС
Третий подход, МСГСБАМ
Рис. 13. Конечная форма скаляра после переноса на структурированных сетках для решения конвективного уравнения переноса
СО С 0.1 С 0.5 С 1
Рис. 14. Конечная форма скаляра после с применением схемы 1ГО при числе Куранта СРЬ 1
5. Практическое применение
Мно1Х)фазные течения достаточно часто встречаются в различных индустриальных и научных приложениях [7,8,25 29]. Такие течения встречаются практически везде, где имеется жидкость, и играют огромную роль в природе и технике: поверхностные волны, капли, струи, пузыри, пленки и т.д. Экспериментальное исследование подобжих) рода задач зачастую сопряжено со значительными трудностями и затратами, поэтому разработка эффективных численных методов и алх'оритмов, способных достоверно описывать такие течения является актуальной задачей. Среди мши'офазных течений особую сложность представляют течения со свободной поверхностью, частным случаем которых являются течения неемешивающихея жидкостей. Основной проблемой при моделировании таких течений является отслеживание поверхности раздела на фиксированной расчетной сетке.
Проведенное исследование позволяет эффективно использовать, описываемые схемы и подходы для решения практических задач. На рис. 15 представлены результаты численно-IX) моделирования задачи об обрушении плотины [34]. Расчет, в отличие от [7], выполнен на произвольной неструктурированной сетке, и из результатов видно, что оптимальный под-
бор параметров, основанный на проведенном исследовании позволяет весьма точно описать эксперимент как качественно, так и количественно (рис. 16, 17).
а) б)
Рис. 15. Обрушение плотины: а) результаты численного моделирования: б) экспериментальные данные [7]
■как-рнмен I я.1М1ьк-
1ЯННЫС рС)у.|ЫЯГЫ Ч11С.1Г111101 о
моделирования:
— ЛОГОС
--Н рябо1Ср]
О
0 12 3 4
Рис. 16. Зависимость высоты столба жидкости на левой стенке резервуара от времени
Сплошной черной линией изображена граница раздела фаз вода и воздух, серым цветом выделена жидкая фаза, черными штрихами показаны векторы скорости. Результаты построены для различных моментов времени, Т теоретический период колебаний но закону косинуса, который в нашем случае составляет 0.3739 с [7,34]. В начальный момент времени граница раздела фаз имеет форму косинуса. Затем ее отпускают и она начинает колебаться под действием силы тяжести и инерции жидкости. В момент времени £ = 1/4 -Т граница раздела фаз принимает вид горизонатальной линии, а затем в Ь = 1/2 ■ Т - снова косинусоиды. Процесс повторяется во времени. На рис. 19 показано колебание уровня жидкости вдоль левой стенки резервуара с течением времени. Погрешность по времени вычисляется следующим образом:
—■100%, и
где £з - период колебаний, полученный по результатам численного моделирования, £^ теоретический период Т [34].
Рис. 17. Зависимость положения породной кромки жидкости от времени
открытая граница
воздух 1
вода
Рис. 18. Задача о колебании жидкости под действием силы тяжести (1- схема экспериментальной установки, 2 - поверхность воды на период £ = 1/4 •Т, 3 - поверхность воды на период £ = 1/2 • Т)
о 0,060 <и %
а
0,040
0,0 0,5 1,0 1,5 2,0 2,5
Время, с
Рис. 19. Колебания уровня жидкости вдоль левой стенки резервуара под действием силы тяжести
Наиболее неустойчивыми с точки зрения моделирования являются задачи на равновесие. Ниже представлены результаты моделирования задачи о колебаниях воды в резервуаре под действием силы тяжести (подробное описание задачи дано в [34]).
Видно, что численные результаты дают хорошее качественное и количественное согласование с теорией и численным решением, приведенным в работе [7].
Работа выполнена при финансовой поддержке АО «КБП» им. Шипунова А.Г. (г. Тула) в рамках НИР «Моделирование гидротормозных устройств обращенного веретенного типа с наличием свободного объема в цилиндре с рабочей жидкостью», гранта Президента Российской Федерации для государственной поддержки ведущих научных школ Российской Федерации НШ-6637.2016.5, а также при финансовой поддержке РФФИ - проект № 16-0100267.
Литература
1. Батлер Т.М. Развитие метода LINC. Численные методы в механике жидкостей. М.: Мир, 1973. С. 146-155.
2. Бура,го Н.Г. Численное решение задач МСС с подвижными границами раздела: дис. на соискание уч. ст. д-ра физ.-мат. наук // Н.Г. Бураго. М.: РАН. 1111 \ 1. 2003. 222 с.
3. Lucy L.B. A numerical approach to the testing of the fission hypothesis // Astron. J. 1977. V. 82 (12). P 1013-1024.
4. Harlow F.H., Welch J.E. Numerical calculation of time-dependent viscous incompressible flow // Phvs. Fluids. 1965. V. 8. P. 2182-2189.
5. Daly B.J. A technique for including surface tension effects in hydrodvnamic calculations // J. Comput. Phvs. 1969. V. 4. P. 97-117.
6. Hirt C. W., Nichols B.D. Volume of fluid (VOF) method for the dynamics of free boundaries // Journal of computational physics. 1981. V. 39. P. 201-226.
7. Ubbink O. Numerical prediction of two fluid systems with sharp interfaces // PhD thesis, Imperial College. University of London. 1997.
8. Козелков А. С., Мелешкина Д.П., Куркип А.А., Тара,сова H.B., Лашкип С.В., Курулин В.В. Полностью неявный метод решения уравнений Навье-Стокса для расчета многофазных течений со свободной поверхностью // Вычислительные технологии. 2016. Т. 21, № 5. С. 54-76.
9. Ferziger J.H., Peric М. Computational Method for Fluid Dynamics. Springer-Verlag. New York, 2002.
10. Флетчер К. Вычислительные методы в динамике жидкостей в двух томах. М.: Мир, 1991.
11. Волков К.Н., Дерюгин Ю.Н., Козелков А. С., Емельянов В.Н., Тет,ерина, И. В. Разностные схемы в задачах газовой динамики на неструктурированных сетках. М.: Физматлит, 2014. С. 416.
12. Leonard В.P. The ULTIMATE conservative difference scheme applied to unsteady onedimensional advection // Comput. Methods Appl. Mech. Engrg. 1991. V. 88. P. 1774.
13. Jasak H., Weller H.G., Gasman A.D. High resolution NVD differencing scheme for arbitrarily unstructured meshes // International journal for numerical methods in fluids. 1999. V. 31. P. 431-449.
14. Muzaferija S., Peric M., Sames P., Schelin T. A twofluid Navier-Stokes solver to simulate water entry //22 Svmp, on Naval Hydrodynamics, 1998.
15. Waclawczyk Т., Koronowicz T. Remarks on prediction of wave drag using VOF method with interface capturing approach // Archives of civil and mechanical engineering. 2008. V. 8. P. 5-14.
16. Jasak H. Error Analysis and Estimation for the finite volume method with applications to fluid flows. Thesis submitted for the degree of doctor // Department of Mechanical Engineering, Imperial College of Science, 1996.
17. Khrabry A.I., Smirnov Е.М., Zaytsev D.K. Solving the Convective Transport Equation with Several High-Resolution Finite Volume Schemes: Test Computations // Computational Fluid Dynamics. 2010. P. 535-540.
18. Rusche H. Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions // PhD thesis, Imperial College, University of London, 2002.
19. Betelin V.B., Shagaliev R.M., Aksenov S. V., Belyakov I.M., Deryuguin Yu.N., Kozelkov A.S., Korchazhkin D.A., Nikitin V.F., Sarazov A.V., Zelenskiy D.K. Mathematical simulation of hydrogen-oxvgen combustion in rocket engines using LOGOS code // Acta Astronautica. 2014. V. 96. P. 53-64.
20. Бойко A.B., Нечепуренко Ю.М., Жучков P.H., Козелков А. С. Блок расчета положения ламинарно-турбулентного перехода для пакета программ ЛОГОС // Теплофизика и аэромеханика. 2014. Т. 21, № 2. С. 201-220.
21. Сафронов A.B., Дерюгин Ю.Н., Жучков Р.Н., Зеленский Д.К., Сара,зов A.B., Козелков A.C., Кудимов Н.Ф., Липницкий Ю.М., Панасенко A.B. Результаты валидации многофункционального пакета программ ЛОГОС при решении задач аэрогазодинамики старта и полета ракет-носителей // Математическое моделирование. 2014. Т. 26,№ 9. С. 83-95.
22. Kozelkov A., Kurulin V., Emelyanov V., Tyatyushkina Е., Volkov К. Comparison of convective flux discretization schemes in detached-eddv simulation of turbulent flows on unstructured meshes // Journal of Scientific Computing. 2015. N 89. P. 1-16.
23. Козелков A.C. Курулин B.B. Численная схема для моделирования турбулентных течений несжимаемой жидкости с использованием вихреразрешающих подходов // Вычислительная математика и математическая физика. 2015. Т. 55, № 7. С. 135-146.
24. Козелков A.C., Куркин A.A., Крутякова О.Л., Курулин В.В., Тятюшкина Е.С. Зонный RANS-LES подход на основе алгебраической модели рейнольдсовых напряжений // Известия РАН Механика жидкости и газа. 2015. Т. 50, № 5. С. 24-33.
25. Козелков A.C., Куркин A.A., Пелиновский E.H., Курулин В.В. Моделирование цунами космогенного происхождения в рамках уравнений Навье-Стокса с источниками различных типов // Известия РАН Механика жидкости и газа. 2015. V. 50, № 2. С. 306-313.
26. Козелков A.C., Куркин A.A., Пелиновский E.H. Влияние угла входа тела в воду на высоты генерируемых волн // Известия РАН. Механика жидкости и газа. 2016. Т. 51, № 2. С. 288-298.
27. Козелков A.C., Куркин A.A., Пелиновский E.H., Курулин В.В., Тятюшкина Е.С. Моделирование возмущений в озере Чебаркуль при падении метеорита в 2013 году // Известия РАН Механика жидкости и газа. 2015. Т.50, № 6. С. 828-840.
28. Kozelkov A.S., Kurkin A.A., Pelinovsky E.N., Tyatyushkina E.S., Kurulin V.V., Tarasova N. V. Landslide-type tsunami modelling based on the Navier-Stokes Equations // Science of tsunami Hazards, Journal of Tsunami Society International. 2016. V. 35, № 3. P. 106-144.
29. Козелков A.C. Методика численного моделирования цунами оползневого типа на основе уравнений Навье-Стокса // Вычислительная механика сплошных сред. 2016. Т. 9, № 2. С. 218-236.
30. Погосян М.А., Савельевских Е.П., Шагалиев P.M., Стрелец Д.Ю., Рябов A.A., Корнев A.B., Дерюгин Ю.Н., Спиридонов В.Ф., Циберев К.В. Применение отечественных суперкомпьютерных технологий для создания перспективных образцов авиационной техники // Журнал ВАНТ, сер. Математическое моделирование физических процессов. 2013. Вып. 2. С. 3-18.
31. Погосян М.А., Савелъских Е.П., Стрелец Д.Ю., Корнев A.B. Отечественный суперкомпьютерные технологии в авиационной промышленности // Наука и технологии в промышленности. 2012. № 2. С. 26-35.
32. Волков К.Н., Дерюгин Ю.Н., Емельянов В.Н., Карпенко А.Г., Козелков A.C., Тетерина, И.В. Методы ускорения газодинамических расчетов на неструктурированных сетках. М.: Физматлит, 2013. С. 536, ISBN 978-5-9221-1542-1.
33. Козелков A.C., Шагалиев P.M., Курулин В.В., Ялозо A.B., Лашкин С.В. Исследование потенциала суперкомпьютеров для масштабируемого численного моделирования задач гидродинамики в индустриальных приложениях // Журнал вычислительная математика и математическая физика. 2016. Т.56, № 8. С. 154-165.
34. Козелков A.C., Куркин A.A., Шарипова И.Л., Курулин В.В., Пелиновский E.H., Тятюшкина Е.С., Мелешкина, Д.П., Лашкин С.В., Тарасова Н.В. Минимальный базис задач валидации методов расчета течений со свободной поверхностью // Труды Нижегородского государственного технического университета им. P.E. Алексеева. 2015. № 2(109). С. 49-69.
References
1. Butler Т.М. Advance of LINC method. Proceeding of the Second International Conference on Numerical Methods in Fluid Dynamics. New York. 1971.
2. Burago N.G. Chislennoe reshenie zadach MSS s podviznvmi granichami razdela: PhD. N.G. Burago. M.: RAN, 2003. 222 p. (in Russian).
3. Lucy L.B. A numerical approach to the testing of the fission hypothesis. Astron. J. 1977. V. 82 (12). P 1013-1024.
4. Harlow F.H., Welch J.E. Numerical calculation of time-dependent viscous incompressible flow. Phvs. Fluids. 1965. V. 8. P. 2182-2189.
5. Daly B.J. A technique for including surface tension effects in hydrodvnamic calculations. J. Comput. Phvs. 1969. V. 4. P. 97-117.
6. Hirt С. W., Nichols B.D. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of computational physics. 1981. V. 39. P. 201-226.
7. Ubbink O. Numerical prediction of two fluid systems with sharp interfaces. PhD thesis, Imperial College. University of London. 1997.
8. Kozelkov A.S., Meleshkina D.P., Kurkin A.A., Tarasova N. V., Lashkin S. V., Kurulin V. V. Fully implicit method for solution of Navier-Stokes equations for simulation of multiphase flows with free surface. Journal Computational Technologues. 2016. V. 21, N 5. P. 54-76. (in Russian).
9. Ferziger J.H., Peric M. Computational Method for Fluid Dynamics. Springer-Verlag. New York, 2002.
10. Fletcher Computation Techniques for Fluid Dynamics. Springer-Verlag: Berlin, Heidelberg, New York, London, Paris, Tokyo.
11. Volkov K.N., Deryugin Yu.N., Emelyanov V.N., Kozelkov A.S., Teterina I.V. Raznostnie shemi v zadachakh gazovoi dinamiki na nestrukturirovannikh setkhakh. M.: Fizmavlit, 2014. P. 416. (in Russian).
12. Leonard B.P. The ULTIMATE conservative difference scheme applied to unsteady onedimensional advection. Comput. Methods Appl. Mech. Engrg. 1991. V. 88. P. 17-74.
13. Jasak H., Weller H.G., Gasman A.D. High resolution NVD differencing scheme for arbitrarily unstructured meshes. International journal for numerical methods in fluids. 1999. V. 31. P. 431-449.
14. Muzaferija S., Peric M., Sames P., Schelin T. A twofluid Navier-Stokes solver to simulate water entry. 22 Svmp, on Naval Hydrodynamics, 1998.
15. Waclawczyk T., Koronowicz T. Remarks on prediction of wave drag using VOF method with interface capturing approach. Archives of civil and mechanical engineering. 2008. V. 8. P. 5-14.
16. Jasak H. Error Analysis and Estimation for the finite volume method with applications to fluid flows. Thesis submitted for the degree of doctor. Department of Mechanical Engineering, Imperial College of Science, 1996.
17. Khrabry A.I., Smirnov E.M., Zaytsev D.K. Solving the Convective Transport Equation with Several High-Resolution Finite Volume Schemes: Test Computations. Computational Fluid Dynamics. 2010. P. 535-540.
18. Rusche H. Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions. PhD thesis, Imperial College. University of London, 2002.
19. Betelin V.B., Shagaliev R.M., Aksenov S. V, Belyakov I.M., Deryuguin Yu.N., Kozelkov A.S., Korchazhkin D.A., Nikitin V.F., Sarazov A.V., Zelenskiy D.K. Mathematical simulation of hvdrogen-oxvgen combustion in rocket engines using LOGOS code. Acta Astronautica. 2014. V. 96. P. 53-64.
20. Boiko A.V., Nechepurenko Yu.M., Zhuchkov R.N., Kozelkov A.S. Laminar-turbulent transition prediction module for LOGOS package. Thermophvsics and Aeromechanics. 2014. V. 21, N 2. P. 201-220. (in Russian).
21. Deryugin Yu.N., Zhuchkov R.N., Zelenskiy D.K., Kozelkov A.S., Sarazov A.V., Kudimov N.F., Lipnickiy Yu.M., Panasenko A.V.,Safronov A.V. Validation results for the LOGOS multifunction software package in solving problems of aerodynamics and gas dynamics for the lift-off and injection of launch vehicles. Mathematical Models and Computer Simulations. 2014. V. 26, N 9. P. 83-95. (in Russian).
22. Kozelkov A., Kurulin V, Emelyanov V, Tyatyushkina E., Volkov K. Comparison of convective flux discretization schemes in detached-eddv simulation of turbulent flows on unstructured meshes. Journal of Scientific Computing. 2015. N 89. P. 1-16.
23. Kozelkov A.S., Kurulin V.V. Eddy resolving numerical scheme for simulation of turbulent incompressible flows. Vvcisl. matem. i mat. fizika - Computational Mathematics and Mathematical Physics. 2015. V. 55, N 7. P. 135-146. (in Russian).
24. Kozelkov A.S., Krutyakova O.L., Kurkin A.A., Kurulin V.V., Tyatyushkina E.S. Zonal RANS-LES Approach Based on an Algebraic Reynolds Stress Model. Fluid Dynamics. 2015. V. 50, N 5. P. 24-33. (in Russian).
25. Kozelkov A.S., Kurkin A.A., Pelinovskii E.N., Kurulin V.V. Modeling the Cosmogenic Tsunami within the Framework of the Navier-Stokes Equations with Sources of Different Types. Fluid Dynamics. 2015. V. 50, N 2. P. 306-313. (in Russian).
26. Kozelkov A.S., Kurkin A.A., Pelinovskii E.N. Effect of the Angle of Water Entry of a Body on the Generated Wave Heights. Fluid Dynamics. 2016. V. 51, N 2. P. 288-298. (in Russian).
27. Kozelkov A.S., Kurkin A.A., Pelinovskii E.N., Kurulin V.V, Tyatyushkina E.S. Modeling the Disturbances in the Lake Chebarkul Caused by the Fall of the Meteorite in 2013. Fluid Dynamics. 2015. V. 50, N 6. P. 828-840. (in Russian).
28. Kozelkov A.S., Kurkin A.A., Pelinovsky E.N., Tyatyushkina E.S., Kurulin V.V, Tarasova N. V. Landslide-type tsunami modelling based on the Navier-Stokes Equations. Science of tsunami Hazards, Journal of Tsunami Society International. 2016. V. 35, N 3. P. 106-144.
29. Kozelkov A.S. A numerical technology for landslide tsunami simulations based on Navier-Stokes equations. Computational continuum mechanics. 2016. V. 9, N 2. P. 218-236. (in Russian).
30. Pogosyan M.A., Savel'evskikh E.P., Shagaliev R.M., Strelets D.Yu,., Ryabov A.A., Kornev A.V, Deryugin Yu.N., Spiridonov V.F., Tsiberev K.V. Application of Russian
Supercomputer technologies to develop the advanced models of aviation technology. J. VANT. Ser. Mat. Mod. Fiz. Proc. 2013. V. 2. P. 3-18. (in Russian).
31. Pogosyan M.A., Savel'evskikh E.P., Strelets D.Yu., Kornev A.V. Supercomputer technologies in aircraft industry. Nauka i tehnologii v promvshlennosti. 2012. N 2. P. 26-35. (in Russian).
32. Volkov K.N., Deryugin Yu.N., Emelyanov V.N., Karpenko A.G., Kozelkov A.S., Teterina I. KMetody uskoreniva gasodinamicheskikh raschetov na nestructurirovannvkh setkakh. M.: Phismatlit* 2013. P. 536, ISBN 978-5-9221-1542-1. (in Russian).
33. Kozelkov A.S., Shagaliev R.M., Kurulin V.V., Yalozo A.V., Lashkin S.V. Investigation of Supercomputer Capabilities for the Scalable Numerical Simulation of Computational Fluid Dynamics Problems in Industrial Applications. USSR Comput. Math, and Math. Phvs. 2016. V. 56, N 8. P. 154-165. (in Russian).
34. Kozelkov A.S., Kurkin A.A., Sharipova I.L., Kurulin V.V., Pelinovsky E.N., Tyatyushki-na E.S., Meleshkina D.P., Lashkin S.V., Tarasova N.V. Minimal basis tasks for validation of methods of calculation of flows with free surfaces. Transactions of Nizhni Novgorod State Technical University. 2015. N 2(109). P. 49-69. (in Russian).
Flocmynujia e pedaK'n'um 05.06.2017
