CC BY
18
2019 Математика и механика № 60
МЕХАНИКА
УДК 532.529.5
DOI 10.17223/19988621/60/6
Е.И. Борзенко, Е.И. Хегай
ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ЗАПОЛНЕНИЯ ЁМКОСТИ НЬЮТОНОВСКОЙ ЖИДКОСТЬЮ С ПРИМЕНЕНИЕМ VOF-МЕТОДА1
Проведено исследование процесса заполнения плоской ёмкости с центральным телом в присутствии силы тяжести. Течение описывается уравнениями Навье - Стокса и неразрывности, на стенках выполняются условия прилипания. Жидкость подается через канал на верхней крышке с заданным расходом. Задача решается численно с использованием метода контрольного объема и корректирующей процедуры SIMPLE. Местоположение свободной поверхности определяется с помощью метода PLIC VoF. Выявлены различные режимы заполнения ёмкости в зависимости от определяющих параметров. Получены кинематические характеристики течения.
Ключевые слова: заполнение, ёмкость, вязкая жидкость, свободная поверхность, численное моделирование, алгоритм SIMPLE, метод PLIC VoF, режим.
Во многих промышленных областях важную роль играет течение жидкости со свободной поверхностью, такое, как заполнение емкостей, растекание жидкости, обтекание объектов и т.п. Для исследования используются методы математического моделирования, которые позволяют получить сведения о поведении свободной поверхности жидкости и особенностях кинематических и динамических характеристик потока.
К настоящему времени известны различные методы численного исследования течения жидкости со свободной поверхностью. Обзор этих методов достаточно подробно представлен в [1, 2]. Одним из широко распространённых методов в вычислительной гидродинамике является метод VoF (Volume of Fluid), предложенный в [3]. Дальнейшее развитие этот метод получил в [4]. В этой работе предложена модификация метода, называемая PLIC VoF (Piecewise-Linear Interface Calculation), которая позволяет представить свободную поверхность в виде наклонных отрезков.
Успешное применение PLIC VoF демонстрируется в [5-10]. В [5] с использованием метода PLIC VoF проводится моделирование процесса обрушения дамбы жидкости. Сравнение результатов моделирования с результатами эксперимента показало хорошее согласование. В [6] исследовалось влияние модели турбулентности на поведение жидкости при обтекании плотины. Достоверность вычислительной методики подтверждена сравнением с многочисленными экспериментами. В [7] приводятся результаты моделирования процесса столкновения жидкости
1 Исследование выполнено за счет гранта Российского научного фонда (проект №18-19-00021).
с препятствием, полученные различными методами отслеживания свободной поверхности. В [8] разработана модификация алгоритма РЫС УоБ, позволяющая проводить расчеты с более высокой точностью. В [9], используя РЫС УоБ, проводится исследование течения ньютоновской жидкости через пористую среду. В [10] РЫС УоБ применяется для моделирования процесса заполнения ёмкости пенополиуретанами.
Методы вычислительной гидродинамики применяются для исследования течения при заполнении емкостей. В работе [11] рассматривается медленное течение вязкой жидкости, реализующееся в технологии метода свободного литья на этапах течения жидкости в массопроводе и заливки пресс-форм различной конфигурации. Показано, что монолитность изделий может нарушаться за счет образования газовых включений и внутренних границ раздела в жидкости. Описаны условия, при которых возможно появление таких дефектов: при течении жидкости в массопроводе, при наличии в нем конструктивного элемента типа «диафрагма» и при заполнении пресс-форм в струйном режиме. В [12] моделируется процесс литья под давлением при заполнении пресс-формы вязкой жидкостью. Показано, что такой процесс может проходить с образованием воздушных полостей или линий спая. Это, в свою очередь, приводит к внутренним дефектам готового изделия. При этом существуют отдельные диапазоны значений числа Стокса, обеспечивающие бездефектный режим заполнения.
Авторы [13, 14] исследовали влияние метода литья топливной композиции на производительность и надежность энергетических установок. Согласно проведенным исследованиям, высказано предположение, что распределение порций топливной композиции внутри пресс-формы влияет на локальную скорость горения. В [15, 16] с помощью математического моделирования проводился расчет рабочих характеристик твердотопливного двигателя при различных вариантах литья
топливных композиций. Показано влияние Г? процесса изготовления топливных зарядов на
баллистические свойства твердотопливных ц{) двигателей.
Целью настоящей работы является изучение процесса заполнения ёмкости вязкой жидкостью и выявление особенностей кинематических и динамических характеристик течения.
У
г,
m
i
H
L
Рис. 1. Область решения Fig. 1. Solution domain
Постановка задачи
Рассматривается течение ньютоновской жидкости, реализуемое при заполнении плоской ёмкости с центральным телом в поле силы тяжести. Область решения показана на рис. 1. Жидкость поступает в ёмкость через входное сечение Г2 с заданным расходом. Вдоль плоскости симметрии имеется твердое тело высотой И < Н, и шириной I < Ь. Считается, что в начальный момент заполнен только входной канал.
h
i
x
Математическая постановка включает уравнения движения и неразрывности, записанные в безразмерном виде
„ . dU dUU dVU Л dP
Re |-+-+-I =--+
dt dx dy ) dx
( d 2U d 2U Л
v йх2 dy2 )
(1)
„ , dV dUV dVV Л dP Re |-+-+-I =--+
dt dx dy ) dy
fd 2V 5 2V Л
v dx2 dy2 )
+ W; (2)
U |p2 = 0, V|r =-2 (1 -(2 x - L )2 ), ^ < x < - (4)
U + = 0; (3)
dx dy
Re = <UD, W =
ц цU о
где U и V - проекции вектора скорости на оси декартовой системы координат x и y, соответственно; P - давление; Re - число Рейнольдса; W - параметр, характеризующий соотношение гравитационных и вязких сил в потоке.
При обезразмеривании выбраны следующие масштабы: D - ширина входного канала; U0 - среднерасходная скорость жидкости во входном сечении; ц - вязкость жидкости; р - плотность жидкости; g - ускорение свободного падения.
В качестве граничных условий во входном сечении Г2 задан параболический профиль, соответствующий установившемуся течению ньютоновской жидкости в плоском бесконечном канале
f(,-(2x - L )2), ^ < x < ^
На твердых стенках Г выполняется условие прилипания
U г1 = 0, V Г, = 0- (5)
На свободной поверхности должны выполняться условия отсутствия касательных напряжений и равенство нормального напряжения внешнему давлению.
В начальный момент времени свободная граница является прямой y = H.
Метод решения
Для численного решения сформулированной задачи используется вычислительная технология VoF-метода. Основные уравнения решаются методом контрольного объема, с применением корректирующей процедуры SIMPLE. Положение свободной поверхности в каждый момент времени определятся с помощью VoF-метода.
Идея метода заключается в том, что вводится функция F, значение которой равно единице во всех точках, занятых жидкостью, и равно нулю в остальных точках. При осреднении по контрольному объёму разностной сетки среднее значение F равно объёмной доле контрольного объёма, занятой жидкостью. В частности, F = 1 соответствует контрольному объёму, полностью заполненному жидкостью, а F = 0 указывает на то, что контрольный объём не содержит жидкости. Контрольные объёмы, в которых F принимает значение от нуля до единицы, содержат свободную поверхность.
Значения этой функции с течением времени находятся из уравнения
dF dF dF
— + U — + V— = 0, (6)
dt dx dy
которое дискретизируется на разностной сетке, и значения потоков через грани рассчитываются, исходя из геометрических соображений и рассчитанных значений скорости на гранях контрольного объёма.
Нормаль к свободной поверхности направлена вдоль градиента функции Е. Когда известно значение Е и направление нормали к свободной поверхности в граничном контрольном объёме, можно провести линию, аппроксимирующую свободную поверхность. Поэтому, помимо определения граничного контрольного объёма, функция Е используется для определения местоположения жидкости внутри него. Классический метод УоБ предполагает существование только двух вариантов расположения свободной поверхности внутри контрольного объёма (вертикально или горизонтально). При решении данной задачи использовалась модификация оригинального метода - РЫС УоБ, которая предполагает произвольный наклон свободной поверхности [19].
Методические расчеты
Для тестирования разработанного алгоритма и программы расчета проведена проверка аппроксимационной сходимости на последовательности сеток. В качестве контролируемых характеристик выбраны форма свободной поверхности и закон сохранения массы жидкости. Используемые значения определяющих параметров Яе = 1 и W = 20 обеспечивают картину течения, представленную на рис. 2 в момент времени t = 4.
2
1
Рис. 2. Форма свободной поверхности на сетках с шагом по пространству h,
в момент времени t = 4 при Re = 1, W = 20: 1 - h = —, 2 - h = —, 3 - h = — и 4 - h = —
10 20 40 80
Fig. 2. Free surface shape calculated using the grid step, h, at time instant of t = 4
for Re = 1 and W = 20: h = (1) —, (2) —, (3) —, and (4) —
10 20 40 80
4
3
0
Ошибка в законе сохранения массы жидкости рассчитывалась по формуле
Мо (Г)- М (Г)
E =
•100%.
Мо «)
Здесь М0 (/) - масса жидкости, которая поступила в ёмкость через входное отверстие за время М (/) - вычисленная масса жидкости в ёмкости в момент времени
Максимальное различие формы свободной поверхности на ряде сеток наблюдается на линии трехфазного контакта жидкости и центрального тела. Различие координаты формы свободной поверхности в узлах, лежащих на твердой стенке центрального тела, вычислялось по формуле Е2 = |уИ - у^2| -100%. уИ и у^2 -
координаты у точек свободной поверхности в сечении х = 2 на сетках с шагами И и И/2 соответственно.
Различие формы свободной поверхности по ширине стекающего слоя рассчитывалось по формуле Е3 = |хИ - хц2| -100%. хИ и хц2 - координаты х точек свободной поверхности в сечении у = 3 на сетках с шагами И и И/2 соответственно.
В таблице приведены значения Е1, Е2 и Е3 на последовательности сеток. Максимальная ошибка в законе сохранения массы жидкости Е1 получалась в момент времени, соответствующий полному заполнению. Е2 и Е3 вычислялись при параметрах, для которых реализуется ситуация, показанная на рис. 2.
Вычисленные значения Е1, Е2 и Е3 на сетках с шагом Н по пространству, %
h 1 10 1 20 1 40 1 80
E, 0.359 0.081 0.030 0.008
E, 0.060 0.038 0.013
E3 0.00207 0.00097 0.00020
Результаты, представленные в таблице, демонстрируют аппроксимационную сходимость по выбранным величинам.
Аналогичная задача о заполнении ёмкости рассматривалась в [11] в приближении ползущего течения, а решение получено методом граничных элементов. На рис. 3 демонстрируется сравнение результатов, полученных методом РЫС УоБ (рис. 3, Ь, сС), с данными [11] (рис. 3, а, с). Наблюдается качественное согласование с небольшими количественными отклонениями.
Рис. 3. Эволюция свободной поверхности при H = 7.5, h = 2.5, L = 5, и l = 0.5: a - ползущее течение, W = 0.04 [11], b - Re = 0.1, W = 0.04, c - ползущее течение, W = 0.4 [11], d - Re = 0.1, W = 0.4 Fig. 3. Free surface evolution at H = 7.5, h = 2.5, L = 5, and l = 0.5: (a) a creeping flow at W = 0.04 [11], (b) Re = 0.1 and W = 0.04, (c) a creeping flow at W = 0.4 [11], and (d) Re = 0.1 and W = 0.4
В [17] исследовалась устойчивость струи высоковязкой жидкости, натекающей на горизонтальную твердую поверхность (рис. 4). Приведен график (рис. 5) зависимости критической высоты входного канала над горизонтальной твердой стенкой Н от параметра W, при превышении которой происходит потеря устойчивости струи, выражающаяся в ее периодическом изгибании. При параметрах, соответствующих точкам, лежащим выше кривой 1, струя неустойчивая (рис. 4, а), ниже -устойчивая (рис. 4, Ь). Сравнение показывает согласование результатов.
H Д -1
20 А Л л 3 4
Лл Д А А
10 Д О** о N
и О о ^N^; <>о°о°
0 1 1 1 1 1 1 1 1 1 ......... 1 1 1 1 1 1 1 1 1
-2
-1
0 lg W
Рис. 4. Эволюция свободной поверхности в случаях неустойчивого и устойчивого режимов Re = 0.1, W = 0.1: а - H = 15 и b - H = 10 Fig. 4. Free surface evolution in the cases of unstable and stable regimes at Re = 0.1, W = 0.1: H = (а) 15 and (b) 10
Рис. 5. Зависимость критической высоты H от параметра W: 1 - данные из [17], 2 - результаты настоящей работы при Re = 0.1, 3 и 4 - экспериментальные данные из [18] и [19]
Fig. 5. Critical height H as a function of the parameter W: (1) data from [17], (2) results from the current work at Re = 0.1, (3) and (4) experimental data from [18] and [19], respectively
Результаты
Все дальнейшие расчеты проводились при следующих размерах: Ь = 3, Н = 5, I = 1, И = 4; Н0 выбиралась достаточной для исключения влияния входного сечения на характер течения внутри ёмкости. Давление внутри ёмкости в начальный момент времени считается равным нулю.
В ходе проведения параметрических исследований выявлены четыре режима заполнения, которые качественно отличаются эволюцией свободной поверхности, и продемонстрированы на рис. 6.
При доминировании вязких эффектов над гравитационными наблюдается режим сплошного заполнения, когда свободная поверхность полностью перекрывает боковые полости (рис. 6, а). Увеличение гравитационных эффектов приводит к режиму стекающего по центральному телу слоя (рис. 6, Ь). В этом случае жидкость стекает по стенке центрального тела, достигает дна, а затем заполняется остальная часть ёмкости. Увеличение инерционных эффектов приводит к формированию струи (рис. 6, с), и, в дальнейшем - к режиму стекания по внешним стенкам боковых полостей (рис. 6, а).
Рис. 6. Четыре основных режима заполнения ёмкости: a - сплошное заполнение (Re = 5, W = 0.1), b - стекание по центральному телу (Re = 1, W = 10), c - струя (Re = 6, W = 10), d - стекание по боковым стенкам (Re = 15, W = 10) Fig. 6. Four basic regimes of a reservoir filling: (a) draining all-over the lateral grooves (Re = 5, W = 0.1), (b) draining down the central rod (Re = 1, W = 10), (c) draining as a jet (Re = 6, W = 10), and (d) draining down the lateral grooves (Re = 15, W = 10)
На рис. 7 показано распределение кинематических и динамических характеристик в момент времени t = 8 в режиме сплошного заполнения. По мере заполнения в местах обтекания центрального тела образуется двумерное течение, а вдоль боковых полостей - одномерное. Давление во входном канале увеличивается по
Рис. 7. Кинематика течения в режиме сплошного заполнения (Re = 1, W = 1): U (a), V(b) и P (c) Fig. 7. Flow kinematics in the regime of fluid draining all-over the lateral grooves (Re = 1, W = 1): (a) U, (b) V, and (c) P
мере заполнения. В боковых полостях, вдали от фронта свободной поверхности, изолинии горизонтальны и параллельны друг другу, что соответствует одномерному течению.
На рис. 8 показано распределение кинематических и динамических характеристик в момент времени t = 8 в режиме стекания по центральному телу. В этом случае при заполнении боковых полостей образуются застойные зоны у дна ёмкости, в которых скорость практически равна нулю и жидкость покоится.
Рис. 8. Кинематика течения в режиме стекания по центральному телу
(Re = 1, W = 10): U (a), V (b) и P (с) Fig. 8. Flow kinematics in the regime of fluid draining down the central rod (Re = 1, W = 10): (a) U, (b) V, and (с) P
Поля скоростей и давления в момент времени t = 8 в режиме струи и в режиме стекания по боковым стенкам ёмкости показаны на рис. 9 и 10 соответственно. Характер распределения полей продольной и поперечной скоростей и давления в этих режимах качественно совпадает с предыдущим.
Исследование массораспределения жидкости в процессе заполнения ёмкости позволяет оценивать неоднородность пространственного распределения свойств материала конечного изделия. Поэтому в настоящей работе проведено исследование распределения порций жидкости, подаваемых через входное отверстие в пресс-форму за единичный интервал времени. Топограммы распределения порций для четырех режимов приведены на рис. 11 для момента времени, соответствующего полному заполнению емкости. Видно, что в первом режиме заполнение имеет фонтанирующий характер [20], а в остальных режимах порции жидкости распределяются параллельными слоями.
Анализ истории деформации элементов жидкости позволяет прогнозировать качество конечного изделия. В связи с этим проведено исследование эволюции и ориентации объёмов жидкости квадратной формы, подаваемых через входной канал.
Рис. 9. Кинематика течения в режиме струи (Re = 6, W = 10): U (а), V (b) и P (с) Fig. 9. Flow kinematics in the regime of fluid draining as a jet (Re = 6, W = 10): (а) U, (b) V, and (с) P
с
с N
i t
—ъо^ f J&: K^fsZ
1— 35---1 L-35—J
1—40-1 1-40-'
Рис. 10. Кинематика течения в режиме стекания по боковым стенкам (Re = 15, W = 10): U(a), V(b) и P (c) Fig. 10. Flow kinematics in the regime of fluid draining down the lateral grooves (Re = 15, W = 10): (a) U, (b) V, and (c) P
Рис. 11. Топограммы распределения порций жидкости: a - сплошное заполнение (Re = 5, W = 0.1), b - стекание по центральному телу (Re = 1, W = 10), c - струя (Re = 6, W = 10), d - стекание по боковым стенкам (Re = 15, W = 10)
Fig. 11. Topograms of the fluid mass distribution: (a) draining all-over the lateral grooves (Re = 5, W = 0.1), (b) draining down the central rod (Re = 1, W = 10), (c) draining as a jet (Re = 6, W = 10), and (d) draining down the lateral grooves (Re = 15, W = 10)
В первом режиме сказывается влияние фонтанирующего характера процесса. При этом объёмы жидкости вытягиваются вдоль боковых полостей. В остальных режимах элементы жидкости меньше подвергаются деформации. Элементы располагаются параллельно. На рис. 12 показаны положения деформированных объёмов жидкости на момент полного заполнения ёмкости.
С/
^ /
1 SJ
1 \J -_>
т
J J J J
Рис. 12. Деформация объёмов жидкости, подаваемых во входной канал через единичный интервал времени поочередно в левую и правую половину: a - сплошное заполнение (Re = 5, W = 0.1), b - стекание по центральному телу (Re = 1, W = 10), c - струя (Re = 6, W = 10), d - стекание по боковым стенкам (Re = 15, W = 10)
Fig. 12. Strain of the fluid portions entering through the inlet channel per unit time alternately into the left and right grooves: (a) draining all-over the lateral grooves (Re = 5, W = 0.1), (b) draining down the central rod (Re = 1, W = 10), (c) draining as a jet (Re = 6, W = 10), and (d) draining down the lateral grooves (Re = 15, W = 10)
На рис. 13 представлена диаграмма режимов заполнения. При малых значениях параметра W процесс заполнения проходит в режиме сплошного заполнения (область I). При малых значениях числа Рейнольдса жидкость стекает по центральному телу (область II). Увеличение Рейнольдса приводит к отрыву потока от центрального тела, процесс переходит в режим струи (область III) и при дальнейшем росте Re в режим стекания по внешним стенкам боковых полостей (область IV).
Re
15
10
Рис. 13. Диаграмма режимов: I - сплошное заполнение, II -стекание по центральному телу, III - струя, IV - стекание по боковым стенкам
Fig. 13. Diagram of the regimes: (I) draining all-over the lateral grooves, (II) draining down the central rod, (III) draining as a jet, and (IV) draining down the lateral grooves
"i—I—I—I——I—I—I—I——I—I—I—I——I—I—I—г
5 10 15 W
5
0
Заключение
Исследован процесс заполнения ёмкости с центральным телом под действием силы тяжести. Алгоритм решения основных уравнений, описывающих течение, основан на совместном использовании метода PLIC VoF и корректирующей процедуры SIMPLE. Параметрические исследования проводились в диапазоне чисел Рейнольдса 0 < Re < 20 и значений параметра W, характеризующего соотношение гравитационных и вязких сил в потоке 0 < W < 20.
В результате продемонстрированы четыре режима заполнения ёмкости. При доминировании вязких эффектов над гравитационными наблюдается режим сплошного заполнения, когда свободная поверхность полностью перекрывает боковые полости. Увеличение гравитационных эффектов приводит к режиму стекающего по центральному телу слоя. В этом случае жидкость стекает по стенке центрального тела, достигает дна, а затем заполняется остальная часть ёмкости. Увеличение инерционных эффектов приводит к формированию струи и в дальнейшем - к режиму стекания по внешним стенкам боковых полостей. Проведен анализ кинематики режимов заполнения. Показано, что в режиме сплошного заполнения реализуется одномерное течение в боковых полостях и двумерное - в окрестности обтекания центрального тела, в остальных режимах образуются застойные зоны в окрестности дна. Исследован характер массораспределения жидкости внутри ёмкости. В режиме сплошного заполнения проявляется фонтанирующий характер течения, а в остальных режимах порции жидкости распределяются параллельными слоями. Проведен анализ деформации элементов жидкости в процессе заполнения. Наибольшим деформациям элементы жидкости подвергаются в режиме сплошного заполнения. Построена диаграмма распределения режимов в зависимости от числа Рейнольдса и параметра W.
ЛИТЕРАТУРА
1. Richard J. McSherry, Ken V. Chua, Thorsten Stoesser. Large eddy simulation of free-surface flows // J. Hydrodynamics. 2017. V. 29(1). P. 1-12. DOI 10.1016/S1001-6058(16)60712-6.
2. Ruben Scardovelli and Stephane Zaleski. Direct numerical simulation of free-surface and interfacial flow // Annu. Rev. Fluid Mech. 1999. V. 31. P. 567-603. DOI 10.1146/annurev. fluid.31.1.567.
3. Nichols B.D., Hirt C.W., Hotchkiss R.S. SOLA-VOF: a solution algorithm for transient fluid flow with multiple free boundaries. Los Alamos Scientific Laboratory Report: LA-8355, 1980. 121 p.
4. Jang W., Jilesen J., Lien F.S., Ji Н. A study on the extension of a VOF/PLIC based method to a curvilinear co-ordinate system // Int. J. Computational Fluid Dynamics. 2008. V. 22(4). P. 241-257. DOI: 10.1080/10618560701880953.
5. Alibek Issakhov, Yeldos Zhandaulet, Aida Nogaeva. Numerical simulation of dam break flow for various forms of the obstacle by VOF method // International Journal of Multiphase Flow. 2018. V. 109. P. 191-206. DOI 10.1016/j.ijmultiphaseflow.2018.08.003.
6. Hargreaves D.M., Morvan H.P., Wright N.G. Validation of the volume of fluid method for free surface calculation: The broad-crested weir // Engineering Applications of Computational Fluid Mechanics. 2007. V. 1(2). P. 136-146. DOI: 10.1080/19942060.2007.11015188.
7. Susann Hansch, Dirk Lucas, Thomas Hohne, Eckhard Krepper, Gustavo Montoya. Comparative simulations of free surface flows using vof-methods and a new approach for multi-scale interfacial structures // Proceedings of the ASME 2013 Fluids Engineering Summer Meeting (FEDSM2013-16104). 2013 DOI: 10.1115/FEDSM2013-16104.
8. Shaswat Saincher, Jyotirmay Banerjee. A redistribution-based volume preserving PLIC-VOF technique // Numerical Heat Transfer. 2014. V. 67. P. 338-362. DOI: 10.1080/ 10407790.2014.950078.
9. Yin X., Zaricos I., Karadimitriou N.K., Raoof A., Hassanizadeh S.M. Direct simulations of two-phase flow experiments of different geometry complexities using Volume-of-Fluid (VOF) method // Chemical Engineering Science. 2018. V. 195. P. 820-827. DOI 10.1016/j.ces.2018.10.029.
10. Mohsen Karimi, Hermes Droghetti, Daniele L. Marchisio. Multiscale modeling of expanding polyurethane foams via computational fluid dynamics and population balance equation // Macromol. Symp. 2016. V. 360. P. 108-122. DOI 10.1002/masy.201500108.
11. Пономарева М.А. Моделирование медленных течений вязкой жидкости со свободной поверхностью / Дис. канд. физ.-мат. наук: 01.02.05 - Механика жидкости, газа и плазмы, 2011. 159 с.
12. Ponomareva MA., Yakutenok VA. Simulation of mold filling by a highly viscous fluid using the 2D indirect boundary element method // Boundary Elements and Other Mesh Reduction Methods XXXVIII. 2015. V. 61. P. 285-296. DOI:10.2495/BEM380231.
13. Toru Shimada, Kazushige Kato, Nobuhiro Sekino, Nobuyuki Tsuboi, Yoshio Seike, Mihoko Fukunaga, Yu Daimon, Hiroshi Hasegawa, Hiroya Asakawa. Advanced Computer Science on Internal Ballistics of Solid Rocket Motors // Trans. JSASS Aerospace Tech. Japan. 2010. V. 8. P. 29-37. DOI: 10.2322/tastj.8.Pa_29.
14. Hasegawa H., Fukunaga M., Kitagawa K., Shimada T. Burning rate anomaly of composite propellant grains // Combustion, Explosion, and Shock Waves. 2013. V. 49. P. 583-592 DOI: 10.1134/S0010508213050109.
15. Le Breton P., Ribereau D. Casting Process Impact on Small-Scale Solid Rocket Motor Ballistic Performance // J. Propulsion and Power. 2002. V. 18. P. 1211-1217. DOI: 10.2514/2.6055.
16. Ribereau D., Le Breton P., Ballereau S. Casting Process Effect on Composite Solid Propellant Burning Rate // AIAA 2001-3946. 2001. DOI: 10.2514/6.2001-3946.
17. Пономарева МА., Шрагер Г.Р., Якутенок В А. Устойчивость плоской струи высоковязкой жидкости, натекающей на горизонтальную твердую плоскость // Механика жидкости и газа. Т. 46(1). 2011. С. 53-61.
18. Cruickshank J.O., Munson B.R. Viscous fluid buckling of plane and axisymmetric jets // J. Fluid Mech. 1981. V. 113. P. 221-239. DOI: 10.1017/S0022112081003467.
19. Cruickshank J.O. Low-Reynolds-number instabilities in stagnating jet flows // J. Fluid Mech. 1988. V. 193. P. 111-127. DOI: 10.1017/S0022112088002071.
20. Борзенко Е.И., Фролов ОЮ, Шрагер Г.Р. Фонтанирующее течение вязкой жидкости при заполнении канала с учетом диссипативного разогрева // Механика жидкости и газа. 2014. Т. 49(1). С. 45-55 DOI: 10.7868/S0040357114060013.
Статья поступила 04.04.2019 г.
Borzenko E.I., Hegaj E.I. (2019) NUMERICAL INVESTIGATION OF THE RESERVOIR FILLING WITH A NEWTONIAN FLUID USING THE VOF-METHOD. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika [Tomsk State University Journal of Mathematics and Mechanics]. 60. pp. 73-86
DOI 10.17223/19988621/60/6
Keywords: filling, reservoir, viscous fluid, free surface, numerical simulation, SIMPLE algorithm, PLIC-VoF method, regimes.
In this paper, the filling of a plane reservoir, having a central rod, with a fluid under gravity is studied. The flow is described by the Navier-Stokes and continuity equations. On the solid walls, the no-slip boundary conditions are satisfied. The fluid with a fixed flow rate enters the reservoir through the channel on the top. The problem is solved numerically using the control volume method and SIMPLE procedure. The free surface location is determined using the PLIC-VOF method. Mathematical formulation of the problem includes two non-dimensional parameters, namely, the Reynolds number (Re) and the ratio of the gravity forces to the viscous forces in a fluid (W).
As a result, four regimes of the reservoir filling have been discovered. When viscous effects dominate over gravitational effects, the regime characterized by draining all-over the lateral grooves is observed, wherein the free surface completely overlaps the lateral grooves. An increase in the gravitational effects leads to the regime characterized by draining down the central rod. In this case, the fluid flows along the central rod walls, then reaches the bottom, and fills in the rest part of the reservoir. Further increase in the inertial effects gives rise to the regime characterized by draining as a jet, and, subsequently, to the regime of draining down the lateral grooves. The kinematic characteristics of the flows are presented. It is shown that in the case of draining allover the lateral grooves, one-dimensional flow occurs in the lateral grooves and two-dimensional flow occurs in the vicinity of the central rod. In other regimes, dead zones are formed in the vicinity of the bottom. The behavior of the fluid mass distribution over the reservoir is studied. In the regime of draining fluid all-over the lateral grooves, a fountain flow is observed, and in other regimes, the portions of the fluid are distributed in parallel layers. The deformation of the fluid portions during the filling is analyzed. The diagrams illustrating the reservoir filling regimes are presented at various values of Reynolds number and parameter W.
Financial support. The research is implemented at the expanses of the Russian Science Foundation (project No. 18-19-00021).
BORZENKO Evgeniy Ivanovich (Candidate of Physics and Mathematics, Associate Professor, Tomsk State University, Tomsk, Russian Federation). E-mail: borzenko@ftf.tsu.ru
HEGAJ Efim Igorevich (Tomsk State University, Tomsk, Russian Federation). E-mail: efim_h@ftf.tsu.ru
REFERENCES
1. McSherry R.J., Chua K.V., Stoesser Th. (2017) Large eddy simulation of free-surface flows. Journal of Hydrodynamics. 29(1). pp. 1-12. DOI: 10.1016/S1001-6058(16)60712-6.
2. Scardovelli R., Zaleski S. (1999) Direct numerical simulation of free-surface and interfacial flow. Annual Review of Fluid Mechanics. 31. pp. 567-603. DOI: 10.1146/annurev.fluid.31.1.567.
3. Nichols B.D., Hirt C.W., Hotchkiss R.S. (1980) SOLA-VOF: a solution algorithm for transient fluid flow with multiple free boundaries. Los Alamos Scientific Laboratory Report: LA-8355.
4. Jang W., Jilesen J., Lien F.S., Ji R (2008) A study on the extension of a VOF/PLIC based method to a curvilinear co-ordinate system. International Journal of Computational Fluid Dynamics. 22(4). pp. 241-257. DOI: 10.1080/10618560701880953.
5. Issakhov A, Zhandaulet Ye., Nogaeva A. (2018) Numerical simulation of dam break flow for various forms of the obstacle by VOF method. International Journal of Multiphase Flow. 109. pp. 191-206. DOI: 10.1016/j.ijmultiphaseflow.2018.08.003.
6. Hargreaves D.M., Morvan H.P., Wright N.G. (2007) Validation of the volume of fluid method for free surface calculation: the broad-crested weir. Engineering Applications of Computational Fluid Mechanics. 1(2). pp. 136-146. DOI: 10.1080/19942060.2007.11015188.
7. Hansch S., Lucas D., Hohne Th., Krepper E., Montoya G. (2013) Comparative simulations of free surface flows using VOF-methods and a new approach for multi-scale interfacial structures. Proceedings of the ASME 2013 Fluids Engineering Summer Meeting. FEDSM2013-16104. DOI: 10.1115/FEDSM2013-16104.
8. Saincher Sh., Banerjee J. (2014) A redistribution-based volume preserving PLIC-VOF technique. Numerical Heat Transfer. 67. pp. 338-362. DOI: 10.1080/10407790.2014.950078.
9. Yin X., Zaricos I., Karadimitriou N.K., Raoof A., Hassanizadeh S.M. (2018) Direct simulations of two-phase flow experiments of different geometry complexities using Volume-of-Fluid (VOF) method. Chemical Engineering Science. 195. pp. 820-827. DOI: 10.1016/j.ces.2018.10.029.
10. Karimi M., Droghetti H., Marchisio D.L. (2016) Multiscale modeling of expanding polyurethane foams via computational fluid dynamics and population balance equation. Macromo-lecular Symposia. 360. pp. 108-122. DOI: 10.1002/masy.201500108.
11. Ponomareva M.A. (2011) Modelirovanie medlennykh techeniy vyazkoy zhidkosti so svobod-noypoverkhnost'yu [Modeling of the creeping flow of viscous fluid with a free surface]. Dissertation in Physics and Mathematics.
12. Ponomareva M.A., Yakutenok V.A. (2015) Simulation of mold filling by a highly viscous fluid using the 2D indirect boundary element method. Boundary Elements and Other Mesh Reduction MethodsXXXVIII. 61. pp. 285-296. DOI:10.2495/BEM380231.
13. Shimada T., Kato K., Sekino N., Tsuboi N., Seike Y., Fukunaga M., Daimon Y., Hasegawa H., Asakawa H. (2010) Advanced computer science on internal ballistics of solid rocket motors. Transactions of The JSASS, Aerospace Technology Japan. 8. pp. 29-37. DOI: 10.2322/tastj.8.Pa_29.
14. Hasegawa H., Fukunaga M., Kitagawa K., Shimada T. (2013) Burning rate anomaly of composite propellant grains. Combustion, Explosion, and Shock Waves. 49. pp. 583-592 DOI: 10.1134/S0010508213050109.
15. Le Breton P., Ribereau D. (2002) Casting process impact on small-scale solid rocket motor ballistic performance. Journal of Propulsion and Power. 18. pp. 1211-1217. DOI: 10.2514/2.6055.
16. Ribereau D., Le Breton P., Ballereau S. (2001) Casting process effect on composite solid propellant burning rate. American Institute of Aeronautics and Astronautics. 2001-3946. DOI: 10.2514/6.2001-3946.
17. Ponomareva M.A., Shrager G.R., Yakutenok V.A. (2011) Stability of a plane jet of a highly viscous fluid impinging on a horizontal solid wall. Fluid Dynamics. 46(1). pp. 51-58. DOI: 10.1134/S0015462811010054.
18. Cruickshank J.O., Munson B.R. (1981) Viscous fluid buckling of plane and axisymmetric jets. Journal ofFluidMechanics. 113. pp. 221-239. DOI: 10.1017/S0022112081003467.
19. Cruickshank J.O. (1988) Low-Reynolds-number instabilities in stagnating jet flows. Journal ofFluid Mechanics. 193. pp. 111-127. DOI: 10.1017/S0022112088002071.
20. Borzenko E.I., Frolov O.Yu., Shrager G.R. (2014) Fountain viscous fluid flow during filling a channel when taking dissipative warming into account. Fluid Dynamics. 49(1). pp. 37-45. DOI: 10.1134/S0015462814010062.
Received: April 4, 2019
