ЖЭТФ, 2020, том 157, вып. 5, стр. 944-956

ВОЛНЫ НАД ИСКРИВЛЕННЫМ ДНОМ:

МЕТОД СОСТАВНОГО КОНФОРМНОГО ОТОБРАЖЕНИЯ

В. П. Рубан^*

Институт теоретической физики им. Л. Д. Ландау Российской академии наук

142432, Черноголовка, Московская обл., Россия

Поступила в редакцию 19 декабря 2019 г.,

после переработки 19 декабря 2019 г.

Принята к публикации 23 декабря 2019 г.

Описан компактный и эффективный численный метод исследования плоских течений идеальной жид-

кости с гладкой свободной границей над искривленным и неоднородно подвижным дном. Используются

точные уравнения движения в терминах так называемых конформных переменных. В дополнение к ра-

нее известному применению для сдвиговых течений с постоянной завихренностью (включая нулевую),

здесь сделано обобщение на случай потенциальных течений в равномерно вращающихся системах ко-

ординат, где к силе тяжести добавляются центробежная сила и сила Кориолиса. Дан краткий обзор

предыдущих результатов использования данного подхода в ряде физически интересных задач, таких

как моделирование волн цунами, обусловленных подвижкой неоднородного дна, динамика брэгговских

(щелевых) солитонов над пространственно-периодическим профилем дна, явление возврата Ферми - Па-

ста - Улама для волн в конечном бассейне, образование аномальных волн на встречном неоднородном

течении, распространение уединенной волны на сдвиговом течении и ее набегание на перепад глубин.

Кроме того, представлен ряд новых численных результатов, относящихся к нелинейной динамике свобод-

ной границы в частично заполненных жидкостью замкнутых вращающихся контейнерах — центрифугах

сложной формы. Уравнения движения в этом случае отличаются некоторыми существенными деталями

от x-периодических систем.

DOI: 10.31857/S0044451020050168

о них мало кто знал. Только с развитием и массо-

вым распространением вычислительной техники в

1. ВВЕДЕНИЕ

90-е гг. прошлого века и с эффективным практи-

В некоторых разделах гидродинамики прибли-

ческим воплощением алгоритмов быстрого дискрет-

жение идеальной несжимаемой жидкости оказыва-

ного преобразования Фурье, уравнения движения

ется достаточно приемлемым. В частности, сюда от-

свободной поверхности в конформных переменных

носится теория поверхностных волн на воде. Описа-

были заново открыты Захаровым и его сотрудни-

ние течений существенно упрощается, если считать

ками [3-6], став затем рабочим инструментом для

их потенциальными. К тому же, поскольку в про-

многих исследователей (см. [7-34] и многочислен-

странственно-двумерном случае уравнение Лапласа

ные ссылки там). Главным образом, моделируются

для потенциала поля скорости конформно инвари-

волны на бесконечно глубокой воде, для которых за

антно, а речь идет об областях с кривыми граница-

два десятка лет получены важнейшие результаты,

ми, на помощь приходит теория конформных отоб-

включая пионерские наблюдения аномальных волн

ражений и соответствующих аналитических функ-

(так называемых волн-убийц) в высокоточных чис-

ций. Поэтому неудивительно, что идея использо-

ленных экспериментах группы Захарова [7-10]. Но

вать конформные отображения для точного пред-

существует ряд интересных задач, в которых прин-

ставления уравнений нелинейной динамики потен-

ципиально важными являются эффекты взаимодей-

циальных поверхностных волн возникла достаточ-

ствия волн с неоднородным профилем дна. С це-

но давно (в работах Овсянникова [1,2]). Долгое вре-

лью их исследования метод конформных перемен-

мя эти уравнения оставались невостребованными, и

ных был обобщен в 2004 г. автором данной статьи на

случай искривленного неподвижного дна [35]. Вол-

^* E-mail: ruban@itp.ac.ru

944

ЖЭТФ, том 157, вып. 5, 2020

Волны над искривленным дном. . .

ны над подвижным дном были затем рассмотрены

ускорениям, а также деформациям формы попереч-

в работе [36], а сдвиговые течения со свободной гра-

ного сечения. При достаточно больших скоростях

ницей над переменной глубиной были промоделиро-

вращения Ω ≫ ν/h² (где ν — кинематическая вяз-

ваны методом конформного отображения в работе

кость) вязкие эффекты становятся несущественны-

[37]. Динамическими объектами в этом описании яв-

ми на временах порядка нескольких десятков Ω-1,

ляются потенциальная составляющая поля скорости

что соответствует многим оборотам. Весьма важно,

и конформное отображение горизонтальной полосы

что течение жидкости при этом вовсе не потенци-

0 ≤ v ≤ 1 из вспомогательной комплексной плоскос-

ально, а близко к твердотельному вращению с угло-

ти w = u + iv на подвижную область в вертикаль-

вой скоростью Ω. Поэтому естественным шагом ока-

ной плоскости z = x + iy, занятую жидкостью. Кон-

зывается переход во вращающуюся систему коорди-

формное отображение z(w, t) = Z(ζ(w, t), t) предста-

нат, где, как известно, к силе тяжести добавляются

вимо в виде композиции двух аналитических функ-

центробежная сила и сила Кориолиса. Возмущения

ций — неизвестной ζ(w, t) и заданной Z(ζ, t), при-

(бездивергентного) поля скорости под воздействием

чем первая из них оставляет действительную ось на

ускорений и деформаций контейнера оказываются

месте (т. е. Im ζ(u, t) = 0), а вторая задает парамет-

чисто потенциальными, так что имеет место обоб-

рически форму дна X(b)(s, t) + iY(b)(s, t) = Z(s, t).

щенное уравнение Бернулли, в котором отдельным

В результате для параметризации формы свобод-

слагаемым фигурирует гармонически сопряженная

ной границы достаточно одной неизвестной веще-

к потенциалу функция тока (из-за силы Кориоли-

ственной функции a(u, t) = Re ζ(u + i, t). Комплекс-

са). Надо сказать, что даже относительно более про-

ный потенциал скорости с учетом кинематического

стая задача о нелинейных волнах в круглом соосном

условия на нижней границе также зависит от од-

контейнере и при отсутствии силы тяжести иссле-

ной неизвестной вещественной функции ψ(u, t). Из

дована мало, хотя является весьма изящной и со-

обобщенного уравнения Бернулли и кинематическо-

держательной, не говоря уже о полном варианте с

го условия на верхней (свободной) границе следуют

некруглыми деформируемыми контейнерами и их

уравнения движения для a(u, t) и ψ(u, t). Уравне-

ускорениями. Поэтому применение конформных пе-

ния содержат линейные операторы, диагональные в

ременных здесь выглядит очень перспективным.

