М ЕХАНИКА ЖИДКОСТИ И ГАЗА № 6 • 2014
УДК 532.517,532.613
ВЛИЯНИЕ ПЕРЕПАДА ПЛОТНОСТИ И ПОВЕРХНОСТНОГО НАТЯЖЕНИЯ НА ПОВЕРХНОСТИ РАЗДЕЛА ТЕКУЧИХ СРЕД НА РАЗВИТИЕ НЕУСТОЙЧИВОСТИ РЭЛЕЯ-ТЕЙЛОРА
© 2014 г. С. Н. ЯКОВЕНКО
Институт теоретической и прикладной механики им. С.А. Христиановича РАН, Новосибирск e-mail:yakovenk@itam.nsc.ru, s.yakovenko@mail.ru
Поступила в редакцию 30.01.2014 г.
Для описания эволюции поверхности раздела двух несмешивающихся сред вводится уравнение для функции объемной фракции, а эффект кривизны поверхности раздела учитывается с помощью "континуальной" модели силы поверхностного натяжения. Приведены результаты численного решения задачи неустойчивости Рэлея—Тейлора с различным перепадом плотности (Pi/P2) на поверхности раздела, в том числе и для реальных случаев, где имеются данные измерений. Эволюция неустойчивости на начальном этапе не зависит от р1/р2 и согласуется с линейной теорией Тейлора, затем (при р1/р2 < 5) наблюдается спиралевидное развитие неустойчивости Кельвина—Гельмгольца. При р1/р2 < 2 картина развития неустойчивости сохраняется симметричной до значительных моментов времени, тогда как с ростом р1/р2 усиливается асимметрия. Поверхностное натяжение, как и вязкость, приводит к демпфированию роста возмущений неустойчивости Рэлея—Тейлора и вторичных мелкомасштабных нерегулярностей поверхности раздела.
Ключевые слова: неустойчивость Рэлея—Тейлора, поверхность раздела сред, численное моделирование.
Одно из наиболее интересных и интригующих явлений в механике жидкости, газа и плазмы — неустойчивость Рэлея—Тейлора (НРТ), возникающая при возмущении поверхности раздела между слоями тяжелой (сверху) и легкой (снизу) текучих сред, — наблюдается от микромасштабов до астрофизических масштабов, например, в Крабо-видной туманности. Распыление межзвездного газа, выталкиваемого из плоскости галактики магнитными полями и космическими лучами, может развиваться по сценарию НРТ. Подобные процессы имеют место при вспышке сверхновых [1, 2], в искусственных и природных объектах с термоядерным синтезом, при управляемых реакциях горения и пожарах, в мантии, океанах и атмосфере Земли, существенно влияя на климат. Один из эффектов НРТ — возникновение турбулентности и интенсификация смешения за счет конвективного переноса. В недавних Международных конференциях по турбулентному смешению эффектам НРТ уделено значительное внимание (см., например, [3]). Примером актуальности и существенного интереса международного научного сообщества также является сотрудничество крупнейших мировых центров по численному исследованию НРТ [4]. Наибольший из известных расчетов трехмерной НРТ [1] был выполнен при использовании кубической области вычислений с 30723 узлами, для чего потребовалось около двух недель и 65536 процессоров (т.е. полная мощность IBM BlueGene/L, крупнейшего на момент расчетов суперкомпьютера).
Большинство работ по экспериментальным и численным исследованиям НРТ выполнено в классической постановке с двумя слоями несмешивающихся сред или в стратифицированной среде со ступенчатым или линейным профилем плотности (см., например, [5]). Эволюция НРТ в условиях, приближенных к реальным течениям окружающей среды, изучена слабее. Недавно показано [6], что при обрушении внутренних волн в устойчиво стратифицированном потоке также проявляются эффекты НРТ, очевидно инициирующие переход к турбулентности. Обрушение внутренних волн — один из основных источников геофизической турбулентности, однако механизмы развития неустойчивости, возникающей при конвективном опрокидывании волн, до сих пор недостаточно изучены. Для более ясного понимания особенностей развития НРТ остаются актуальными исследования эволюции возмущений на поверхности раздела между двумя слоями разной плотности.
Для описания развития НРТ необходимы адекватные методы разрешения поверхности раздела между областями несмешивающихся сред. Очевидными условиями при этом являются малость вычисляемой конечной толщины поверхности раздела по сравнению с другими характерными масштабами длины, а также эффективность и устойчивость численного алгоритма.
Один из широко используемых методов разрешения поверхности раздела — концепция фракционного объема среды VOF [7], определяющая функцию f объемной фракции, которая равна единице в области, занимаемой более плотной средой, и нулю — в области с менее плотной средой. Вводимое уравнение для f соответствует движению функции объемной фракции вместе с текучей средой. Для аппроксимации адвекции в [7] была развита специальная процедура донор—акцептор, включающая схемы против потока или по потоку в зависимости от ориентации поверхности. Метод VOF применялся в разных работах, например, при моделировании [8] твердотельного вращения и НРТ. Обеспечивая резкую поверхность раздела, этот подход приводит, однако, к генерации ложных изолированных образований вида плавающих осколков [8].
Для уменьшения погрешностей и улучшения описания поверхности раздела предлагались различные модификации VOF-метода и другие подходы, например, применение MUSCL-схемы [9] с интерполянтами высокого порядка. В таких схемах для контроля ложных осцилляций, возникающих около разрывов в методах высокого порядка, вводятся функции-ограничители "уменьшения полных вариаций" (TVD). В частности, TVD-ограничитель в виде "сжимающего" оператора, предложенный в [9], дает приемлемые результаты для распределений как скорости и давления, так и функции f. С другой стороны, поверхность раздела уже не получается такой же резкой, как в VOF-методе, и может иметь толщину порядка нескольких размеров ячеек. Тем не менее, и в этом случае положение поверхности раздела можно однозначно определить, находя точки, где f = 0.5.
Для описания поверхности раздела при малых линейных масштабах движения становится существенным влияние кривизны поверхности раздела и поверхностного натяжения. Один из эффективных в численной реализации подходов учета поверхностного натяжения — "континуальная" модель [10].
В настоящей работе исследуется эволюция НРТ на поверхности раздела текучих не-смешивающихся сред с различным перепадом плотности тяжелой и легкой сред (р1/р2), характеризуемым числом Атвуда A = (р1 — р2)/(р1 + р2), при помощи численных методов, развитых в [9—12]. Вычисленные характеристики сопоставляются с данными линейной теории [13], измерений на поверхностях раздела вода—бензол, вода-воздух [14] и численного моделирования, представленными в других работах. Развитие НРТ - классическая задача механики жидкости и газа, взятая в ряде исследований [7-11, 15-21] для верификации методов разрешения поверхности раздела и учета поверхностного натяжения. В частности, в [15-20] было отмечено влияние перепада плотности двух сред на эволюцию НРТ. Однако моделирование в реальных условиях и
сравнение с данными опыта выполнено лишь в малой части работ, и результаты сравнения (например, в [15]) оказывались зачастую неудовлетворительными.
1. Основные уравнения и численные схемы. В качестве определяющих используются уравнения неразрывности, Навье—Стокса и объемной фракции / для несжимаемого течения с учетом поверхностного натяжения [11]
^ = 0 (1.1)
дх.
дщ + дщик _ 1 dp + 1 д dt dxk рдх/ p dxk
L V
/дщ + дик
дхк дх{)_
- g5a + F (1.2)
P = Plf + P2 (1 - f), ^ = + Ц2 (1 - f)
f + дМ/ = 0 (1.3)
dt dxj
где и, — компоненты вектора скорости, x, — координаты, g — ускорение силы тяжести, ц — динамическая вязкость, р — плотность, p — давление, t — время. Значения р1 и ц соответствуют более тяжелой жидкости при f = 1, тогда как значения р2 и ц2 — более легкой жидкости (газу) при f = 0.
Для решения основных уравнений применена, как и в [11, 12], довольно простая процедура одновременных итераций находимых из (1.1)—(1.2) значений скорости и давления, основанная на методе искусственной сжимаемости, с релаксацией по псевдовремени для удовлетворения уравнения неразрывности. Дифференциальные уравнения (1.1)—(1.3) дискретизируются на смещенной сетке для предотвращения рассогласования полей скорости и давления: компоненты вектора скорости определяются в центрах граней ячеек расчетной сетки, а скалярные величины p, f, р и ц — в центрах ячеек. Для производных по времени и псевдовремени взяты явные схемы первого порядка. Для членов молекулярной диффузии в (1.2) использована центрально-разностная схема.
В [12] приведена иерархия пяти различных способов аппроксимации членов адвекции в (1.2), (1.3) и их верификация в задаче разрушения плотины. В частности, процедура донор—акцептор, применяемая в VOF-методе [7], обеспечивала резкую поверхность раздела, но демонстрировала [12] появление ложных образований в виде "плавающих осколков", как и в [8]. Наиболее приемлемые результаты с достаточно малой толщиной поверхности раздела и ее гладкой формой обеспечивала MUSCL-схема с QUICK-интерполянтами и TVD-ограничителями типа сжимающего оператора, предложенная в [9]. Эта схема использована в настоящей работе для аппроксимации адвекции. Более подробно алгоритмы решения (1.1)—(1.3) (без учета поверхностного натяжения) приведены в [12], где результаты расчетов показали хорошее согласие с данными опытов.
Эффекты поверхностного натяжения вводятся в виде объемных сил вида
F , K = nm =, n =f, [/] = /2 - /1 = 1 (1.4)
P [f ] dxm yjnknk dxt
в уравнения Навье—Стокса (1.2) согласно континуальной модели CSF [10, 11], где ст — коэффициент поверхностного натяжения, к — кривизна поверхности раздела, nt — вектор нормали к поверхности.
Наиболее корректное представление сглаженной функции объемной фракции fs предполагает ее плавное изменение поперек поверхности раздела конечной толщины (в несколько ячеек сетки) за счет свертки f с интерполяционной функцией ядра [10].
Однако функция fs может быть взята в качестве первого приближения равной и самой функции объемной фракции, найденной численно [11] (т.е. fs = f в упрощенной версии CSF-модели), в связи с возможным размазыванием распределений f(x) за счет действия схемной вязкости. Конечно-разностный аналог (1.4) для двумерных течений с границами в виде стенок и плоскостей симметрии дан в [11], вместе с
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.