CC BY
9Математические заметки СВФУ Октябрь—декабрь, 2022. Том 29, № 4
УДК 519.6
МОДЕЛИРОВАНИЕ ДИНАМИКИ ЖИДКОСТИ СО СВОБОДНОЙ ПОВЕРХНОСТЬЮ В ГРАВИТАЦИОННОМ ПОЛЕ СХЕМОЙ КАБАРЕ
Н. А. Афанасьев, В. М. Головизнин, П. А. Майоров, А. В. Соловьев
Аннотация. Предложен явный балансно-характеристический метод КАБАРЕ для расчета динамики жидкости со свободной поверхностью в гравитационном поле в слабосжимаемом приближении. Разработанный метод обладает вторым порядком аппроксимации по времени и пространству, минимальным вычислительным шаблоном в одну пространственно-временную ячейку и минимальной численной вязкостью. Разностная схема тестируется на задачах с различными значениями коэффициента поверхностного натяжения и ускорением свободного падения с разными знаками, в том числе на задаче о развитии неустойчивости Рэлея — Тейлора. Учет сил поверхностного натяжения позволяет избавиться от высокочастотных колебаний на свободной поверхности при расчете неустойчивых задач и регуляризует решение.
DOI: 10.25587/SVFU.2023.23.70.007
Ключевые слова: уравнения гиперболического типа, балансно-характеристиче-ские схемы, смешанные эйлерово-лагранжевы переменные, слабосжимаемая жидкость, неустойчивость Рэлея — Тейлора, поверхностное натяжение.
1. Введение
В последние годы развитие балансно-характеристических схем семейства КАБАРЕ [1] было посвящено решению систем гиперболических уравнений на подвижных расчетных сетках в смешанных эйлерово-лагранжевых переменных и, в частности, решению задач со свободной поверхностью. Потребность в таких методах возникает при решении задач океанологии [2, 3] и задач FSI (fluid-structure interaction) [4,5]. При этом разработанная ранее схема КАБАРЕ на подвижных сетках для уравнений динамики идеального газа [5] позволяла моделировать только малые колебания свободной поверхности газа.
С целью увеличения робастности схемы КАБАРЕ для задач со свободной поверхностью была разработана модификация предложенного в [5] метода на
Работа выполнена при поддержке РНФ (грант 18-11-00163, часть работ по разработке алгоритма для уравнений динамики слабосжимаемой жидкости) и РФФИ (грант 20-3190037, часть работ по учету сил поверхностного натяжения и расчету неустойчивости Рэлея — Тейлора).
© 2022 Афанасьев Н. А., Головизнин В. М., Майоров П. А., Соловьев А. В.
систему уравнений динамики слабосжимаемой жидкости в приближении Бус-синеска. Это позволило построить численную негидростатическую модель динамики жидкости для задач океанологии [6], которая была верифицирована на тестах с умеренными амплитудами колебаний свободной поверхности. Тем не менее для полной верификации предложенного метода требуется провести тесты на задачах, при решении которых ячейки расчетной сетки сильно деформируются. Наиболее подходящим тестом для этого является задача о развитии неустойчивости Рэлея — Тейлора [7].
Рэлей-тейлоровская неустойчивость в основном моделируется в несжимаемом приближении с помощью большого количества различных численных методов, среди которых: методы конечного объема [8], методы гладких частиц [9], вихревые методы в смешанных эйлерово-лагранжевых переменных [10] и многие другие. Среди подходов в слабосжимаемом приближении следует отметить методы ЬаШсе-Во^тапп [11]. С помощью балансно-характеристических методов неустойчивость Рэлея — Тейлора свободной поверхности до сих пор не моделировалась.
В данной работе мы предлагаем адаптацию метода [6] на задачи динамики свободной поверхности в слабосжимаемом приближении с учетом сил поверхностного натяжения. Разработанный метод обладает вторым порядком аппроксимации по времени и пространству, минимальным вычислительным шаблоном в одну пространственно-временную ячейку и минимальной численной вязкостью. Именно малая диссипативность метода требует введения в модель сил поверхностного натяжения, которые позволяют подавить развитие высокочастотных неустойчивостей на свободной поверхности при расчетах неустойчивых задач и регуляризировать решение. Отметим, что для расчетов задач об устойчивых колебаниях свободной поверхности вводить силы поверхностного натяжения, вообще говоря, не надо [6]. Разностная схема тестируется на задачах с различными значениями коэффициента поверхностного натяжения и ускорением свободного падения с разными знаками. Показано, что при моделировании неустойчивости Рэлея — Тейлора ускорение свободной поверхности стремится к ускорению свободного падения. При расчетах устойчивых колебаний поверхности учет сил поверхностного натяжения сглаживает свободную поверхность, при этом почти не изменяя амплитуду колебаний.
Текст организован следующим образом. В разд. 2 обсуждается модель динамики жидкости в приближении Буссинеска для слабосжимаемого случая. В разд. 3 кратко описывается балансно-характеристический метод для моделирования задач о развитии неустойчивости Рэлея — Тейлора. В разд. 4 как в дифференциальную, так и в разностную модель вводятся силы поверхностного натяжения. Тестовые расчеты для различных значений коэффициента поверхностного натяжения и ускорения свободного падения с разными знаками приводятся в разд. 5. Статью завершает разд. 6 с заключительными замечаниями.
2. Модель динамики слабосжимаемой жидкости
Рассмотрим систему двумерных уравнений гидродинамики в приближении Буссинеска с учетом слабосжимаемости жидкости в смешанных эйлерово-лагранжевых переменных [6]:
1 д,7в дви дв(ш - г) _
J д£ дж дз '
1 д-1ви дви2 дв{уо - г)и 1 дбР _
J дж дз р0 дж '
1 д^ад двида — 2) 1 д^Р /р \
+ +-^-- +--= - — - 1 .9, С1)
J д£ дж д^ р0 д^ \р0
1 <9 ^(9 <9р(9и дрв{%и — г)
J д£ дж дз '
= № = с2(0-1),
где (ж, ^) — эйлерова система координат, которая связана с лагранжевой системой координат (ж', 2') якобианом перехода J = д(ж, ^)/д(ж', 2') = дя/дя', р — плотность жидкости, р0 — средняя начальная плотность жидкости, в — безразмерный параметр, показывающий отклонение объема лагранжевой частицы от первоначального, 5Р — приращение давления относительно гидростатического: Р(ж, £) = р0д(Н0 — г) + 5Р(ж, 2, £), Н0 — средняя начальная высота жидкости, (и, ад) — компоненты вектора скорости по направлениям ж и 2, с — искусственная скорость звука, д — ускорение свободного падения. Значение с подбирается таким образом, чтобы |5в| = |в — 1| < 0.01.
Будем предполагать, что система уравнений (1) описывает динамику жидкости в области, ограниченной слева и справа по оси ж и снизу по оси 2 стенками, параллельными данным осям, а сверху по оси 2 располагается свободная поверхность, что накладывает на систему следующие граничные условия:
и(жт1п,-М) = 0, и(жтах,М) = 0, ■ш(ж,2т;п,-£) = 0, (2)
5Р(ж, 2тах, = — р0д(Н) — ¿тах), (3)
<4)
Будем также предполагать, что плотность в начальный момент времени всюду одинакова: р(ж, 2, 0) = р0. Так как плотность в системе уравнений (1) по сути играет роль некоторой пассивно переносимой величины, то и во все моменты времени будет выполняться р(ж, 2, £) = р0, что позволяет не рассматривать закон сохранения массы (четвертое уравнение в (1)), а правую часть в законе сохранения импульса по оси 2 (третье уравнение в (1)) считать равной нулю.
Система уравнений (1) является гиперболической. По каждому из направлений ж, 2 можно выписать ее характеристическую форму (без учета закона
сохранения массы):
д Iх dIx dIz dIz
^T+AX^r = GX> 757 +Лг7Г = С2 dt дх dt dz
IX = töf = (
1
= I u-\--——öü,u--——öü,w
T
y/p0[0]loc ' y/PÖ[6]lo
c c
Ax = diag(Ax), A l=u + ——, A %=u, (5)
Vpö VPÖ v >
Iz = (I/ )T = (■
1
т/Р0Щ1ос"",~ VPÖ\0]lc
= I w H--r^i—öe.w--r^i—öe.u
T
Az = diag(Az), Af = w - z H--\% = w-z--\l=w-z,
Vpö VPÖ
где Iх, Iz — векторы локальных инвариантов Римана [12], Ах, Az — собственные значения по соответствующим направлениям, Gx, Gz — векторы правых частей, выражения для которых нам неважны. Величины с индексами loc в (5) — некоторые локальные величины. Например, при введении разностной сетки эти локальные величины можно считать постоянными в пределах каждой пространственно-временной ячейки. Система уравнений (1) позволяет построить консервативные фазы балансно-характеристической схемы, аппроксимирующие законы сохранения, а система (5) — характеристическую фазу, аппроксимирующую перенос локальных инвариантов Римана системы.
3. Схема КАБАРЕ для представленной модели
Кратко опишем алгоритм схемы КАБАРЕ для системы уравнений (1) с граничными условиями (2)—(4). Данный алгоритм является развитием предложенной ранее схемы для двумерных уравнений газовой динамики [5]. Наиболее полное описание алгоритма представлено в [6].
Пусть имеется некоторая структурированная четырехугольная сетка с прямыми ребрами по оси z и косыми по оси х: u>h = Zi^) \ г = 0, Nx, j = 0, Nz}, Хо = Xmix, xl+i -хг = hl+1/2, zho = ¿i.max, zhnz = ¿i,min = 0. Пусть также имеется неравномерная сетка по времени uiT = {tn | tn+i — tn =r„, n = 0, К — 1}. Зададим в центрах ячеек сетки Uh х ит так называемые консервативные пере-
п п+1/2
менные: р"+1/2 к+1/2 — на целых слоях по времени, ^i+1/2 к+1/2 — на полуцелых слоях по времени. В центрах ребер сетки зададим так называемые потоковые переменные на целых слоях по времени: фП+1/2 к и ФПк+i/i' В качестве консервативных и потоковых переменных зададим полный набор неизвестных рассматриваемой системы: {р, ф} = {в, u, w, SP}.
Общий алгоритм схемы КАБАРЕ при переходе со слоя по времени n на слой n + 1 выглядит следующим образом:
1. Вычисление положения сетки и площадей ячеек на полуцелом слое по времени n +1/2.
Сначала с помощью аппроксимации граничного условия (4) вычисляются скорости центров ребер сетки на свободной границе (первая «балансная» фаза
для свободной границы):
7п _ у п
7П -1ПП -11П 1+1'°
г+1/2,0 ~~ шг+1/2,0 "¿+1/2,0 и ' ^ '
"г+1/2
Затем, учитывая, что уравнение (4) можно рассматривать как нелинейное уравнение переноса с ненулевой правой частью для инварианта Римана 2, проводится экстраполяция ¿-координат узлов на свободной поверхности (вторая характеристическая фаза для свободной границы):
22п—1//2,0 — если °.5(иП-1/2,0 + <+1/2,0) > 0,
2П+1 = { 22 "+1/2
i+1 /2,0 - z"+1,0. если °-^иГ-1/2,0 + <+1/2,j) < 0. (7)
0-5 [2Г-+11//22,0 + 2Г++1/220] ИНаЧе-
Экстраполяция (7) также дополняется монотонизацией на основе принципа максимума [13]. Скорости узлов на свободной границе вычисляются так, чтобы узлы прошли половину пути до координаты, полученной по экстраполяции (7):
zn+1 _
¿■0 = (8)
Скорости внутренних узлов и центров ребер выбираются так, чтобы высота ячеек в каждом слое по оси z была равномерной:
АГ — k N — h _
■п _1_2V_Z;/n -п _ ilV__Zz,n Li (V /Q\
zi,k — jy zi, 0> zi+l/2,k — дг zi+1/2,0' 1 — J-,JVZ. ^
Окончательно с помощью вычисленных скоростей находятся z-координаты сетки на полуцелом слое по времени:
п+1/2 _ п | Т" -п «+1/2 _ „ Т~п .„ п„ч
Zi,k — zi,k 2 hk' Zi+l/2,k ~ Zi+l/2,k "Т" 2 i+l/2,fe'
и по координатам центров ребер сетки на полуцелом слое вычисляются новые
С.п+1/2 „ / п+1/2 п+1/2 \
плоЩади ячеек 5^+1/2^+1/2 = 0.5^г+1/Н2г+1/2,к + ^+1/2,к+ 1) .
2. Первая (балансная) фаза КАБАРЕ: аппроксимация дивергентных форм законов сохранения (1) на ячейках сетки на слое п и нахождение консервативных переменных ^П+11//22к+1/2 на полуцелом слое (см. обозначения и шаблон на рис. 1):
(Sf)П+1/2 - (Sf)n
с 4- яп (-/п\ a"i7™
+ aMZ2 — Z1 J — aMZ3 ~~ Z4 J
Тп/2
— а£ (^ — ¿п) + аВ (^п — 4) + — ёВ = 0. (11)
Данные разностные уравнения записаны в векторной форме и используют следующие векторы:
f = (в, ви, вад)т, а = (ви, ви2 + ¿Р/р0, виад)т,
Ы, 22 )
(2:4,24)
Гв
(11,21)
Рис. 1. Шаблон консервативных фаз (11) и (15).
ё = (в('ш - г), ви(ад - г), - г) + бР/ро)1
3. Вторая (характеристическая) фаза КАБАРЕ: вычисление локальных инвариантов Римана по направлениям х и г, их экстраполяция и нахождение потоковых переменных ^П+1/2 к и ^"++1/2 на следующем целом слое по времени.
Для каждой ячейки сетки с помощью консервативных значений, полученных на первой фазе (11), вычисляются собственные значения (Ах):+1/2 и (А|)"+1/2 на полуцелом слое по времени по формулам (5). Затем на основании знаков полученных собственных значений определяются направления переноса локальных инвариантов Римана и проводится их экстраполяция на ребра сетки на слое п + 1. Так, например, экстраполяция по направлению х осуществляется следующим образом (см. шаблон и обозначения на рис. 2):
( ТХ)П + 1
(13 )г,к+1/2
хЧ:+1/2 ( тх\:
- (1з )ЬЬ,
2(1*)сь
Ч п+1/2
х \ : X )НН,
2(1*):г2 - ах) от ):+1/2+(/х):г/21
если 0.5 [(Ах если 0.5[(Ах
X
Л сЬ :+1/2 X ) сЬ
п+1/2 + (А^^21 > 0,
xv + (ах):+1/2 1
иначе
чин \<0 (12)
Рис. 2. Шаблон характеристической фазы для экстраполяции инвариантов по оси х (12).
Значения требуемых для экстраполяции локальных инвариантов при этом вычисляются по формулам (5), в которых величины с индексами 1ос берутся
в середине той пространственно-временной ячейки, со стороны которой осуществляется перенос. Получившиеся в результате (12) значения дополнительно подвергаются монотонизации на основе принципа максимума [12].
Экстраполяция по направлению г осуществляется аналогично. С помощью найденных значений инвариантов на ребрах сетки на слое п +1 находятся значения потоковых переменных ^Г+1% к и ^Г++1/2. Особого внимания требует алгоритм нахождения потоковых переменных на граничных ребрах сетки (как на границах-стенках, так и на свободной границе). В данной работе мы не будем останавливаться на вопросах реализации граничных условий, подробно это описано в [6].
4. Вычисление положения сетки и площадей ячеек на следующем целом слое по времени п + 1. Скорости узлов на свободной границе принимаются равными скоростям, использованными при предыдущем передвижении сетки (8), чтобы узлы дошли до положений, предписанных им экстраполяцией (7):
гп+1 = г"к. (13)
Скорости центров ребер на свободной границе находятся с помощью аппроксимации граничного условия (4) по потоковым значениям скоростей и, ад, найденным на характеристической фазе алгоритма (третья «балансная» фаза для свободной границы):
гГ+1 _ гг+1
¿п+1 1 _,,"+! г+1'° г'° (]Л)
2+1 /2,0 2+1 /2,0 г+1/2,0 +1/2 '
Окончательно с помощью вычисленных скоростей находятся г-координаты сетки на следующем целом слое по времени:
= + (15)
и по координатам центров ребер сетки на слое п+1 вычисляются новые площади ячеек 5,"+11/21к+1/2 = 0.5^г+1/2(г2"+1/2,к + гг"++11/2,к+1).
5. Третья (балансная) фаза КАБАРЕ: аппроксимация дивергентных форм законов сохранения (1) на ячейках сетки на слое п + 1 и нахождение консервативных переменных ^Г+1% к+1/2 на следующем слое (обозначения и шаблон совпадают с первой фазой (11)):
/сиЛг+1 /о р\Г+1/2 I*31 )с_I*31 )с__п+1/ п+1
^72 я ( з
+ аЯ+Чг2+1 - гГ'+1) - аь+Чг ' ~ - г
ЪГ^1 (г"+1 - г^1) + Ъ^1 (г"+1 - г^1) + ёГ^Ч - кс = 0. (16)
6. Вычисление следующего шага по времени тГ+1 по заданному числу Куранта СРЬ:
, гГ+1 _ г Г+1
"г+1/2 г+1 /2,/с г+1/2, к+1
гп+1 = т-тш 1ШП —^-г, Т /;п+1 ) ). (17)
г,к,А \|(АХ/2,к+1/2 I )!+1 /2,к+1/2 I
Предложенный алгоритм (6)—(17) представляет из себя две соединенных схемы КАБАРЕ: (11), (12), (16) — схемы КАБАРЕ для системы уравнений (1), и (6), (7), (14) — схемы КАБАРЕ для уравнения движения свободной границы (4). Таким образом, данный алгоритм обладает всеми типичными свойствами схем семейства КАБАРЕ: вторым порядком аппроксимации по пространству и времени, обратимостью по времени при отключенных процедурах монотонизации (за исключением, разве что, алгоритма выбора направления переноса на свободной границе (7)), минимальным вычислительным шаблоном в одну пространственно-временную ячейку, минимальной численной вязкостью.
4. Силы поверхностного натяжения
При моделировании неустойчивого перемещения свободной границы в поле гравитационных сил в первую очередь развиваются высокочастотные возмущения свободной границы, рост которых очень быстро приводит к аварийному останову программы. Для того чтобы изучить динамику неустойчивости самой низкой гармоники, надо ввести в модель силы поверхностного натяжения, которые бы не давали развиваться высокочастотным возмущениям.
Введем силы поверхностного натяжения сначала в дифференциальную модель (1). Энергия поверхностного натяжения пропорциональна площади свободной поверхности:
а
где ц € [а, Ь] — параметр в уравнениях свободной поверхности х = х(д), г = г(д), а — коэффициент поверхностного натяжения.
Компоненты вектора сил поверхностного натяжения можно получить из вариации энергии (18):
а
откуда следуют выражения для компонент удельных сил:
рх _ d ( f dx\ ^ f dz\ ' дхЛ
дч\ [\dqj \dqJ \ dq)'
Fz - -a—i\(—\2 ( — У]1/2— \
\d(lJ \ dqj'
-1/2,- - (19)
Определим теперь силы поверхностного натяжения для разностного случая. Рассмотрим следующий разностный аналог энергии поверхностного натяжения (18), в котором интеграл по свободной поверхности аппроксимируется с помощью координат центров ребер сетки:
А^-1
Е1епз = а ^ V(Жг+1/2 - Жг_1/2)2 + (^+1/2,0 - ^-1/2,0)2
+ (ж0+1/2 - Ж0)2 + (¿0+1/2,0 - zO,o)2
+ (TyJ(ждг^ - Жлг^-1/2)2 + {zNn,о - ZNx-1/2,О)2- (20)
Тогда удельные силы поверхностного натяжения, действующие перпендикулярно ребру с центром в точке (xi+1/2,zi+1/200\ определяются так:
IPX dEhenS / xi+3/2 — xi+1/2 xi+1/2 — xi-1/2 г a I 1 /о — - — —(7
г+1'2 дхм/о I l:+1 i
ci+1/2 \ 4+1
„z _ dEhenS _ / zi+3/2,0 — zi+1/2,0 zi+1/2,0 — zi-1/2,0 \ (91)
l+1/2 - 77--- - I -7---7- 1 ' V ;
dzi+1/2,0 \ 4+1 4
Ь — у (хг+1/2 — хг-1/2)2 + (^г+1/2,0 — ^г-1/2,0)2-
Отметим, что формулы (21) справедливы для ребер на свободной поверхности, не граничащих с левой и правой стенками. Формулы для ребер у этих границ можно получить, продифференцировав (20) по соответствующим переменным. Очевидно, разностные выражения (21) аппроксимируют дифференциальные операторы (19) со вторым порядком.
Введенные силы (21) добавляются в правые части балансных фаз алгоритма (11) и (16) для слоя ячеек, граничащих со свободной поверхностью. Так, например, первая фаза (11) модифицируется следующим образом:
(ер )п+1/2 _(СГ)п
V01 Л+1/2,0+1/2 V /г+1/2,0+1/2 „п {„п «л / п «Л --г аН\г2 ~ г1 ) ~ аЬ\г3 ~~ )
г/2
- bn (z2 - z3) + (zn - zn) + hi+1/2 - hi+1/2 = F!+1 /2^1+1/2,0 , (22)
^ — (0, РХ ^ * У , 1П+1/2 = V (*П+1 — *П) 2 + «+1,0 — <о)2.
Третья фаза (16) модифицируется аналогично.
Отметим, что учет сил поверхностного натяжения в уравнениях (22) накладывает дополнительные, вообще говоря, более жесткие условия на шаг по времени, чем (17). Чтобы их учесть, приведенные в следующем разделе расчеты проводились с относительно малым числом Куранта СЕЬ.
5. Тестовые расчеты
Рассмотрим следующую модельную задачу для уравнений (1), (2). Пусть слабосжимаемая жидкость находится в начальный момент времени в области (ж, г) € [0, х [0,Н0], длина области = 40, средняя высота жидкости Н0 = 20. Компоненты скорости и(ж,г,4 = 0) = ад(ж, г, 4 = 0) = 0, объемы не возмущены: 6(ж, г, 4 = 0) = 1.
Возмутим верхнюю (свободную) границу следующим образом:
Тогда при силе натяжения, направленной вниз (д > 0), свободная граница начнет устойчиво колебаться. А при силе натяжения, направленной вверх (д < 0), начнет развиваться неустойчивость по типу Рэлея — Тейлора [7]. Последнюю задачу можно трактовать по-разному. С одной стороны, она описывает процесс стекания закрепленной с одного конца плоской капли по наклонной крыше без учета сил трения. С другой стороны, она описывает движение жидкости вверх под воздействием некоторого магнитного поля [14,15]. В любом случае именно такая постановка задачи о развитии неустойчивости позволяет оценить робастность предложенного алгоритма в случае расчетов с сильно деформируемыми расчетными ячейками. Отметим, что мы рассматриваем только линейные стадии развития неустойчивости Рэлея — Тейлора, когда ускорение свободной поверхности экспоненциально стремится к |д| [16].
В этом разделе мы приводим результаты работы алгоритма (6)-(17) на серии задач с возмущенной границей (23) для е = 10-5 с различными значениями коэффициента поверхностного натяжения а и знаками ускорения свободного падения д. Все расчеты проводились на сетке 40 х 20 расчетных ячеек с равномерным шагом по ж и равномерным в каждом вертикальным столбце ячеек шагом по г. Число Куранта во всех расчетах полагалось достаточно малым и равным СРЬ = 0.05 для того, чтобы соблюсти условия устойчивости, связанные с введением в модель сил поверхностного натяжения (21). Как было показано в [6], алгоритм (6)-(17) дает хорошие качественные и количественные результаты при искусственной скорости звука с > ^\д\Но- Приведенные ниже расчеты проводились при с = 5^/\д\Но.
5.1. Развитие неустойчивости без сил поверхностного натяжения.
Рассмотрим задачу (23) с параметрами д = —1 и а = 0. Результаты расчетов по данной задаче на момент времени 4 = 6.864 (последний момент времени перед аварийным остановом программы) представлены на рис. 3. Как уже упоминалось ранее, при отсутствии сил поверхностного натяжения на свободной поверхности начинают расти высокочастотные возмущения, что очень быстро приводит к аварийному останову и не позволяет оценить ускорение узлов свободной поверхности. Данный расчет демонстрирует важность учета в модели сил поверхностного натяжения. Дополнительно отметим, что величины 66 на последний момент времени достигают значений, немного больших допустимого
(23)
для приближения слабой сжимаемости порога в 0.01. Тем не менее авост наступает не из-за превышения данного порога, а из-за слишком сильной деформации расчетных ячеек. Это замечание справедливо и для всех последующих тестов, заканчивающихся аварийным остановом.
□□ □□
□□ □□ □□
(а)
(б)
Рис. 3. Результаты расчетов задачи (23) с параметрами д = -1 и а = 0 (развитие неустойчивости без сил поверхностного натяжения) на момент времени 4 = 6.864. (а) профиль возмущения х-координаты свободной поверхности 8х = х — Но . (б) расчетная сетка на указанный момент времени и профиль 89 на ней.
/ПЬез С: -0.02 -0.01 0 0.01 0.02
(а)
(б)
Рис. 4. Результаты расчетов задачи (23) с параметрами д = —1 и а = 2.5 (развитие неустойчивости с недостаточными силами поверхностного натяжения) на момент времени 4 = 26.4882. (а) профиль возмущения х-координаты свободной поверхности 8х = х — Но . (б) расчетная сетка на указанный момент времени и профиль 89 на ней.
5.2. Развитие неустойчивости с недостаточными силами поверхностного натяжения. Рассмотрим задачу (23) с параметрами д = —1 и
а = 2.5. Результаты расчетов по данной задаче на момент времени 4 = 26.4882 (последний момент времени перед аварийным остановом программы) представлены на рис. 4. Приведенные результаты демонстрируют динамику свободной поверхности в случае, когда силы поверхностного натяжения компенсируют неустойчивое развитие не всех высокочастотных гармоник. При этом расчет может быть осуществлен до более далекого времени, чем в случае с а = 0, но нескомпенсированные высокие гармоники возмущений на достаточно далеком времени слишком сильно деформируют расчетные ячейки и «обрушивают» расчет.
Отметим, что при расчетах задач о развитии неустойчивости в какой-то момент могут получаться несимметричные результаты (что и происходит в приведенном расчете), несмотря на то, что исходная задача и сам алгоритм являются симметричными. Это объясняется возрастающей при неустойчивых расчетах ролью ошибок в операциях над числами с плавающей точкой. Все расчеты устойчивых задач обладают свойством симметрии, что было показано в работе [6] и будет продемонстрировано в некоторых последующих тестах.
5.3. Развитие неустойчивости с достаточными силами поверхностного натяжения. Рассмотрим задачу (23) с параметрами д = —1 и а = 3.2. Указанное значение а подобрано экспериментальным образом и является минимальным значением коэффициента поверхностного натяжения, при котором компенсируется неустойчивый рост всех высокочастотных возмущений свободной поверхности. Результаты расчетов по данной задаче на момент времени 4 = 29.3428 (последний момент времени перед аварийным остановом программы) представлены на рис. 5. Отметим, что некоторые вертикальные слои расчетных ячеек растянулись или сжались почти в 2 раза, прежде чем расчет был «обрушен» формированием слишком косых ячеек. Это говорит о высокой устойчивости представленного алгоритма (6)—(17) к сильному растяжению и сужению ячеек.
Неустойчивое развитие только самого низкочастотного возмущения позволяет исследовать ускорение свободной границы. Согласно теории ускорение границы по оси г должно стремиться к | д| , при этом на начальном этапе расти экспоненциально. Графики для ускорения среднего узла свободной поверхности г№х/2,0 и ¿п(г№х/2,0) представлены на рис. 6. График на рис. 6(а) позволяет заключить, что примерно на момент времени 4 = 2.8 усредненное значение ускорения максимально приближается к предельному значению, после чего, когда силы поверхностного натяжения становятся слишком велики, начинает уменьшаться. При этом значения ускорения испытывают колебания, увеличивающиеся в амплитуде с течением времени (из-за этого иногда наблюдаются значения ускорения, большие предельного |д|). Эти колебания связаны с конечностью скорости звука в предложенной математической модели (1) и минимальной численной вязкостью используемого алгоритма: любые ошибки по сравнению с дифференциальной моделью практически не затухают, а в силу неустойчивости только растут и со временем сказываются на движении узлов
(а) (б)
Рис. 5. Результаты расчетов задачи (23) с параметрами д = —1 и ст = 3.2 (развитие неустойчивости с достаточными силами поверхностного натяжения) на момент времени 4 = 29.3428. (а) Профиль возмущения ^-координаты свободной поверхности 8х = х — Но . (б) Расчетная сетка на указанный момент времени и профиль 89 на ней.
(а) (б)
Рис. 6. Результаты расчетов задачи (23) с параметрами д = —1 и ст = 3.2 (развитие неустойчивости с достаточными силами поверхностного натяжения). (а) Зависимость ускорения среднего узла свободной поверхности ¿N^/2,0 от времени. (б) Зависимость 1п(Х_^/2,0) от времени.
на свободной границе.
График на рис. 6(б) показывает, что в течение всего расчета, кроме начала и конца, ускорение увеличивалось экспоненциально. Нелинейность 1п(,г) в конечные момент времени связана с образованием слишком больших сил поверхностного натяжения, а в начальные моменты времени — с периодом «под-страивания» начального условия (23) под рассматриваемую область. В самом деле, возмущение в виде синуса не является элементарным решением уравнения динамики свободной границы рассматриваемой области, и алгоритму нуж-
(а)
(б)
Рис. 7. Колебания свободной поверхности в поле сил тяжести. (а) д = -1, а = 50. (б) д =1, а = 3.2. 0 — начальный профиль. 1 — половина периода первого колебания. 2 — полный период колебаний. 3 — половина периода второго колебания. 4 — два полных периода колебаний.
но некоторое количество шагов по времени, чтобы поверхностное натяжение подавило в синусе высокочастотные возмущения.
5.4. Подавление неустойчивости силами поверхностного натяжения. При достаточно больших значениях коэффициента поверхностного натяжения (а > <тсгй) свободная поверхность вместо экспоненциального роста начнет устойчиво колебаться. Пограничное значение <тсгй можно оценить из условия равенства энергии поверхностного натяжения и приращения потенциальной энергии в начальный момент времени.
Результаты расчетов по задаче (23) с параметрами д = -1 и а = 50 представлены на рис. 7(а). На графике представлены профили приращений ¿-координат узлов свободной поверхности на разные моменты времени, когда средний узел границы принимает самое высокое и самое низкое положения за период колебаний. Приведенные результаты позволяют заключить, что учет в алгоритме сил поверхностного натяжения практически не привносит в алгоритм дополнительной численной вязкости: поверхность продолжает колебаться с примерно одинаковой амплитудой (на графиках представлены только профили свободной поверхности на первые 2 периода колебаний, но заметного затухания не возникает и при дальнейших расчетах). Кроме того, наличие существенного поверхностного натяжения способствует более быстрой трансформации начального возмущения (синуса) в элементарное решение уравнения динамики свободной поверхности.
5.5. Поверхностное натяжение в поле гравитационных сил, направленных вниз. Продемонстрируем, что силы поверхностного натяжения (21) можно использовать и в типичных для океанологии задачах, когда силы тяжести направлены в сторону нижний границы области. На таких задачах с
X
а = 0 алгоритм (6)—(17) уже был протестирован в [6] при достаточно больших отклонениях свободной границы.
Рассмотрим задачу (23) об устойчивых колебаниях свободной поверхности с параметрами g = 1 и а = 3.2. Выбранный коэффициент поверхностного натяжения равен минимальному значению, при котором в аналогичной неустойчивой задаче не образуется высокочастотных неустойчивостей. Результаты расчетов по данной задаче на разные моменты времени, когда средний узел границы принимает самое высокое и самое низкое положения за период колебаний, представлены на рис. 7(б). Приведенные результаты хорошо согласуются с полученными ранее графиками для системы уравнений газовой динамики [5]. При этом минимумы и максимумы свободной границы на разных периодах колебаний могут отличаться более существенно, чем в предыдущем тесте. Это объясняется тем, что а в этой задаче относительно мало, из-за чего окрестности локальных экстремумов колеблются с отличной от всего остального профиля частотой. При увеличении а результаты будут все больше походить на рис. 7(а).
6. Заключение
В данной работе была предложена адаптация метода [6] на задачи динамики свободной поверхности в слабосжимаемом приближении с учетом сил поверхностного натяжения. Разработанный метод обладает вторым порядком аппроксимации по времени и пространству, минимальным вычислительным шаблоном в одну пространственно-временную ячейку и минимальной численной вязкостью. В модели были учтены силы поверхностного натяжения, которые позволили подавить высокочастотные колебания свободной поверхности и провести исследование линейной стадии развития неустойчивости Рэлея — Тейлора.
Проведенные тесты показали возможность моделирования течения слабо-сжимаемой жидкости балансно-характеристическими методами на расчетных сетках в смешанных эйлерово-лагранжевых координатах. При этом при проведении неустойчивых расчетов некоторые ячейки растягивались или сужались почти в 2 раза до «обрушения» расчета, что говорит о высокой робастности разработанного метода. В случае расчетов устойчивых колебаний свободной поверхности введенные силы поверхностного натяжения не привносили в модель существенной численной вязкости.
ЛИТЕРАТУРА
1. Karabasov S. A., Goloviznin V. M. Compact accurately boundary-adjusting high-resolution technique for fluid dynamics // J. Comput. Phys. 2009. V. 228. N 19. P. 7426-7451.
2. Головизнин В. М., Майоров Павел А., Майоров Петр А., Соловьев А. В. Новый численный алгоритм для уравнений многослойной мелкой воды на основе гиперболической декомпозиции и схемы КАБАРЕ // Морской гидрофиз. журн. 2019. Т. 35, № 6. С. 600620.
3. Goloviznin V. M., Maiorov Pavel A., Maiorov Petr A., Solovjev A. V. Validation of the low dissipation romputational algorithm CABARET-MFSH for multilayer hydrostatic flows with a free surface on the lock-release experiments //J. Comput. Phys. 2022. V. 463. 111239.
4. Головизнин В. М., Афанасьев Н. А. Бесшовный балансно-характеристический метод решения задач взаимодействия жидкости и газа с деформируемыми объектами // Мат. моделирование. 2021. T. 33, № 10. С. 65-82.
5. Афанасьев Н. А., Майоров П. А. Схема КАБАРЕ на подвижных сетках для двумерных уравнений газовой динамики и динамической упругости // Вычисл. методы и программирование. 2021. Т. 22, № 4. С. 306-321.
6. Головизнин В. М., Майоров Петр А., Майоров Павел А., Соловьев А. В., Афанасьев Н. А. Негидростатическая модель динамики жидкости CABARET-NH // Мат. моделирование. (Принята в печать).
7. Drazin P. G., Reid W. H. Hydrodynamic stability. Cambridge: Camb. Univ. Press, 2004.
8. Yilmaz I., Edis F. O., Saygin H. Application of an all-speed implicit finite-volume algorithm to Rayleigh-Taylor instability // Int. J. Comput. Methods. 2015. V. 12, N 3. 1550018.
9. Shadloo M. S., Zainali A., Yildiz M. Simulation of single mode Rayleigh-Taylor instability by SPH method // Comput. Mech. 2013. V. 51. P. 699-715.
10. Tryggvason G. Numerical simulations of the Rayleigh-Taylor instability //J. Comput. Phys. 1988. V. 75, N 2. P. 253-282.
11. He X., Chen S., Zhang R. A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability //J. Comput. Phys. 1999. V. 152, N 2. P. 642-663.
12. Головизнин В. М., Зайцев М. А., Карабасов С. А., Короткин И. А. Новые алгоритмы вычислительной гидродинамики для многопроцессорных вычислительных комплексов. М.: Изд-во Моск. ун-та, 2013.
13. Головизнин В. М., Карабасов С. А. Нелинейная коррекция схемы КАБАРЕ // Мат. моделирование. 1998. Т. 10, № 12. С. 107-123.
14. Гасилов В. А., Головизнин В. М. Численное решение одной модельной задачи о релей-тейлоровской неустойчивости // Препринт ИПМ АН СССР № 19. М.: ИПМ АН СССР, 1977. 23 c.
15. Hillier A., Berger T., Isobe H., Shibata K. Numerical simulations of the magnetic Rayleigh-Taylor instability in the Kippenhahn-Schluter prominence model. I. Formation of upflows // Astrophys. J. 2012. V. 746. P. 120.
16. Duchemin L., Josserand C., Clavin P. Asymptotic behavior of the Rayleigh-Taylor instability // Phys. Rev. Lett. 2005. V. 94, N 22. 224501.
Поступила в редакцию 29 октября 2022 г. После доработки 27 ноября 2022 г. Принята к публикации 29 ноября 2022 г.
Афанасьев Никита Александрович
Московский государственный университет им. М. В. Ломоносова, Ленинские горы, 1, Москва 119991 vmnaf @cs.msu.ru
Головизнин Василий Михайлович Московский государственный университет им Ленинские горы, 1, Москва 119991 gol@ibrae.ac.ru
Майоров Петр Александрович Московский государственный университет им Ленинские горы, 1, Москва 119991; ИБРАЭ РАН,
Б. Тульская ул., 52, Москва 115191 maiorov.peter@gmail.com
Соловьев Андрей Валерьевич ИБРАЭ РАН,
Б. Тульская ул., 52, Москва 115191 solovjev@ ibrae.ac.ru
. М. В. Ломоносова,
. М. В. Ломоносова,
Математические заметки СВФУ Октябрь—декабрь, 2022. Том 29, № 4
UDC 519.6
SIMULATING THE DYNAMICS OF A FLUID WITH
A FREE SURFACE IN A GRAVITATIONAL
FIELD BY A CABARET METHOD
N. A. Afanasiev, V. M. Goloviznin, P. A. Maiorov, and A. V. Soloviov
Abstract: An explicit conservative-characteristic CABARET method is proposed for calculating the dynamics of a fluid with a free surface in a gravitational field in a weakly compressible approximation. The developed method has the second order of approximation in time and space, minimum computational template of one space-time cell and minimum numerical viscosity. A difference scheme is tested on problems with various values of the surface tension coefficient and gravitational acceleration with various signs, including the problem of the development of the Rayleigh—Taylor instability. Taking into account the forces of surface tension makes it possible to get rid of high-frequency oscillations on the free surface when calculating unstable problems and regularizes the solution.
DOI: 10.25587/SVFU.2023.23.70.007
Keywords: equations of hyperbolic type, balance-characteristic schemes, mixed Euler— Lagrangian variables, weakly compressible fluid, Rayleigh—Taylor instability, surface tension.
REFERENCES
1. Karabasov S. A. and Goloviznin V. M., "Compact accurately boundary-adjusting high-resolution technique for fluid dynamics," J. Comput. Phys., 228, No. 19, 7426-7451 (2009).
2. Goloviznin V. M., Maiorov Pavel A., Maiorov Petr A., and Solovjev A. V., "New numerical algorithm for the multi-layer shallow water equations based on the hyperbolic decomposition and the CABARET scheme," Phys. Oceanography, 26, No. 6, 528-546 (2019).
3. Goloviznin V. M., Maiorov Pavel A., Maiorov Petr A., and Solovjev A. V., "Validation of the low dissipation computational algorithm CABARET-MFSH for multilayer hydrostatic flows with a free surface on the lock-release experiments," J. Comput. Phys., 463, Article ID 111239 (2022).
4. Goloviznin V. M. and Afanasiev N. A., "Monolithic balance-characteristic method for solving problems of the interaction of a liquid and gas with deformable objects," Math. Models Comput. Simulations, 14, No. 3, 398-410 (2022).
5. Afanasiev N. A. and Maiorov P. A., "CABARET scheme on moving grids for two-dimensional equations of gas dynamics and dynamic elasticity [in Russian]," Vychisl. Metody i Program., 22, No. 4, 306-321 (2021).
6. Goloviznin V. M., Maiorov Pavel A., Maiorov Petr. A., Solovjev A. V., and Afanasiev N. A., "Nonhydrostatic fluid dynamics model CABARET-NH," Mat. Model. (2023) [accepted].
7. Drazin P. G. and Reid W. H., Hydrodynamic Stability, Camb. Univ. Press, Cambridge (2004).
8. Yilmaz I., Edis F. O., and Saygin H., "Application of an all-speed implicit finite-volume algorithm to Rayleigh-Taylor instability," Int. J. Comput. Methods, 12, No. 3, Article ID 1550018 (2015).
© 2022 N. A. Afanasiev, V. M. Goloviznin, P. A. Maiorov, and A. V. Soloviov
9. Shadloo M. S., Zainali A., and Yildiz M., "Simulation of single mode Rayleigh-Taylor instability by SPH method," Comput. Mech., 51, 699-715 (2013).
10. Tryggvason G., "Numerical simulations of the Rayleigh-Taylor instability," J. Comput. Phys., 75, No. 2, 253-282 (1988).
11. He X., Chen S., and Zhang R., "A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability," J. Comput. Phys., 152, No. 2, 642-663 (1999).
12. Goloviznin V. M., Zaitsev M. A., Karabasov S. A., and Korotkin I. A., Novel Algorithms of Computational Hydrodynamics for Multicore Computing [in Russian], Mosk. Gos. Univ., Moscow (2013).
13. Goloviznin V. M. and Karabasov S. A., "Nonlinear correction of Cabaret scheme [in Russian]," Mat. Model., 10, No. 12, 107-123 (1998).
14. Gasilov V. A. and Goloviznin V. M., "Numerical solution of a model problem of Rayleigh-Taylor instability [in Russian]," Preprint IPM AN SSSR, No. 19, IPM AN SSSR, Moscow (1977), 23 p.
15. Hillier A., Berger T., Isobe H., and Shibata K., "Numerical simulations of the magnetic Rayleigh-Taylor instability in the Kippenhahn-Schluter prominence model. I. Formation of upflows," Astrophys. J., 746, 120 (2012).
16. Duchemin L., Josserand C., and Clavin P., "Asymptotic behavior of the Rayleigh-Taylor instability," Phys. Rev. Lett., 94, No. 22, 224501 (2005).
Submitted October '27, 2022 Revised November 27, 2022 Accepted November 29, 2022
Nikita A. Afanasiev Lomonosov Moscow State University, 1 Leninskie Gory, Moscow 119991, Russia vmnaf@cs.msu.ru
Vasily M. Goloviznin Lomonosov Moscow State University, 1 Leninskie Gory, Moscow 119991, Russia gol@ibrae.ac.ru
Petr A. Maiorov
Lomonosov Moscow State University, 1 Leninskie Gory, Moscow 119991, Russia; Nuclear Safety Institute of RAS, 52 B. Tulskaya Street, Moscow 115191, Russia maiorov.peter@gmail.com
Andrey V. Soloviov
Nuclear Safety Institute of RAS,
52 B. Tulskaya Street, Moscow 115191, Russia
solovjev@ibrae.ac.ru