фурье-представлении, так что псевдоспектральный

Дальнейшее изложение организовано следую-

численный метод с использованием быстрого пре-

щим образом. В разд. 2 дается достаточно подроб-

образования Фурье оказывается наиболее естествен-

ное описание используемого метода, включая вывод

ным и удобным.

точных уравнений движения и обсуждение числен-

За прошедшее с тех пор время был получен еще

ной схемы. В разд. 3 кратко обсуждаются наибо-

ряд результатов, которые подтвердили практичес-

лее интересные предыдущие результаты, получен-

кую пользу и высокую эффективность метода. В

ные автором путем использования составного кон-

частности, были исследованы: динамика брэгговс-

формного отображения. Раздел 4 посвящен новому

ких (щелевых) солитонов над пространственно-пе-

применению конформных переменных для описания

риодическим профилем дна [38,39]; явление возвра-

поверхностных волн во вращающихся системах. На-

та Ферми - Паста - Улама для волн в конечном бас-

конец, разд. 5 содержит краткое заключение.

сейне [40, 41]; формирование аномальных волн на

2. ОПИСАНИЕ МЕТОДА

встречном неоднородном течении [42]. Цель пред-

лагаемой работы — сделать краткий обзор на эту

2.1. Точные уравнения движения

тему, а также представить некоторые новые резуль-

таты для ранее не обсуждавшейся в связи с кон-

Вывод точных и явных уравнений движения сво-

формными переменными классической задачи о ди-

бодной границы проделаем для случая волн в цен-

намике свободной поверхности в частично заполнен-

трифугах, поскольку бесконечные вдоль x системы

ных замкнутых вращающихся контейнерах (см., на-

во многом аналогичны и были подробно рассмотре-

пример, [43-45] и ссылки там). Практически речь

ны в цитированных выше работах [35-37].

может идти о квазидвумерных центрифугах (непра-

Как известно, уравнение Эйлера для плоских те-

вильной формы цилиндрах с характерным попереч-

чений идеальной несжимаемой жидкости во враща-

ным размером R и длиной h ≪ R, либо о круг-

ющейся (против часовой стрелки) системе коорди-

лых эксцентрических), которые вращаются вокруг

нат содержит дополнительные силы (покомпонент-

горизонтальной оси (ось параллельна образующим)

но) (Ω²x, Ω²y) + (2ΩV(y), -2ΩV(x)) и допускает ре-

и возможно подвержены угловому и/или линейному

дукцию V = (V(x), V(y)) = (ϕ_x, ϕ_y) = (θ_y , -θ_x), где

945

ЖЭТФ, вып. 5

В. П. Рубан

ЖЭТФ, том 157, вып. 5, 2020

нижние индексы обозначают частные производные,

u и u + i, они связаны между собой линейным пре-

функция ϕ(x, y, t) представляет собой потенциал по-

образованием (см. [36])

ля скорости, а θ(x, y, t) — соответствующую гармо-

нически сопряженную функцию тока. Обобщенное

Φ(s)(u, t) = e^-kΦ(b)(u, t),

(5)

уравнение Эйлера поэтому имеет вид

где e^-k

≡ exp(

∂_u). Это означает соотношение

(s)

Φ_k

(t) = e^-kΦ(b)k(t) для соответствующих фурье-об-

ϕ_t + (ϕ2x + ϕ2y)/2 - Ω²(x² + y²)/2 + 2Ωθ +

разов. Для дальнейшего изложения нам потребуют-

+ g(t)(y cosΩt + xsinΩt) +

P = 0,

(1)

ся еще несколько линейных операторов. Эти опера-

торы

^R, и

T =

R-1 также диагональны в фу-

где

P — давление, деленное на плотность жидкости.

рье-представлении:

Эффективное поле тяжести g(t) зависит от времени,

если ось вращения испытывает вертикальные уско-

S_k = 1/ch(k), R_k = i th(k), T_k = -i cth(k).

(6)

рения.

Важно, что в результате действия любого из этих

Пусть область течения в каждый момент време-

операторов на действительную функцию получает-

ни представляет собой деформированный круг с од-

ся также чисто действительная функция.

ной «дырой» — свободной поверхностью. Рассмот-

Уравнение (1), переписанное в новых перемен-

рим скалярную функцию v(x, y, t), удовлетворяю-

ных, на свободной границе имеет вид (динамическое

щую уравнению Лапласа v_xx + v_yy = 0 и принимаю-

граничное условие)

щую фиксированные граничные значения: v = 0

на внешней границе области (т. е. на стенке контей-

(

)

нера — дне), и v = 1 на свободной (внутренней)

Re Φ(s)t - Φ(s)uZ(s)t/Z(s)

+ |Φ(s)u/Z(s)u|²/2 -

границе. Такая функция существует и единствен-

− Ω²|Z(s)|²/2 + 2 Ω ImΦ(s) +

на. Построим теперь гармонически сопряженную к

+ g(t)Im[Z(s) exp(iΩt)] = 0.

(7)

ней функцию u(x, y, t). Эта функция многозначна

(и определена с точностью до аддитивной констан-

К нему добавляются два кинематических гранич-

ты). Ее приращение L(t) при обходе вдоль свободной

ных условия:

поверхности (конформный модуль) зависит, вообще

(

)

говоря, от времени. Составим комплексную комби-

Z(s)t

Z(s)u

= -ImΦ(s)u,

(8)

нацию w = u + iv, которая является аналитической

(

)

функцией комплексной переменной z = x + iy. Яс-

Z(b)t

Z(b)u

= -ImΦ(b)u,

(9)

но, что такая криволинейная замена координат со-

где

Z означает комплексно-сопряженную величину.

ответствует конформному отображению z(w, t) пря-

Для нас будет удобно представить Φ(b)(u, t) в виде

моугольника с размерами L × 1, причем переменная

v является аналогом радиальной координаты, а пе-

Φ(b)(u, t) =

