CC BY
20
ТЕОРЕТИЧЕСКАЯ И МАТЕМАТИЧЕСКАЯ ФИЗИКА
Численное моделирование волн Фарадея на основе уравнений гидродинамики в приближении мелкой воды
Т. Г. Елизарова1,0, Д. С. Сабурин2,6
1 Институт прикладной математики имени М. В. Келдыша РАН.
Россия, 125047, Москва, Миусская пл., д. 4.
2 Московский государственный университет имени М.В. Ломоносова, физический факультет, кафедра математики. Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2.
E-mail: а telizar@mail.ru, 6 saburin@physics.msu.ru
Статья поступила 17.10.2014, подписана в печать 11.11.2014.
В рамках уравнений гидродинамики, записанных в приближении мелкой воды, проведено численное моделирование генерации волн Фарадея. Вычислительный алгоритм построен на основе регуляризованных уравнений мелкой воды. Полученные численные результаты соответствуют имеющимся теоретическим и экспериментальным данным. В частности, получены частотные диапазоны возбуждения первой и второй моды колебаний, в рамках которых наблюдаются регулярные, нерегулярные и разрушающиеся волны в зависимости от величины вызывающей их перегрузки.
Ключевые слова: волны Фарадея, уравнения мелкой воды, численное моделирование.
УДК: 517.958:533.7. PACS: 47.35.Bb, 02.60.Cb, 47.11.Df.
Введение
Волнами Фарадея называют стоячие волны, возбуждаемые в сосуде с жидкостью, который совершает периодические колебания в направлении силы тяжести. Частота этих волн кратна половине частоты колебаний сосуда с жидкостью. Впервые данное явление было обнаружено в экспериментах Фарадея в 1831 г. и с тех пор интенсивно изучается. Экспериментально установлено, что даже при малых амплитудах движения сосуда высота возбуждаемых волн достигает больших значений в нескольких узких диапазонах частот внешнего воздействия. Подобные колебания наблюдаются как в природных процессах, так и в технических устройствах и при транспортировке жидкостей.
В качестве примеров исследований волн Фарадея укажем, в частности, [6, 7]. В настоящее время изучение колебательных движений типа волн Фарадея проводятся на основе экспериментов и теоретических моделей [5-8]. Теоретические подходы основаны на анализе решений упрощенных модельных уравнений, либо на численном решении нестационарных уравнений Навье-Стокса, дополненных уравнением для описания свободной поверхности жидкости (см., например, [4] и библиографию к этой статье). Оба эти подхода представляют существенные трудности. Модельные нелинейные уравнения ограничены в постановках и не позволяют учитывать многомерность задачи, форму сосуда и другие особенности течения. Решение полных уравнений гидродинамики достаточно трудоемко и примеры численного моделирования указанных задач весьма ограничены. Поэтому поиск новых эффективных подходов к численному моделированию волн Фарадея является актуальным.
В настоящей работе авторы впервые проводят численное моделирование возбуждения колебаний жидкости, возникающих под действием периодического внеш-
него воздействия в вертикальном направлении, с помощью уравнений гидродинамики в приближении мелкой воды (МВ) [9, 10]. Модель МВ описывает течения вязкой несжимаемой жидкости со свободной поверхностью в сосудах произвольной формы. Приближение уравнений МВ существенно проще полной системы уравнений Навье-Стокса, но имеет свои ограничения: в модели МВ предполагается, что вертикальная компонента скорости течения несущественна и по ней выполнено осреднение.
В большинстве расчетов мы рассматриваем движение жидкости в плоскости симметрии сосуда. Задача решается в неинерциальной системе координат. Расчеты проводятся на основе регуляризованных уравнений мелкой воды [11, 12], которые в свою очередь формально могут рассматриваться как вариант баротропного приближения изученных ранее квазигазодинамических уравнений [13].
1. Постановка задачи и метод численного решения
Рассмотрим сосуд с жидкостью, совершающий вертикальные колебания (рис. 1). Закон смещения сосуда по времени имеет вид S(t) = s cos(Üt), где s — амплитуда колебаний сосуда. В неинерциальной системе координат, жестко связанной с сосудом, вертикальное ускорение g, действующее на жидкость, примет вид
g(t)= go(l + SЮ" cos(^))
(1)
где §о=980 см/с2 — ускорение силы тяжести. Величину Я = sП2/g0 называют перегрузкой. Гармонические колебания сосуда с частотой П могут вызвать незатухающие колебания жидкости с частотой ип, если выполняется соотношение
шП
ип = пр (2)
я сое СИ
Рис. 1. Иллюстрация к постановке задачи о формировании волн Фарадея
где п — номер моды гармонической волны, т — целое число. Область частот, в которой возбуждаются колебания при частотах шп ~ П/2, соответствует основному частотному диапазону возмущения. Таким образом, первая мода с частотой ш1 соответствует частоте колебаний сосуда П = 2ш1. Вторая мода соответствует частоте П = ш2. Из экспериментов известно, что для образования волн Фарадея амплитуда колебаний сосуда 5 не должна быть слишком малой, но в то же время величина перегрузки Я не должна превышать единицу, Я < 1.
Длина волны Хп, волновое число £ и частота колебаний шп для приближения МВ связаны следующим образом:
пп
Хп = 2Ь/п, £ = , 2пс ппс ¡—г
шп = — = , с = ^/gh,
(3)
(4)
где Ь — горизонтальный размер сосуда, h — высота уровня жидкости в сосуде, с — скорость распространения малых возмущений.
В экспериментах [1-3] наблюдаются три категории поверхностных волн Фарадея:
1) регулярные волны, т. е. волны, сохраняющие периодичность во времени и пространственную симметрию;
2) нерегулярные волны, т. е. волны, сохраняющие связность объема и имеющие нарушение временной и пространственной симметрии;
3) разрушающиеся волны: волны, сопровождающиеся срывом капель и струй со свободной поверхности.
Будем описывать генерацию волн Фарадея в рамках приближения МВ с помощью следующих уравнений:
дН .. , _ — + шУ hu = 0, д1
д(Ни) , ,, ч „gh2 , ,
-д^- + й\у(Ни ® и) + = -(и\и\.
(5)
(6)
Неизвестными величинами в системе (5)-(6) являются Н(х, — высота уровня жидкости и и(х, — скорость. Здесь (л — безразмерный коэффициент гидравлического трения. Задача рассматривается в неинерциальной системе координат, связанной с сосудом, в которой ё = ё(1) — переменное значение ускорения силы тяжести (1). На границах расчетной области ставятся условия непротекания для скорости и отражения для
высоты слоя жидкости
и = 0, дх1 = 0.
ОХ1
Система уравнений (5)-(6) дополняется соответствующими начальными условиями.
Для численного решения уравнений МВ (5)-(6) будем использовать регуляризованный вид этих уравнений (см., например, [11, 12]). Для плоского одномерного течения система регуляризованных уравнений МВ имеет вид
дН + дт=0
31 дх
дНи д]ти д (дП д1 дх
+ дх(ёг) = дх- (ии\,
(7)
(8)
где
]т = Ни - т),
т
= Т-(дди + , (9)
Н V дх
' дх )
П = тиЩ и
( ди дН\
[и дХ + ё~]
дх
. . дН . ди\ ёН2 ди
дхУ + тёЛи дх + Н дх) + т~Г дх
(10)
Здесь т — параметр регуляризации, или сглаживания. При т = 0 система уравнений (7)-( 10) переходит в классическую систему уравнений мелкой воды (5)-(6). При проведении численных расчетов слагаемые с коэффициентом т рассматриваются как регуляризи-рующие добавки.
Для численного решения регуляризованных уравнений МВ используется явная по времени разностная схема с аппроксимацией всех пространственных производных центральными разностями. Значения искомых переменных Н(х, {) и и(х, отнесены к узлам пространственной сетки. Устойчивость численного алгоритма обеспечивается слагаемыми с коэффициентом т, величина которого связывается с шагом пространственной сетки Дх и вычисляется в виде
д х _
т = а^, с =у/ёНх1), (11)
где 0 < а <1 — численный коэффициент, выбираемый из условий точности и устойчивости счета. Шаг по времени выбирается в соответствии с условием Куранта, которое для данной задачи имеет вид
Дt = в(Дх/с)тп . (12)
Здесь число Куранта 0 < в <1 зависит от величины параметра регуляризации т и подбирается в процессе вычислений для обеспечения монотонности численного решения. Разностный алгоритм решения системы уравнений (7)-(10) приведен, например, в [11]. В расчетах использовались равномерные пространственные сетки с числом ячеек 100, 500 и 1000. Дальнейшее сгущение пространственной сетки не приводило к заметному изменению численных результатов. Выбор оптимальных коэффициентов а и в для близкой задачи о возбуждении колебаний жидкости в резервуаре обсуждался в [12]. В приведенных далее расчетах а = 0.3, в = 0.1, число пространственных ячеек равно 500. Изменение коэффициента а в пределах от 0.1 до 0.5 практически не влияет на результаты расчетов.
2. Результаты численного моделирования
Данные для расчетов выбирались в соответствии с параметрами экспериментов, описанных в [1, 3]. Отметим, что обнаружить резонанс в системе — достаточно трудоемкая задача. В связи с этим были произведены расчеты серии задач для трех начальных значений Н0 уровня заполнения жидкостью прямоугольного сосуда длины Ь = 50 см: Н0 = 15, 30 и 50 см в диапазоне частот возмущения П от 10 до 60 рад/с с шагом ДП = 0.5 рад/с для различных значений перегрузки Я < 1 и с различными величинами амплитуды колебаний сосуда s. Значения s варьировались в пределах от 0.01 до 5 см с шагом 0.1 см. Коэффициент гидравлического трения / выбран равным 0.005. Изменение / в 10 раз не влияло на результаты расчетов. В качестве начального условия задавалась покоящаяся жидкость и(х, I) = 0 с начальным возмущением свободной поверхности в виде стоячей волны первой волновой моды: Н(х, 0) = Й0(1 + 0.25 ^(пх/Ь)).
В результате численных экспериментов для начального уровня заполнения Н0 = 15 см была обнаружена первая мода волн Фарадея в диапазоне частот возмущений от П = 14.5 до 17.5 рад/с и вторая мода в диапазоне частот возмущений от П = 29 до 33.5 рад/с. Для колебаний характерен резкий срыв на интервале ДП = 0.5 рад/с. Эти данные совместно с результатами расчетов для уровней заполнения Н0 = 30 и 50 см приведены в таблице. Резонансные частоты соответствуют оценкам в приближении МВ (4).
Коридор резонансных частот, полученных в численных
экспериментах при Яшах ^ 0.12, и теоретическое значение резонансной частоты в приближении МВ (4)
Эксперимент Н0, см; мода рад/с Пшах 9 рад/с П по ф-ле (4), рад/с s Яшах
15, 1-я 15.1 15.7 15.21 0.40 0.100
15, 2-я 30.3 31.2 30.49 0.12 0.119
30, 1-я 21.3 22.1 21.56 0.20 0.099
30, 2-я 42.9 44.2 43.10 0.06 0.113
50, 1-я 27.5 28.6 27.83 0.12 0.100
50, 2-я 55.3 57.1 55.66 0.035 0.116
тов 1 и 4 подтвердили закономерности, полученные при одномерных расчетах.
Графики наглядно показывают, что в ограниченных частотных диапазонах образуются стационарные колебательные движения жидкости. Вне этих диапазонов заданное начальное возмущение поверхности жидкости быстро затухает. При этом амплитуды колебаний для первой и второй мод (см., например, кривые 1 и 4, 2 и 6 и 3 и 8) существенным образом зависят от величины перегрузки (см. кривые 4 и 5). Амплитуды колебаний зависят также от начального уровня жидкости в резервуаре.
тт
тах
35 30 25 20 15 10 5 0
10
20
30
40
50
60
О, с-1
На рис. 2 приведены полученные в численных экспериментах значения амплитуды колебаний Ншах, возникающих под действием колебаний сосуда с частотой П . Данные приведены для начального уровня жидкости Н0 = 15, 30 и 50 см при различных амплитудах возмущения s, указанных в таблице. Девять графиков построены следующим образом: шесть графиков (1, 2, 3, 5, 7, 9) построены для амплитуд возмущений, при которых значение перегрузки Яшах = Я(Пшах) < 0.11. Эти графики соответствуют первой моде колебаний. Графики 4, 6 и 8 соответствуют вторым модам, построенных для тех s, при которых для первой моды Яшах < 0.12. Таким образом, для каждого значения начального уровня жидкости Н0 приведено три графика: первая и вторая моды, при которой Яшах < 0.12, и вторая мода, при которой Яшах для первой моды <0.12. Например, для Н0 = 15 см, это графики 1, 5 и 4 при амплитудах возмущения s, равных 0.4, 0.12 и 0.4 см соответственно. Двумерные расчеты вариан-
Рис. 2. Полученные в численном эксперименте частотные диапазоны возмущений и соответствующие амплитуды стационарных волн: линии 1, 4 соответствуют первой и второй модам при Н0 = 15 см, s = 0.4 см; линии 2, 6 — первой и второй модам при Н0 = 30 см, s = 0.2 см; линии 3, 8 — первой и второй модам при ^0 = 50 см, s = 0.1 см; линия 5 — ^0 = 15 см, s = 0.12 см; линия 7 — ^0 = 30 см, s = 0.06 см; линия 9 — Н0 = 50 см, s = 0.035 см
Соответствующие диаграммы устойчивости приведены на рис. 3. Для каждой амплитуды возмущения s колебания образуются внутри определенного диапазона частот, выделенного на рисунке. Таким образом, в системе образуются незатухающие волны Фарадея только при условии, что частота возмущения находится внутри границ для своей моды.
Из графиков следует, что для каждой моды существует пороговое значение амплитуды возмущений s* и соответствующее значение перегрузки Я, ниже которой данная мода не может образоваться на всем диапазоне частот. Так, для первой моды при Н0 = 15 см пороговое значение составляет s* = 0.23 см, для второй моды s* = 0.08 см.
Зависимость амплитуды колебаний жидкости от времени для регулярных волн первой моды ( П = = 15.6 рад/с, s = 0.4 см, Я = 0.1, Н0 = 15 см.) показана на рис. 4, 5. Параметры колебательного течения соответствуют линии 1 на рис. 2. С помощью этого графика была подсчитана частота колебаний для волн Фарадея,
к, см
15.4 15.2 15.0
14.8
см
Рис. 3. Полученные в численном эксперименте диаграммы устойчивости: 1 — к = 0.15 (1-я мода), 5* =0.23 см; 2 - к = 0.3 (1-я мода), 5* =0.12 см; 3 — к = 0.5 (1-я мода), 5* = 0.09 см; 4 - к = 0.15 (2-я мода), 5* = 0.08 см; 5 - к = 0.3 (2-я мода), 5* = 0.04 см; 6 — к = 0.5 (2-я мода), 5* = 0.03 см (к — начальный уровень жидкости, 5* — минимальная амплитуда, при которой реализуются стационарные волны. Черные квадраты соответствуют расчетным точкам
наблюдающаяся в численных экспериментах. На всем временном диапазоне она составила ш\ = 7.5-7.8 рад/с. Таким образом, в наших расчетах было получено соотношение П « .
Форма колебаний первой моды меняется в зависимости от величины перегрузки й. При малых значениях й колебания максимально близки к стоячим волнам (регулярные волны, рис. 4, 6), при больших й они имеют вид бегущих волн в замкнутом сосуде (нерегулярные волны, рис. 7). В обоих случаях наблюдается один узел и две пучности, что соответствует первой моде колебаний.
Внешний вид колебаний для второй моды волн Фарадея представлен на рис. 8, 9. При малых й их форма близка к регулярным стоячим волнам (рис. 8). При больших значениях перегрузках й < 1 появляются разрушающиеся волны (рис. 9). Как и в натурном эксперименте, они имеют вид бегущих с двух сторон навстречу друг другу волн, после столкновения которых образуются срывы струй и капель со свободной поверхности. На рис. 9 экстремумы величины /г в силу используемого описания связной свободной поверхности выглядят как пики, однако их ширина составляет порядка 10 шагов сетки, что говорит о том, что механизм образования этих экстремумов не связан с численной неустойчивостью.
Таким образом, основные черты получающихся волновых движений согласуются с теоретическими и экспериментальными данными, описывающими волны
14.6
40
60
80
с
Рис. 4. Временная зависимость амплитуды колебаний поверхности жидкости с возмущающей частотой П = 15.6 рад/с, 5 = 0.4 см и й = 0.1 для к0 = 15 см. Общий вид
60.00 60.41
61.00 61.22
Рис. 5. Фрагмент временной зависимости амплитуды колебаний поверхности жидкости с возмущающей частотой П = 15.6 рад/с, 5 = 0.4 см и /? = 0.1 для ко = 15 см
к, см
15.4 ¿ = 37.8 с ? = 26.1
15.2
15.0
14.8
14.6
0
10
20
30
40 Ь, см
Рис. 6. Поверхности жидкости для первой моды волн Фарадея, ко = 15 см, П = 15.635 рад/с, 5 = 0.4 см, Я = 0.1
ТЕОРЕТИЧЕСКАЯ И МАТЕМАТИЧЕСКАЯ ФИЗИКА
7
А, см
18
16
14
1 - < = 37.5 с
2342 < = 37.7 с < = 37.9 с < = 38.1 с 4
- / / •ч ^ \ / \
/ / \ ч
О
10
20
30
40 Ь, см
Рис. 7. Поверхности жидкости для четырех последовательных моментов времени, первая мода волн Фарадея, ко = 15 см, П = 15.7 рад/с, 5 = 1.6 см, Я = 0.4
1г, см
15.5 -
15.0 -
14.5 -
14.0
Рис. 8. Поверхности жидкости для второй моды волн Фарадея, ко = 50 см, П = 29.4 рад/с, 5 = 0.25 см, Я = 0.22
А, см
150
Л
I- < = 25 с
2 - < = 25.03 с
3 - < = 25.05 с 4- < = 25.12 с
100 -
40 Ь, см
Рис. 9. Поверхности жидкости для четырех последовательных моментов времени, вторая мода волн Фарадея, ко = 50 см, П = 62.5 рад/с, 5 = 0.25 см, Я = 0.98
Фарадея. А именно в численном эксперименте получено, что образуются незатухающие колебания в частотных диапазонах, соответствующих приближению МВ (4). Обнаружено возникновение первой и второй мод волн Фарадея в фиксированных частотных диапазонах, при этом частота колебаний в два раза меньше частоты возмущения. На границе частотных диапазонов колебания имеют резкий срыв <0.5 рад/с. Имеется пороговое значение амплитуды колебаний сосуда 5*, ниже которого стационарные колебания в жидкости не образуются. В зависимости от величины перегрузки й наблюдаются регулярные, нерегулярные и разрушающиеся волны.
Однако частоты возбуждения колебаний в расчетах оказываются выше, чем в эксперименте. В частности, в [1] диапазон частот возмущений П, при которых образуются колебания второй моды для /?о = 15 см при перегрузках й < 0.3, составляет от 20 до 24 рад/с, в то время как, согласно нашим результатам, частота возбуждения лежит в диапазоне от 29 до 33 рад/с (рис. 2, линия 4). По мнению авторов, это объясняется простотой выбранной модели, в которой не учитываются пространственные эффекты. Кроме того, в используемой модели не исследовано влияние сил трения о стенки и дно сосуда.
Заключение
В настоящей работе впервые показана возможность численного моделирования волн Фарадея с помощью уравнений мелкой воды. Движение жидкости предполагалось одномерным. Расчеты были выполнены с помощью регуляризованных уравнений мелкой воды.
Основные черты рассчитанных волновых движений согласуются с теоретическими и экспериментальными данными, описывающими волны Фарадея. В частности, в численном эксперименте было получено формирование незатухающих параметрических колебаний в частотных диапазонах, соответствующих приближению мелкой воды, при этом частота колебаний в два раза меньше частоты возбуждающей силы. Обнаружено возникновение первой и второй мод волн Фарадея в ограниченных частотных диапазонах. В зависимости от величины перегрузки получено формирование регулярных, нерегулярных и разрушающихся волн. Отметим, что полученные в численных расчетах диапазоны частот внешнего воздействия, вызывающих возбуждение колебаний, оказываются выше значений, зафиксированных в [1]. По мнению авторов, улучшение согласия данных расчета и эксперимента может быть достигнуто при использовании более полной вычислительной модели, включающей в себя более адекватное описание эффектов трения о стенки.
Несомненные плюсы предложенного метода — простота его программной реализации, быстрота численного расчета, а также более широкие по сравнению с экспериментом возможности для изучения волн Фарадея в сосудах различной формы с различными вариантами внешнего воздействия. Таким образом, численное моделирование волн Фарадея с помощью уравнений мелкой воды является перспективным подходом к изучению данного явления.
Работа выполнена при финансовой поддержке РФФИ (грант 13-01-00703-а).
Список литературы
1. Калиниченко В.А., Нестеров С.В., Секерж-Зенько-вич С.Я., Чайковский А.А. // Изв. РАН. Механика жидкости и газа. 1995. № 1. С. 122.
2. Калиниченко В.А. // Изв. РАН. МЖГ. 2009. № 4. С. 112.
3. Калиниченко В.А., Сереж-Зенкович С.Я. // Изв. РАН. Механика жидкости и газа. 2010. № 1. С. 128.
4. Perinet N., Juric D., Tuckerman L.S. // J. Fluid Mech. 2009. 635. P. 1.
5. Wright J, You S, Pozrikidis C. // J. Fluid Mech. 2000. 402. P. 1.
6. Болотин В.В. // ПММ. 1956. 20, № 2. C. 61.
7. Моисеев Н.Н. // Матем. сб. 1953. 32(74), № 1. С. 61.
8. Ibrahim R.A. Liquid Sloshing Dynamics: Theory and Applications. Cambridge, 2005.
9. Сретенский Н.Л. Теория волновых движений жидкости. М., 1977.
10. Куликовский А.Г., Погорелое Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М., 2001.
11. Булатов О.В., Елизарова Т.Г. // Журн. вычисл. матем. и матем. физ. 2011 51, № 1. C. 170.
12. Елизарова Т.Г., Сабурин Д.А. // Матем. моделирование. 2013. 25, № 3. C. 75.
13. Елизарова Т.Г. Квазигазодинамические уравнения и методы расчета вязких течений. М., 2007.
The numerical simulation of Faraday waves on the basis of hydrodynamic equations in the shallow water approximation
T.G. Elizarova1a, D.S. Saburin2 6
1M. V. Keldysh Institute for Applied Mathematics, Russian Academy of Sciences, Moscow 125047, Russia. 2 Department of Mathematics, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia.
E-mail: a telizar@mail.ru, 6 saburin@physics.msu.ru.
Faraday waves were numerically simulated in the shallow-water approximation. The algorithm is based on the regularized shallow-water equations. The results correspond to the existing theoretical and experimental data. In particular, the parametric resonance is found for the first and second wave modes, in which regular and irregular waves with respect to the overload value are observed.
Keywords: Faraday waves, shallow-water equations, numerical modeling. PACS: 47.35.Bb, 02.60.Cb, 47.11.Df. Received 17 October 2014.
English version: Moscow University Physics Bulletin 1(2015). Сведения об авторах
1. Елизарова Татьяна Геннадьевна — доктор физ.-мат. наук, гл. науч. сотрудник ИПМ РАН, профессор; e-mail: telizar@mail.ru.
2. Сабурин Дмитрий Сергеевич — аспирант; e-mail: saburin@physics.msu.ru.