Sψ(u, t) - i(1 - iR

∂-1uf(u, t),

(10)

ременная u — аналогом угловой координаты.

где ψ(u, t) и f(u, t) — некоторые неизвестные дейст-

Форма свободной границы теперь будет задана

вительные функции. Тогда, благодаря соотношению

параметрически формулой

(5) и операторным равенствам

X(s)(u, t) + iY(s)(u, t) ≡ Z(s)(u, t) = z(u + i, t),

(2)

e^-kS = (1 + i R), e-k(1 - i R) = ˆ,

(11)

тогда как профиль дна — формулой

для Φ(s)(u, t) мы будем иметь формулу

X(b)(u, t) + iY(b)(u, t) ≡ Z(b)(u, t) = z(u, t).

(3)

Φ(s)(u, t) = (1 + iR)ψ(u, t) -

∂-1uf(u, t).

(12)

Весьма важно, что комплексный потенциал

Далее мы используем в уравнениях функцию

ϕ(u, v, t) + iθ(u, v, t) = φ(w, t) представляет собой

аналитическую функцию. Обозначим граничные

Ψ(u, t) ≡ (1 + iR)ψ(u, t).

(13)

значения этой функции следующим образом:

Теперь необходимо принять во внимание, что

φ(u + i, t) ≡ Φ(s)(u, t), φ(u, t) ≡ Φ(b)(u, t).

(4)

функция z(w, t) может быть представлена как ком-

позиция двух функций (см. [35, 36]), т. е. z(w, t) =

Поскольку Φ(s)(u, t) и Φ(b)(u, t) являются значения-

= Z(ζ(w,t),t), где известная аналитическая функ-

ми одной и той же аналитической функции в точках

ция Z(ζ, t) = X(ζ, t) + iY (ζ, t) определяет форму

946

ЖЭТФ, том 157, вып. 5, 2020

Волны над искривленным дном. . .

дна отображением действительной оси. Конформ-

Далее рассуждаем так: поскольку

ное отображение Z(ζ, t) не должно иметь особенно-

стей в пределах достаточно широкой горизонталь-

ξ_t/ξ_u = ζ_t(w, t)/ζ_w(w, t)|w=u+i

ной полосы над действительной осью в плоскости

представляет собой взятое на u + i значение ана-

ζ, чтобы было «достаточно места» для волн боль-

шой амплитуды в плоскости z (в частном случае

литической фукции, которая действительна на дей-

ствительной оси, то имеет место соотношение между

у Z(ζ,t) могут отсутствовать особенности во всей

верхней полуплоскости ζ). Разложение конформно-

действительными и мнимыми частями: Im(ξ_t/ξ_u) =

= R Re(ξ_t/ξ_u), так что Im(ξ_t/ξ_u) = -Q означает ξ_t =

го отображения z(w, t) в композицию не единствен-

но, поскольку если Z(ζ, t) = Z₁(β(ζ, t), t), где анали-

= -ξ_u

T + i)Q. Это дает нам уравнение, определя-

тическая функция β(ζ, t) оставляет действительную

ющее a_t:

ось на месте и не имеет особенностей слишком близ-

a_t = - Re[ξ_u

T + i)Q].

(22)

ко от нее, то

Теперь мы подставляем нужные выражения в урав-

нение (7), для того чтобы найти ψ_t. Простые преоб-

Z(ζ(w, t), t) = Z₁(β(ζ(w, t), t), t) = Z₁(ζ₁(w, t), t).

разования приводят к следующему уравнению:

Промежуточная аналитическая функция ζ(w, t)

принимает действительные значения на действи-

⁽Ψ_uZ_t(ξ, t))-

ψ_t = - Re[Ψ_u

T + i)Q] + Re

тельной оси, и поэтому справедливо равенство

Z_ξ(ξ, t)ξ_u

∫

|Ψ_u|² -

Sf)²

a_k(t)

- g(t)Im[Z(ξ, t)exp(iΩt)] +

ζ(w, t) =

e^ikw

a_-k = a_k,

(14)

2|Z_ξ(ξ, t)ξ_u|²

ch(k)

2π

Ω²

|Z(ξ, t)|² - 2Ω(Im Ψ -

∂-1uf).

(23)

где a_k(t) есть фурье-преобразование некоторой дей-

ствительной функции a(u, t). На дне ζ(u, t)

Sa(u, t) и потому

Таким образом, выведены точные уравнения движе-

ния для волн на свободной поверхности идеальной

Z(b)(u, t) = Z

Sa(u, t), t).

(15)

жидкости в центрифугах. Что касается сдвиговых

течений в бесконечных вдоль x системах, то для

На свободной поверхности мы имеем соотношения

них все делается аналогично, но используется дру-

гое обобщенное уравнение Бернулли (см. [37]).

ζ(u + i, t) ≡ ξ(u, t) = (1 + iR)a(u, t),

(16)

Z(s) = Z(ξ, t), Z(s)u = Z_ξ(ξ, t)ξ_u,

(17)

2.2. Численное моделирование

Z(s)t = Z_ξ(ξ, t)ξ_t + Z_t(ξ, t).

(18)

Для практического применения системы (19),

(21), (22) и (23) в численном моделировании следует

Поскольку нашей целью является вывести урав-

принять во внимание два существенных обстоятель-

нения движения для неизвестных функций ψ(u, t),

ства, обсуждаемые ниже.

f (u, t) и ξ(u, t), нам следует подставить выражения

Во-первых, период L по переменной u (ранее упо-

(10), (12), (15), (16), (17) и (18) в уравнения (7), (8)

минавшийся конформный модуль отображения) за-

и (9). Уравнение (9) не требует никаких усилий и

висит от времени в силу сингулярности оператора

сразу же дает явное соотношение

T ≈ 1/(ik) =

∂-1u на малых k, что приводит к нали-

(

)



чию непериодической по u составляющей -a_u〈Q〉u в

f =Im Z_t(s,t

Z_s(s, t

Sa_u

a = Reξ.

(19)

правой части уравнения (22) и аналогичной состав-

〈Q〉u в правой части уравнения (23),

ляющей -ψ_u

Уравнение (8), деленное на |Zus)|², принимает вид

где 〈Q〉 — среднее значение величины Q по периоду

)

(

)

L. Поэтому вместо переменной u более удобно ис-

⁽ξ_t

Z_t(ξ, t)

- Im Ψ_u +

+ Im

(20)

пользовать новую переменную ϑ = [2π/L(t)]u, для

ξ_u

Z_ξ(ξ, t)ξ_u

|Z_ξ(ξ, t)ξ_u|²

которой период фиксирован и равен 2π. При этом

уравнение движения для величины α(t) = 2π/L(t)

Таким образом, мы имеем Im(ξ_t/ξ_u) = -Q, где

получается исходя из требования, чтобы при под-

(

)

Im Ψ_u -

Z_t(ξ, t)

становке в уравнения (22) и (23) непериодические по

Q≡

+ Im

(21)

|Z_ξ(ξ, t)ξ_u|²

Z_ξ(ξ, t)ξ_u

ϑ члены αua_ϑ и αuψ_ϑ в левых частях сокращались

947

12*

В. П. Рубан

ЖЭТФ, том 157, вып. 5, 2020

{

с соответствующими непериодическими составляю-

∂

⁽Ψ_ϑZ_t(ξ, t)⁾

ψ_ϑt =

- Re[Ψ_ϑ(Tα + i)Q] -

щими в правых частях. Очевидно, что сокращение

∂ϑ

Z_ξ(ξ, t)ξ_ϑ

происходит при α = -α〈Q〉.

|Ψ_ϑ|² - (^ŜF )²

Во-вторых, если среднее по периоду значение ве-

−

- g(t)Im[Z(ξ, t)exp(iΩt)] +

2|Z_ξ(ξ, t)ξ_ϑ|²

личины f отлично от нуля (что неизбежно имеет

}

место при деформации стенок центрифуги с изме-

Ω²

|Z(ξ, t)|²

- 2Ω(ImΨ_ϑ - ^ŜF),

(29)

нением площади A_c ее поперечного сечения), то в

системе появляется дополнительное потенциальное

круговое течение с циркуляцией Γ(t), так что 〈ψ_ϑ〉 =

где ξ = ϑ + iα + (1 + iRα)ρ, Ψ_ϑ = (1 + iRα)ψ_ϑ,

= Γ/2π, а сам потенциал оказывается многознач-

(

)

ной функцией. Полезно иметь в виду, что в силу



F =Im Z_t(s,t

Z_s(s, t)(1 +^Ŝρ_ϑ)

(30)

сохранения циркуляции скорости вдоль свободной

s=ϑ+Ŝρ

(

)

границы в лабораторной системе координат, имеет

Im Ψ_ϑ -^ŜF

Z_t(ξ, t)

+ Im

(31)

место интеграл движения Γ(t) + 2ΩA_c(t) = const.

|Z_ξ(ξ, t)ξ_ϑ|²

Z_ξ(ξ, t)ξ_ϑ

Данный закон сохранения следует обязательно про-

верять при отладке численного кода.

Эти новые операторы диагональны в дискретном

фурье-представлении:

Таким образом, можно отделить в a(ϑ, t) фикси-

рованную непериодическую часть, написав a(ϑ, t) =

R_α(m) = i th(αm), S_α(m) = 1/ ch(αm),

= ϑ + ρ(ϑ,t), где ρ(ϑ,t) — 2π-периодическая функ-

ция:

T_α(m) = -i cth(αm)

∑

ρ(ϑ, t) =

ρ_m(t)exp(imϑ).

(24)

для m = 0, тогда как T_α(0) = 0. Заметим, что опера-

m=-∞

ˆα не имеет особенности. Система замыкается

тор

следующим условием для α(t), которое обеспечива-

Аналогично для производной ψ_ϑ имеем

ет сокращение непериодических членов в уравнени-

ях (22) и (23):

∑

ψ_ϑ(ϑ, t) =

D_m(t)exp(imϑ),

(25)

∫

m=-∞

α(t) = -

Q(ϑ) dϑ.

(32)

2π

причем D₀ = Γ/2π, ρ_-m = ρ_-m, D_-m =

D_-m. При

этом формулы для соответствующих аналитических

У данной системы уравнений имеется интеграл

функций на свободной границе имеют вид

движения, соответствующий сохранению площади,

занятой жидкостью (площадь центрифуги A_c минус

∑

2ρ_m(t)exp(imϑ)

ξ(ϑ, t) = ϑ + iα(t) +

(26)

площадь полости A_h),

1 + exp(2mα(t))

m=-∞

∮ (

)

X(b)dY(b) - Y(b)dX(b)

∑

∮ (

)

2D_m(t)exp(imϑ)

Ψ_ϑ(ϑ, t) =

(27)

X(s)dY(s) - Y(s)dX(s)

(33)

1 + exp(2mα(t))

m=-∞

Кроме того, если контейнер вращается равномерно

С учетом сделанных выше замечаний, уравнения

с некоторой угловой скоростью Ω + Δ без измене-

движения для функций ρ и ψ_ϑ в переменных (ϑ, t)

ния своей формы, т. е. Z(ζ, t) = exp(iΔt)Z(ζ), то

выглядят аналогично уравнениям (22) и (23), но все

при g = 0 сохраняется сумма кинетической и цент-

u-производные должны быть заменены ϑ-производ-

робежной энергий в соответствующей системе коор-

ными, а прежние операторы

^R,

S, и

T должны быть

динат (завихренность жидкости в лабораторной си-

всюду заменены новыми операторами соответствен-

стеме при этом по-прежнему предполагается равной

но

R_α,

Ŝ_α, и

T_α:

2Ω). Если обозначить X+iY = Z(s) exp(-iΔt) и X_b+

= Z(b) exp(-iΔt), то сохраняющаяся энергия

+ iY_b

ρ_t = - Re[ξ_ϑ(Tα + i)Q],

(28)

есть

948

ЖЭТФ, том 157, вып. 5, 2020

Волны над искривленным дном. . .

∮

штабные волновые структуры на свободной поверх-

E_Δ =

(Ω²+2ΩΔ) (X²+Y²)(X dY - Y dX) +

ности. Как результат такого адаптивного увеличе-

∮

ния N, правые части эволюционных уравнений мо-

(X2b + Y2b)Rα(X2b + Y2b)_ϑ dϑ +

гут быть вычислены с примерно одинаковой числен-

∮

ной ошибкой δ₀ < N10-18, соответствующей типу

[X²+Y²]ψ_ϑ dϑ-

[X2b+Y2b]^Ŝψ_ϑ dϑ -

complex. Поскольку шаг интегрирования по време-

∮ (

)

Γϑ

Γ²α

ни уменьшается при этом как τ ∼ 1/N, погрешность

ψ-

R_αψ_ϑ dϑ +

(34)

в вычислении позиции свободной поверхности при

2π

4π

t ∼ 1 можно оценить как δ_s ≲ N²10-18. Практичес-

Если же Δ = -Ω и g = const (контейнер неподвижен

ки, интегралы движения сохраняются с точностью

в лабораторной системе координат), то сохраняется

до 10-12 десятичных знаков на протяжении большей

сумма кинетической энергии и потенциальной энер-

части эволюции. На конечной стадии, чем большее

гии в поле тяжести:

число N_final используется, тем позднее наступает

∮

момент, когда высокая точность теряется.

E-Ω +

Y²dX = const.

В конце этого раздела нужно сказать несколь-

ко слов о том, как можно задать начальную кон-

При тестировании компьютерной программы все

фигурацию свободной поверхности в конформных

эти законы сохранения также должны проверяться.

переменных, если известна ее зависимость, напри-

Иначе, как показала практика вычислений, можно

мер, в полярных координатах (r, χ) в виде явного

допустить такую ошибку, которая проявляется не во

выражения r = r(χ). Для этого следует положить

всех случаях и поэтому может остаться незамечен-

ρ = 0 и подобрать значение α₀ так, чтобы объем по-

ной.

лости соответствовал функции r(χ). Затем нужно

Уравнения (28)-(32) очень удобны и легки для

формальным образом временно убрать из обобщен-

численного моделирования, если функция Z(ζ, t) да-

ного уравнения Бернулли силу Кориолиса и видоиз-

ется достаточно компактной явной формулой, как

менить при этом потенциал центробежной силы:

это имеет место для многих практически интерес-

Ω²|Z(s)|²/2 → Ω²(|Z(s)|² - r²(χ))/2.

ных профилей дна. Кроме того, в языке програм-

мирования C имеется тип данных complex, а в ма-

Кроме того, следует положить g = 0 и добавить

тематической библиотеке присутствуют элементар-

в уравнение Бернулли линейное однородное зату-

ные функции комплексной переменной, такие как

хание, заменив там ψ_t → (ψ_t + γψ). В результате

cexp, clog и др., что делает процесс написания ко-

эволюции такой видоизмененной системы конфигу-

да довольно простым и приятным делом. Создан-

рация свободной границы быстро отрелаксирует к

ная автором численная схема естественным обра-

нужной форме, после чего можно будет вернуться к

зом опирается на дискретное преобразование Фу-

исходным уравнениям и начать счет.

рье, поскольку все линейные операторы в урав-

нениях могут быть эффективно вычислены в m-

представлении с помощью современного программ-

3. КРАТКИЙ ОБЗОР ПРЕДЫДУЩИХ

ного обеспечения (фактически используется библио-

РЕЗУЛЬТАТОВ

тека FFTW [46]), тогда как все нелинейные опера-

ции локальны в ϑ-представлении. В качестве основ-

Первые численные примеры, приведенные в ра-

ных динамических переменных берутся величины

ботах [35-37], ставили своей главной целью, скорее,

α(t), ρ_m(t) и D_m(t), с 0 ≤ m < M (для отрица-

продемонстрировать применимость метода к опи-

тельных m учитываются соотношения ρ_-m = ρ_m

санию динамики очень крутых (и даже опрокиды-

и D_-m =

D_m). С целью обеспечения устойчиво-

вающихся) волновых профилей над неоднородным

сти численной схемы при больших m, после каждой

дном, нежели проверить какие-либо теоретические

процедуры Рунге - Кутта четвертого порядка остав-

предсказания. Типичная постановка численных экс-

ляются только спектральные компоненты с |m| <

периментов — это задать начальную уединенную

< M_eff, где M_eff ≈ (1/4)N, M ≈ (3/8)N, а N =

волну над относительно глубоким дном (или создать

= 212...19 есть размер массивов для быстрого фу-

небольшую группу волн в результате подвижки дна,

рье-преобразования. В течение каждого численного

моделируя тем самым явление цунами), а затем на-

эксперимента число N удваивается несколько раз

блюдать за ее эволюцией при набегании на более

по мере надобности, если формируются мелкомас-

мелкий участок. Как и следовало ожидать исходя

949

В. П. Рубан

ЖЭТФ, том 157, вып. 5, 2020

из повседневного опыта, при выходе на мелково-

зывается почти полностью сосредоточенной в исход-

дье волна увеличивала свою амплитуду и опроки-

ной гармонике. Применительно к волнам в конеч-

дывалась. В некоторых симуляциях попутно наблю-

ном бассейне данное явление было впервые промо-

далось формирование блокированных волн на тех

делировано в работах [40, 41]. Хотя дно бассейна в

местах, где крупномасштабное течение переходило

данном случае предполагается горизонтальным, все

с меньшей глубины на большую.

остальные существенные ингредиенты метода рабо-

Более нетривиальное явление при взаимодейст-

тают в полной мере. Например, без учета времен-

вии волн с неоднородным дном — так называемые

ной зависимости конформного модуля высокоточ-

брэгговские (или щелевые) солитоны. Если в про-

ные вычисления были бы невозможны. В этих чис-

странственно-одномерной волновой системе (а та-

ленных экспериментах в качестве начального усло-

ковой является свободная граница над двумерным

вия бралась неподвижная жидкость с вертикаль-

течением) волны взаимодействуют с периодической

ным отклонением свободной границы в виде поло-

структурой, то частотный спектр линейных волн де-

вины косинусоидальной волны с максимумом на од-

лится на зоны, между которыми появляются щели.

ном конце бассейна (при x = 0) и с минимумом —

Если учесть нелинейность, то в некоторых случаях

на другом (при x = L/2). Далее волновая система

оказываются возможными долгоживущие локализо-

эволюционировала по типичному для ФПУ-явления

ванные структуры в виде солитонов огибающей для

сценарию — несколько медленных колебаний в ре-

стоячих волн. При этом частота солитона оказыва-

жиме стоячей волны, затем формирование более ко-

ется в пределах щели. Применительно к волнам над

ротких когерентных структур — солитонов, их дли-

периодическими профилями дна (дно имеет пери-

тельная нетривиальная динамика, в течение кото-

од Λ вдоль x, а промодулированная стоячая волна

рой на короткое время появлялись приблизительно

имеет длину 2Λ) такие когерентные структуры бы-

стоячие волны с кратными волновыми числами. За-

ли рассмотрены в работах автора [38, 39]. Теорети-

тем время как бы поворачивало свой ход вспять, и

ческие предсказания, основанные на известных ра-

из менее упорядоченного состояния формировалось

нее решениях приближенных модельных уравнений,

более упорядоченное в виде все той же начальной

были успешно подтверждены прямым численным

косинусоиды. После нескольких стоячих колебаний

моделированием точных уравнений в рамках обсуж-

процесс приблизительно повторялся. В зависимости

даемого здесь метода. Брэгговские солитоны, задан-

от относительной длины бассейна L/2h и началь-

ные в начальный момент аналитической формулой,

ной амплитуды косинуса A₀/h количество взаимо-

не разрушались иногда в течение сотен волновых

действующих солитонов и период возврата были су-

периодов, хотя, конечно же, неучтенные слабонели-

щественно различными. Приблизительно время воз-

нейной теорией эффекты давали о себе знать в слу-

врата подгоняется формулой

чае больших амплитуд. В целом это явление выгля-

T_FPU (g/h)1/2 ≈ C(L/h)²(h/A₀)1/2,

дит довольно любопытно и парадоксально: локали-

зованная стоячая волна на воде не разбивается с те-

где C = 0.148 + 0.096(A₀/h). При этом точность воз-

чением времени на две группы волн, бегущих впра-

врата к (почти) неподвижному состоянию была на

во и влево, хотя, казалось бы, поверхность свободна

удивление высокой для не слишком длинных бассей-

и этому ничто не мешает. Весьма существенно, что

нов L/h ∼ 60-100 и не слишком больших амплитуд

наиболее благоприятный режим для наблюдения та-

A₀/h ∼ 0.12-0.15.

ких щелевых солитонов получается при очень силь-

Кроме того, в работе [41] брались начальные со-

ных неоднородностях дна (по сути, это «гребенка»

стояния в виде нескольких солитонов и моделирова-

из узких и высоких барьеров), которые едва ли мо-

лись экстремальные волны над различными профи-

гут быть столь же адекватно и легко обработаны

лями дна, включая длинноволновые и коротковол-

другми методами.

новые неоднородности. В случае начального состо-

Другой пример парадоксальной нелинейной ди-

яния из нескольких разделенных солитонов над го-

намики — знаменитый возврат Ферми - Паста - Ула-

ризонтальным дном качество возврата было в целом

ма (ФПУ-возврат), когда начальное состояние вол-

значительно хуже, чем для косинуса. Что касается

новой системы в виде одной возбужденной про-

экстремальных волн, то плоское дно и соответству-

странственной гармоники сперва трансформирует-

ющий квазиинтегрируемый режим динамики не спо-

ся нелинейностью в совокупность солитонов, а за-

собствовали их появлению. Если же приблизитель-

тем после долгого периода нетривиальных взаимо-

ная интегрируемость была нарушена неоднородно-

действий между ними волновая энергия снова ока-

стями дна, то система переходила в состояние слу-

950

ЖЭТФ, том 157, вып. 5, 2020

Волны над искривленным дном. . .

, r

max

min

лем дна, тогда как при относительно коротко-кор-

0.80

δ = 0

релированных неоднородностях дна экстремальные

0.1

волны были менее высокими, но зато имели более

0.70

-0.1

резкий гребень. Похожие эффекты наблюдались как

0.60

для волн между двумя вертикальными стенками,

так и для волн с периодическими граничными усло-

0.50

виями без стенок.

m = 1

0.40

Тематика аномальных волн была затронута и в

работе [42], но в несколько иной постановке. Ана-

, r

Время

max

min

литические оценки для волн на глубокой воде в

0.80

δ = 0

присутствии крупномасштабного неоднородного те-

0.1

чения показали, что встречное усиливающееся те-

0.70

-0.1

чение делает волну более модуляционно неустойчи-

0.60

вой и тем самым способствует формированию ано-

мальных волн. Чтобы пронаблюдать этот эффект

0.50

в прямом численном эксперименте, автор включил

m = 2

потенциальное течение (аналог Γ) над неоднород-

0.40

ным профилем дна и промоделировал относитель-

Время

, r

но короткие бегущие волны. Прямое взаимодействие

max

min

таких волн с дном пренебрежимо мало, но суще-

0.80

δ = 0

0.1

ственным было присутствие течения — медленно-

0.70

-0.1

го над глубоким дном и более быстрого над мел-

ким дном. В начальный момент времени в области с

0.60

медленным течением запускался достаточно длин-

0.50

ный промодулированный волновой пакет, который

m = 3

затем распространялся против течения. При выхо-

0.40

де на сильное встречное течение развивалась мо-

Время

дуляционная неустойчивость и образовывались ано-

rmax^{, r}min

мальные волны, в соответствии с теорией. Было так-

0.80

δ = 0

же отмечено, что если на встречное течение выхо-

0.1

дит не слабопромодулированный длинный волновой

0.70

-0.1

пакет, а квазислучайная последовательность групп

0.60

волн, то могут возникать гораздо более высокие вол-

ны-убийцы, чем предсказывает формула, основан-

0.50

ная на решении нелинейного уравнения Шрединге-

m = 4

ра в виде так называемого бризера Ахмедиева. При-

0.40

чина появления более высоких аномальных волн за-

Время

ключается в притягивающем взаимодействии квази-

Рис. 1. Временные зависимости максимального и мини-

солитонных когерентных структур, в каковые пре-

мального значений радиальной координаты свободной по-

вращаются типичные волновые группы при выходе

верхности в слегка нециркулярных центрифугах для раз-

на быстрое встречное течение, тогда как бризер Ах-

личных азимутальных чисел неровности дна m и различ-

медиева и основанная на нем формула не учитыва-

ных отстроек δ угловой скорости вращения от соответству-

ют возможных процессов слияния квазисолитонов.

ющих резонансных значений

чайного волнового поля, где присутствовали квази-

4. НОВОЕ ПРИМЕНЕНИЕ: ВОЛНЫ В

солитонные когерентные структуры различных ам-

ЦЕНТРИФУГАХ

плитуд, причем некоторые из них были значительно

сильнее, чем начальные солитоны. Когда наиболее

Прежде чем перейти к новым численным ре-

сильные противоположно бегущие солитоны сталки-

зультатам, рассмотрим линеаризованные уравнения

вались, возникали экстремальные волны. Наиболее

движения для поверхностных волн в частично за-

высокие волны наблюдались над плавным профи-

полненной, идеально круглой центрифуге (единич-

951

В. П. Рубан

ЖЭТФ, том 157, вып. 5, 2020

y_lab

1.0

t = 42.79

t = 19.84

0.5

-0.5

m = 1

m = 2

-1.0

–1.0

-0.5

0.5

1.0

-1.0

-0.5

0.5

1.0

x_lab

y_lab

1.0

t = 15.77

t = 16.61

0.5

-0.5

m = 3

m = 4

-1.0

-0.5

0.5

1.0

-1.0

-0.5

0.5

1.0

x_lab

Рис. 2. Формирование особенностей на свободной поверхности при δ = 0 для различных m

ного радиуса и с единичной угловой скоростью)

ψ_t = cη - 2^Rψ - cg(t)sin(χ + t),

(35)

при условии относительной малости силы тяжести

η_t = c-1 Rψ_χ,

(36)

(g/Ω²R ≪ 1 в размерных переменных). Будем ис-

пользовать полярные координаты r и χ, так что x +

где оператор

R (аналог оператора

R_α) диагонален в

+iy = r exp(iχ). Неизвестными функциями при этом

дискретном фурье-представлении:

будут малое отклонение свободной границы η(χ, t) =

= [r(χ, t) - c] от равновесного радиуса 0 < c < 1 и

(

)

1-c2m

граничное значение ψ(χ, t) потенциала поля скорос-

R_m = i

= ith mln

≡ iσ_m.

(37)

1+c2m

ти. Обобщенное уравнение Бернулли и кинемати-

ческое граничное условие на свободной границе в

Соотношение θ(χ, t) =

^Rϕ(χ, t) при r = c возник-

главном приближении выглядят следующим обра-

ло из того факта, что выполнение кинематического

зом (при Γ = 0):

условия на стенке контейнера обеспечивается равен-

952

ЖЭТФ, том 157, вып. 5, 2020

Волны над искривленным дном. . .

rmax^{, r}min

ством

0.80

δ = 0

φ_m = A_m(t)r^me^imχ +

A_m(t)r^-me^-imχ.

(38)

0.75

ωvibr=δ+ 1-ω -1

0.1

0.70

-0.1

В случае g = const частное решение неоднородной

0.65

линейной системы (35), (36) имеет вид

0.60

0.55

ψ_st =

cos(χ + t),

(39)

0.50

1+σ₁

0.45

-σ₁g

0.40

η_st =

sin(χ + t).

(40)

100

1+σ₁

Время

Оно соответствует стационарному состоянию в ла-

Рис. 3. Максимальное и минимальное значения ради-

бораторной системе координат (см. подробности в

альной координаты свободной поверхности в зависимос-

работе [43]). Из (35), (36) также следует выражение

ти от времени при вертикальной вибрации оси круг-

для собственных частот («закон дисперсии»; попут-

лой центрифуги с эффективной силой тяжести g(t) =

но заметим, что в работе [43] дисперсионное уравне-

= 0.02 - 0.01 cos(ω_vibr t)

ние выводилось в невращающейся системе, и, веро-

ятно, поэтому ответ не был доведен до столь прос-

y_lab

той формулы)

√

ω_m = σ_m +

σ2m + mσ_m ,

(41)

1.0

t = 53.31

причем положительным m соответствуют опережа-

ющие волны, а отрицательным m — отстающие (в

лабораторных координатах). Ясно, что если в систе-

0.5

му внести возмущения с азимутальным числом m и

относительной скоростью вращения Δ ≈ ω_m/m, то

должен будет наблюдаться резонансный рост соот-

ветствующей гармоники. Такие возмущения проще

всего осуществить, слегка исказив форму центрифу-

ги и заставив ее вращаться не с частотой Ω = 1, а с

-0.5

частотой вблизи 1 + Δ. В численных примерах, при-

веденных на рис. 1 и 2, такой режим обеспечивался

функцией

-1.0

Z(ζ, t) = exp[i(-ω_-m/m + δ)t]e^iζ [1 + ϵe^imζ ]

(42)

-1.0

-0.5

0.5

1.0

x_lab

при ϵ = 0.02 для m = 1, 2, 3, 4. Отстройка от резо-

нанаса фиксировалась параметром δ. В качестве на-

Рис. 4. Образование особенности на свободной поверх-

чальных условий брались функции ρ = 0 и ψ = 0, а

ности в круглой центрифуге с вертикально вибрирующей

начальное значение α₀ = - ln(0.6) давало приблизи-

осью при резонансном значении частоты вибрации

тельно круглую форму свободной границы со сред-

ним радиусом c ≈ 0.6. В этой серии численных экс-

периментов использовалось значение безразмерного

Наиболее эффективно такие возмущения воздейст-

поля тяжести g = 0.02.

вуют на гармоники с m = ±1, причем даже в слу-

На рис. 1 и 2 видно, что действительно имел ме-

чае идеально круглого контейнера, как это следует

сто резонансный рост волн, сопровождавшийся уси-

из приближенных уравнений (35), (36). Резонансные

лением их нелинейности и закончившийся образова-

частоты при этом даются формулой ω_vibr = 1 ± ω±1.

нием особенности на свободной поверхности. В ре-

Соответствующие примеры представлены на рис. 3

альных условиях это означало бы начало перехода

и 4.

течения в трехмерный турбулентный режим.

Чтобы продемонстрировать возможности мето-

Другая возможность привнесения резонансных

да в полной мере, на рис 5 и 6 показаны примеры

возмущений — вертикальные вибрации оси враще-

эволюции свободной поверхности в деформируемой

ния, приводящие к зависимости

центрифуге с изменяющейся (уменьшающейся) со

g(t) = g₀ + A cos(ω_vibr t).

временем площадью сечения. В этих экспериментах

953

В. П. Рубан

ЖЭТФ, том 157, вып. 5, 2020

t = 2.00

t = 14.00

t = 2.00

t = 12.00

c = 0.6

c = 0.7

t = 24.00

t = 24.735

t = 20.00

t = 21.294

t = 2.00

t = 10.00

c = 0.8

Рис. 5. Образование особенностей на свободной поверхнос-

t = 17.00

t = 17.422

ти в центрифуге, деформируемой согласно формуле (43) с

m = 1, C1 = 0.9, C2 = 0.7, для трех значений начального

радиуса

g = 0.1, а форма контейнера соответствовала анали-

поверхности c. Затем часть дна контейнера как бы

тической функции

«поднималась» вдоль радиуса, образуя при этом m

участков с «меньшей глубиной». Уменьшение пло-

Z(ζ, t) =

щади контейнера здесь сопровождалось появлени-

[

]

ем дополнительного кругового потенциального те-

⁽1+C₁e^imζ )

= exp iζ+id(t)ln

-0.5πd(t) ,

(43)

чения с параметром Γ(t), направленного против ча-

1-C₂eimζ

совой стрелки. Особенности (опрокидывающиеся уг-

причем мера отклонения от окружности определя-

лы) формировались на гребне растущей волны вбли-

лась выражением

зи места выхода течения с меньшей глубины на

большую.

d(t) = 0.1[1 - exp(-0.04t)].

Наконец, на рис 7 показана деформация изна-

На рис. 5 параметры m = 1, C₁ = 0.9, C₂ = 0.7. На

чально круглого контейнера к форме, промежу-

рис. 6 было взято m = 2, C₁ = 0.9, C₂ = 0.9. В на-

точной между кругом и квадратом (при g = 0.1,

чальный момент времени это была циркулярно сим-

c = 0.6). В этом численном эксперименте исполь-

метрическая конфигурация с радиусом свободной

зовалось разложение соответствующего эллиптиче-

954

ЖЭТФ, том 157, вып. 5, 2020

Волны над искривленным дном. . .

(

)

μ²

5μ³

35μ⁴

E₄(μ) =

(45)

16 · 13

128 · 17

чтобы построить функцию

)

t = 2.00

t = 14.00

Z(ζ, t) = [1 - 0.

d(t)]e^iζ E₄

d(t) exp(4iζ)

(46)

где

d(t) = 0.8[1-exp(-0.04t)]. И в этом случае имело

место образование особенности. Но если вместо мно-

жителя [1 - 0.

d(t)] брался множитель [1 - 0.

d(t)],

c = 0.6

который при деформации давал более существенное

уменьшение площади контейнера и полости в нем,

то (при c = 0.6) волны на свободной поверхности в

течение длительного времени оставались гладкими

t = 20.00

t = 21.858

и не проявляли тенденции к формированию резких

опрокидывающихся гребней (не показано). Полного

понимания причин такого поведения пока нет.

5. ЗАКЛЮЧЕНИЕ

Рис. 6. Образование особенности на свободной поверхно-

сти в центрифуге, деформируемой согласно формуле (43)

Таким образом, основной результат данной рабо-

с m = 2, C1 = 0.9, C2 = 0.9

ты заключается в разработке и применении метода

составного конформного отображения для описания

динамики волн на свободной поверхности идеаль-

ной несжимаемой жидкости в частично заполнен-

ных плоских центрифугах сложной формы. Как и в

предыдущих применениях, метод продемонстриро-

вал свою высокую точность и эффективность. Прав-

t = 5.00

t = 30.00

да, пока нельзя сказать, что все вопросы исчерпаны.

Например, пока не ясно, как следует модифициро-

вать метод в том случае, когда относительно малая

полость (полый вихрь при Γ = 0) уходит в процессе

c = 0.6

своей динамики далеко от начала координат и она

оказывается занята жидкостью. Очевидно, что кон-

формное отображение типа Z(ζ) = exp(iζ), у кото-

рого при Z = 0 (или в любой другой заранее фикси-

рованной точке) имеется особенность, при этом пе-

t = 45.00

t = 50.68

рестает работать.

Но даже в уже имеющемся варианте метод спо-

собен решать множество интересных задач. Совер-

шенно ясно, что приведенные здесь примеры дале-

ко не исчерпывают всего разнообразия возможных

Рис. 7. Образование особенности на свободной поверхно-

волновых структур и их динамики во вращающихся

сти в центрифуге, деформируемой от круга к квадрату по

системах с нетривиальной геометрией.

формуле (46)

ЛИТЕРАТУРА

ского интеграла до четвертого порядка,

1. Л. В. Овсянников, Динамика сплошной среды, Но-

восибирск, Наука (1973), вып. 15, с. 104.

∫

√

≈ zE₄(z⁴),

(44)

2. L. V. Ovsjannikov, Arch. Mech. 26, 407 (1974).

1-z⁴

3. A. I. Dyachenko, E. A. Kuznetsov, M. D. Spector,

and V. E. Zakharov, Phys. Lett. A 221, 73 (1996).

955

В. П. Рубан

ЖЭТФ, том 157, вып. 5, 2020

A. I. Dyachenko, Y. V. Lvov, and V. E. Zakharov,

25.

M. R. Turner and T. J. Bridges, Adv. Comput. Math.

Physica D 87, 233 (1995).

43, 947 (2016).

А. И. Дьяченко, В. Е. Захаров, Е. А. Кузнецов,

26.

M. R. Turner, J. Fluids Struct. 64, 1 (2016).

Физика плазмы 22, 916 (1996).

27.

M. G. Blyth and E. I. Parau, J. Fluid Mech. 806, 5

А. И. Дьяченко, ДАН 376, 27 (2001).

(2016).

V. E. Zakharov, A. I. Dyachenko, and O. A. Vasilyev,

28.

P. M. Lushnikov, S. A. Dyachenko, and D. A. Silan-

Eur. J. Mech. B/Fluids 21, 283 (2002).

tyev, Proc. R. Soc. A 473 (2202), 2017019 (2017).

A. I. Dyachenko and V. E. Zakharov, Письма в

29.

A. I. Dyachenko, P. M. Lushnikov, and V. E. Zakha-

ЖЭТФ 81, 318 (2005).

rov, J. Fluid Mech. 869, 526 (2019).

V. E. Zakharov, A. I. Dyachenko, and A. O. Proko-

30.

A. I. Dyachenko, S. A. Dyachenko, P. M. Lushnikov,

fiev, Eur. J. Mech. B/Fluids 25, 677 (2006).

and V. E. Zakharov, J. Fluid Mech. 874, 891 (2019).

10.

A. I. Dyachenko and V. E. Zakharov, Письма в

31.

S. A. Dyachenko, J. Fluid Mech. 860, 408 (2019).

ЖЭТФ 88, 356 (2008).

32.

T. Gao, A. Doak, J.-M. Vanden-Broeck, and Zh.

11.

W. Choi and R. Camassa, J. Engng Mech. 125, 756

Wang, Eur. J. Mech./B Fluids 77, 98 (2019).

(1999).

33.

M. V. Flamarion, P. A. Milewski, and A. Nachbin,

12.

Y. A. Li, J. M. Hyman, and W. Choi, Stud. Appl.

Stud. Appl. Math. 142, 433 (2019).

Math. 113, 303 (2004).

34.

D. Kachulin, A. Dyachenko, and A. Gelash, Fluids 4,

13.

Р. В. Шамин, ДАН 406, 112 (2006).

83 (2019).

14.

Р. В. Шамин, ДАН 418, 603 (2008).

35.

V. P. Ruban, Phys. Rev. E 70, 066302 (2004).

15.

Р. В. Шамин, ДАН 432, 458 (2010).

36.

V. P. Ruban, Phys. Lett. A 340, 194 (2005).

16.

В. Е. Захаров, Р. В. Шамин, Письма в ЖЭТФ 91,

37.

V. P. Ruban, Phys. Rev. E 77, 037302 (2008).

68 (2010).

38.

V. P. Ruban, Phys. Rev. E 77, 055307(R) (2008).

17.

В. Е. Захаров, Р. В. Шамин, Письма в ЖЭТФ 96,

39.

V. P. Ruban, Phys. Rev. E 78, 066308 (2008).

68 (2012).

40.

V. P. Ruban, Письма в ЖЭТФ 93, 213 (2011).

18.

D. Chalikov and D. Sheinin, J. Comput. Phys. 210,

247 (2005).

41.

V. P. Ruban, ЖЭТФ 141, 387 (2012).

19.

D. Chalikov, Phys. Fluids 21, 076602 (2009).

42.

В. П. Рубан, Письма в ЖЭТФ 95, 550 (2012).

20.

А. В. Слюняев, ЖЭТФ 136, 785 (2009).

43.

O. M. Phillips, J. Fluid Mech. 7, 340 (1960).

21.

W. Choi, Math. Comp. Simulat. 80, 29 (2009).

44.

A. A. Ivanova, V. G. Kozlov, and A. V. Chigrakov,

22.

S. A. Dyachenko, P. M. Lushnikov, and A. O. Korot-

Fluid Dynamics 39, 594 (2004).

kevich, Письма в ЖЭТФ 98, 767 (2013).

45.

A. A. Ivanova, V. G. Kozlov, and D. A. Polezhaev,

23.

Н. М. Зубарев, Е. А. Кочурин, Письма в ЖЭТФ

Fluid Dynamics 40, 297 (2005).

99, 729 (2014).

24.

P. M. Lushnikov, J. Fluid Mech. 800, 557 (2016).

46.

http://www.fftw.org.

956