Известия РАН. Механика жидкости и газа, 2019, № 3, стр. 83-97

Процесс формирования внутренних волн, инициированный началом движения тела в стратифицированной вязкой жидкости

П. В. Матюшин *

Институт автоматизации проектирования РАН
Москва, Россия

* E-mail: pmatyushin@mail.ru

Поступила в редакцию 11.04.2018
После доработки 16.10.2018
Принята к публикации 18.10.2018

Полный текст (PDF)

Аннотация

Впервые детально рассмотрен механизм формирования пространственных гравитационных внутренних волн, инициированный началом движения диска заданного диаметра и конечной толщины в горизонтальном направлении вдоль оси симметрии диска справа налево с заданной постоянной скоростью в линейно стратифицированной по плотности несжимаемой вязкой жидкости благодаря численному решению системы уравнений Навье–Стокса в приближении Буссинеска и визуализации пространственной вихревой структуры рассчитанного течения. Полученные поля векторов скоростей и возмущения давления обладают горизонтальной и вертикальной плоскостями симметрии, проходящими через ось симметрии диска. Описан процесс формирования течения в верхнем полупространстве, обусловленный сдвиговой и гравитационной неустойчивостями, – сначала между тыльным торцом диска и местом импульсного старта формируются две горизонтальные вихревые нити, трансформирующиеся в ножки шпилькообразной вихревой петли, головная часть которой расположена правее места старта. Далее над местом старта периодически в течение половины периода плавучести жидкости формируются вихревые кольца, левые половины которых преобразуются в полуволны, заполняющие пространство между диском и местом старта.

Ключевые слова: стратифицированная вязкая жидкость, диск, гравитационные внутренние волны, математическое моделирование

Понимание физики сложных трехмерных процессов генерации пространственных гравитационных внутренних волн движением тел в горизонтальном направлении в устойчиво стратифицированной вязкой сплошной среде океана и атмосферы важно как с практической, так и с теоретической точек зрения. Известно, что начало движения тела сопровождается излучением пучка нестационарных гравитационных внутренних волн, которые распространяются вдоль радиус-векторов от места старта Q тела [1, 2].

В настоящей работе при помощи математического моделирования течений линейно стратифицированной по плотности несжимаемой вязкой жидкости и визуализации пространственной вихревой структуры рассчитанных течений впервые детально рассмотрен сложный трехмерный механизм формирования волн, инициируемый началом движения диска диаметром d и толщиной h = 0.76d в горизонтальном направлении вдоль своей оси симметрии Z справа налево со скоростью U при 0.5 < Fr < 4 и Re = 50, где Fr = UTb/(2πd) – внутреннее число Фруда, Re = Ud/ν – число Рейнольдса, Tb и ν – период плавучести и коэффициент кинематической вязкости жидкости соответственно. В [3] описана начальная стадия этого механизма при T < 0.5 в вертикальной плоскости симметрии течения для Fr = 1, Re = 500; T – обезразмеренное на Tb время, прошедшее с начала старта диска.

Традиционно экспериментаторы, исследуя течения стратифицированной вязкой жидкости, рассматривают картину гравитационных внутренних волн и след около двигающихся тел (сфер и цилиндров) [2, 47] или обтекаемых препятствий на поверхности нашей планеты [8] только в вертикальной плоскости, тогда как математическое моделирование дает трехмерную вихревую структуру внутренних волн [3, 911]. Механизм формирования этой структуры экспериментально не исследован. В [10] приведена классификация течений стратифицированной вязкой жидкости около сферы при 0.005 < Fr < 100 и 1 < Re < 500, которая хорошо согласуются с экспериментом [6]. Результаты [12] для сферы при Re = 200, 0.125 < Fr < 100 хорошо согласуются с экспериментом [6] при 0.25 < Fr < 100. При Fr < 0.25 в [6] наблюдается нестационарное периодическое течение в следе за сферой, а в [12] – стационарное течение. Классификация режимов течений около диска толщиной h = 0.76d при 0.05 < Fr < 100 и 50 < Re < 500 приведена в [11].

Для получения начального представления о трехмерном механизме формирования пространственных гравитационных внутренних волн может служить задача о течении, индуцированном диффузией на сфере, помещенной в покоящуюся непрерывно стратифицированную вязкую жидкость [13]. Рассматривалось течение, осесимметричное относительно вертикальной прямой q, проходящей через центр Q сферы, при симметрии поля векторов скорости еще и относительно горизонтальной плоскости, проходящей через точку Q. При T ≤ 0.5 в верхнем полупространстве генерируется осесимметричное вихревое кольцо, заполняющее все это полупространство. При T > 0.5 в течение каждого ∆T = 0.5 в окрестности прямой q над сферой генерируется новое вихревое кольцо, которое уменьшает вертикальные размеры ранее сгенерированных колец. Каждая пара колец формирует одну гравитационную внутреннюю волну. Групповая скорость этих волн перпендикулярна их фазовой скорости и направлена по радиус-вектору от каждого из двух эффективных центров волнообразования – полюсов сферы [13]. Волновая энергия распространяется в радиальных направлениях от полюсов сферы с групповой скоростью параллельно гребням волн. Со временем на больших удалениях от тела фазовые поверхности с постоянной угловой скоростью стремятся к горизонтальным плоскостям, проходящим через центры волнообразования. При T > 500 около этих горизонтальных плоскостей останутся только два вихревых кольца, сильно сплуснутых в вертикальном направлении и наблюдаемых в эксперименте [13].

Для понимания пространственного механизма формирования гравитационных внутренних волн служит и двухмерная задача о равномерном движении бесконечно длинного горизонтально ориентированного цилиндра в горизонтальной направлении справа налево, перпендикулярно своей оси симметрии Z (при Re < 200 и Fr < 1) [2, 3, 9, 10]. Здесь около места Q старта цилиндра появляются два центра волнообразования, как и в эксперименте [13]. Поле векторов скорости здесь будет симметричным относительно горизонтальной прямой, проходящей через точку Q. В течение каждого ∆T = 1, справа и слева от вертикальной прямой q, проходящей через Q, в верхней и нижней полуплоскостях генерируется по одной новой внутренней волне [3, 9, 10]. Волны левее прямой q двигаются влево вместе с цилиндром, а волны правее прямой q остаются на месте и сжимаются в вертикальном направлении под действием вновь нарождающихся волн. Картины изолиний горизонтального градиента плотности около квадратного цилиндра, полученные в результате математического моделирования в [3, 9, 10] при Fr = 0.1, очень хорошо согласуются с теневыми картинами “вертикальная щель – нож Фуко” течений около кругового цилиндра, полученными в эксперименте [2] при Fr = 0.094.

Настоящая статья является логическим продолжением работ автора [3, 9, 10, 13], посвященным анализу механизма формирования гравитационных внутренних волн.

1. ПОСТАНОВКА ЗАДАЧИ И МЕТОД ЕЕ РЕШЕНИЯ

Рассмотрим равномерное течение линейно стратифицированной по плотности несжимаемой вязкой жидкости в горизонтальном направлении Z слева направо со скоростью U. В некоторый момент времени в это течение мгновенно вносится диск с диаметром d, толщиной h = 0.76d и с горизонтальной осью симметрии Z (рис. 1а), и ставится задача об изменении картины течения. Для ее решения в геометрическом центре диска помещается начало неподвижной декартовой системы координат (СК1) (X, Y, Z), где ось X – вертикальна. Плотность жидкости ρ(X, Y, Z) = = $1 - 0.5 \cdot {X \mathord{\left/ {\vphantom {X {\text{A}}}} \right. \kern-0em} {\text{A}}} + S(X,Y,Z)$ обезразмерена на плотность на уровне центра диска ρ0, а координаты X, Y, Z – на d/2; N = 2π/Tb и $\Lambda = g{\text{/}}{{N}^{2}}$ – частота и масштаб плавучести жидкости; ${\text{A}} = \Lambda {\text{/}}d$ – отношение масштабов; g – ускорение свободного падения; S – обезразмеренное на ρ0 возмущение солености, которое в начале расчета равно нулю. Значения параметра A > 100 выбрались для выполнения реализуемых в эксперименте условий N ≈ 1 с–1, 0.1 см < d < 10 см.

Рис. 1.

(а) Постановка задачи об исследовании течения жидкости около диска. (б-н) Визуализация установившегося со временем течения при Fr = 0.5, Re = 50, А = 981.6, T = 11.46 (Сd = 2.649) в пространстве (б–в), в вертикальной плоскости XZ (г–к) и на диске (л–н): (б–в) – изоповерхности β = 0.02; г–ж – изолинии β с шагом 0.01 (г), Sz с шагами 4 × 10–5 и 10–5 (д и е) и S с шагом 10–5 (ж); з–к – линии тока в СК1 (з-и) и СК2 (к); л–н – линии трения передней (л), боковой (м) и тыльной (н) поверхностях диска в СК1

Для математического моделирования поставленной задачи решается система уравнений Навье–Стокса в приближении Буссинеска, записанная в цилиндрической СК (Z, R, φ): Z = Z, $X = R\cos \varphi $, $Y = R\sin \varphi $. В декартовой СК1 эта система уравнений выглядит следующим образом

(1.1)
$\frac{{\partial S}}{{\partial t}} + ({\mathbf{v}} \cdot \nabla {\text{)}}S = \frac{2}{{{\text{ScRe}}}}\Delta S + \frac{{{{\text{v}}_{X}}}}{{2{\text{A}}}}$
(1.2)
$\frac{{\partial {\mathbf{v}}}}{{\partial t}} + ({\mathbf{v}} \cdot \nabla {\text{)}}{\mathbf{v}} = - \nabla p + \frac{2}{{{\text{Re}}}}\Delta {\mathbf{v}} + \frac{{\text{A}}}{{2{\text{F}}{{{\text{r}}}^{2}}}}S\frac{{\mathbf{g}}}{g}$
(1.3)
$\nabla \cdot {\mathbf{v}} = 0$
где ${\mathbf{v}} = \left( {{{v}_{X}},{{v}_{Y}},{{v}_{Z}}} \right)$ – вектор скорости; p – возмущение давления, обезразмеренное на ${{\rho }_{0}}U_{{}}^{2}$; t – время (обезразмеренное на $f = d{\text{/}}\left( {2U} \right) = 1{\text{/}}\left( {2{\text{Fr}}N} \right)$); Sc = ν/κ = 709.22 – число Шмидта; κ – коэффициент диффузии соли; ∇ и Δ – операторы Гамильтона и Лапласа.

Для решения этой задачи использовался численный метод расщепления по физическим факторам МЕРАНЖ [14], который успешно применялся для моделирования течений несжимаемой вязкой жидкости около сфер, цилиндров и дисков [3, 911, 13, 15, 16], а также для течений со свободной поверхностью [14].

Использовалась цилиндрическая расчетная сетка [Z, R, φ] = $\left[ {I\, \times \,J \times K} \right]$ = $\left[ {240\, \times \,110 \times 40} \right]$, занимающая область [–13 < Z < 50, 0 < R < 30, 0 < φ < 2π]. Сетка сгущается ко всем поверхностям диска и к оси Z так, чтобы на скоростной пограничный слой приходилось 5 ячеек сетки. При удалении от поверхности обтекаемого тела вдоль направлений Z или R длины сторон ячеек сетки вдоль Z и R монотонно увеличиваются по полиномиальному закону до некоторого заданного максимального значения и далее не меняются. Минимальные (максимальные) размеры ячеек сетки для Re = 50 вдоль направлений Z и R равны 0.036 и 0.035 (0.41 и 0.67) соответственно.

Применялся разнесенный сеточный шаблон, когда переменные S и p определяются в центрах расчетных ячеек, а компоненты скорости – в центрах их граней; переменные S и p на поверхности тела не определялись.

Приведем алгоритм расчета для декартовой СК1. В начальный момент времени ${{t}_{0}} = {\text{0}}$ задавалось некоторое начальное течение около диска. Пусть при ${{t}_{n}} = n{\tau }$, где τ – величина шага по времени, n – номера шагов (n = 0, 1, 2, 3, …), известны значения S, ${\mathbf{v}} = \left( {{{v}_{X}},{{v}_{Y}},{{v}_{Z}}} \right)$ и p. Тогда схему нахождения неизвестных функций S, v и p в момент времени ${{t}_{{n + 1}}} = \left( {n + 1} \right) \cdot {\tau }$ можно представить в виде

(1.4)
$\frac{{{{S}^{{n + 1}}} - {{S}^{n}}}}{{\tau }} = - ({{{\mathbf{v}}}^{n}} \cdot \nabla {\text{)}}{{S}^{n}} + \frac{2}{{{\text{ScRe}}}}\Delta {{S}^{n}} + \frac{{v_{X}^{n}}}{{2{\text{A}}}}$
(1.5)
$\frac{{{\mathbf{w}} - {{{\mathbf{v}}}^{n}}}}{{\tau }} = - ({{{\mathbf{v}}}^{n}} \cdot \nabla {\text{)}}{{{\mathbf{v}}}^{n}} + \frac{2}{{{\text{Re}}}}\Delta {{{\mathbf{v}}}^{n}} + \frac{{\text{A}}}{{2{\text{F}}{{{\text{r}}}^{2}}}}{{S}^{{n + 1}}}\frac{{\mathbf{g}}}{g}$
(1.6)
$\tau \,\Delta p = \nabla \cdot {\mathbf{w}}$
(1.7)
$\frac{{{{{\mathbf{v}}}^{{n + 1}}} - {\mathbf{w}}}}{\tau } = - \nabla p$
где w – вспомогательная промежуточная скорость.

Сложение уравнений (1.5) и (1.7) дает уравнение (1.2). Уравнение (1.6) получается при скалярном умножении оператора Гамильтона на уравнение (1.7), учитывая уравнения неразрывности (1.3).

Для аппроксимации конвективных членов уравнений (1.4–1.5) использовалась гибридная конечно-разностная схема, для которой характерны второй порядок аппроксимации по пространственным переменным, минимальная схемная вязкость и дисперсия, работоспособность в широком диапазоне Re и Fr и монотонность [14]. Для аппроксимации других пространственных производных уравнений (1.4–1.7) – центральные разности.

Рассмотрим аппроксимацию уравнения Пуассона (1.6) на входной части внешней границы цилиндрической расчетной области (Z = –13)

(1.8)
$\tau \frac{{\left( {{{p}_{{{\text{2,}}j,k}}} - {{p}_{{{\text{1,}}j,k}}}} \right) - \left( {{{p}_{{{\text{1,}}j,k}}} - {{p}_{{{\text{0,}}j,k}}}} \right)}}{{{{h}^{2}}}} + \Lambda = \frac{{{{w}_{{z,{\text{2,}}j,k}}} - {{w}_{{z{\text{,1,}}j,k}}}}}{h} + \Lambda ,$
где h – шаг сетки, ${{p}_{{0{\text{,}}j,k}}}$ и ${{p}_{{1{\text{,}}j,k}}}$ – значения давления в центрах ячеек “0, j, k” и “1, j, k”, расположенных левее и правее по отношению к границе соответственно (j = 1, 2,…, J; k = 1, 2, …, K). В (1.8) граничные значения ${{p}_{{0{\text{,}}j,k}}}$ и ${{w}_{{z,1{\text{,}}j,k}}}$ неизвестны.

Выпишем выражение для ${{w}_{{z,1{\text{,}}j,k}}}$, полученное из уравнения (1.7)

${{w}_{{z,{\text{1,}}j,k}}} = v_{{z{\text{,1,}}j,k}}^{{n + 1}} + \tau \frac{{\left( {{{p}_{{{\text{1,}}j,k}}} - {{p}_{{{\text{0,}}j,k}}}} \right)}}{h}$

Подставим это выражение для ${{w}_{{z,1{\text{,}}j,k}}}$ в уравнение (1.8)

(1.9)
$\tau \frac{{\left( {{{p}_{{{\text{2,}}j,k}}} - {{p}_{{{\text{1,}}j,k}}}} \right)}}{{{{h}^{2}}}} + \Lambda = \frac{{{{w}_{{z,{\text{2,}}j,k}}} - v_{{z,1,j,k}}^{{n + {\text{1}}}}}}{h} + \Lambda .$

Таким образом, на входной части внешней границы цилиндрической расчетной области (Z = –13) вместо уравнения (1.8) можно использовать (1.9), где уже нет неизвестных граничных значений ${{p}_{{0{\text{,}}j,k}}}$ и ${{w}_{{z,1{\text{,}}j,k}}}$. Аналогично можно показать, что для всех оставшихся границ граничные условия для давления не нужны.

На твердой поверхности диска ставятся граничные условия прилипания и непротекания: v = $({{v}_{z}},{{v}_{R}},{{v}_{\varphi }}) = (0,0,0)$, $\partial {\rho /}\partial n = 0$, где n – вектор нормали к поверхности диска (на передней и тыльной сторонах диска: $\partial S{\text{/}}\partial z = 0$; на боковой стороне диска: $\partial S{\text{/}}\partial R = 0.5\cos \varphi {\text{/A}}$). На входной части внешней границы цилиндрической расчетной области (Z = –13): ${\mathbf{v}} = \left( {1,\,0,\,0} \right)$, $\partial S{\text{/}}\partial z = 0$. На боковой части внешней границы (R = 30): ${\mathbf{v}} = \left( {1,\,0,\,0} \right)$, $\partial S{\text{/}}\partial R = 0$. На выходной части внешней границы (Z = 50): ${{v}_{z}} = 1$, ${{v}_{R}} = 0$, $\partial {{v}_{\varphi }}{\text{/}}\partial z = 0$, $\partial S{\text{/}}\partial z = 0$. На бесконечном удалении от тела S = 0, но в силу ограниченных размеров расчетной сетки на ее внешней границе ставятся “свободные” граничные условия: $\partial S{\text{/}}\partial n = 0$. На выходной части внешней границы естественно поставить “свободные” граничные условия $\partial {{v}_{\varphi }}{\text{/}}\partial z = 0$ для третьей компоненты вектора скорости. К сожалению, постановка аналогичных “свободных” граничных условий для всех компонент скорости при Z = 50 приводит к аварийной остановке процесса вычислений.

Созданный программный комплекс математического моделирования и визуализации пространственных течений стратифицированной вязкой жидкости около диска был всесторонне протестирован [3, 11]. Для Re = Fr = 50, A = 9816 коэффициент сопротивления диска равен Сd = 1.923, длина рециркуляционной (застойной) области D1 за диском, отсчитываемая от тыльной критической точки диска, равна L/d = 0.741, что хорошо согласуется с экспериментом [17]. Расчеты проводились на вычислительных ресурсах Межведомственного суперкомпьютерного центра Российской академии наук (МСЦ РАН).

2. ВИЗУАЛИЗАЦИЯ РЕЗУЛЬТАТОВ РАСЧЕТА

При Fr > 10 течение около диска будет эквивалентно течению однородной вязкой жидкости и представляет в СК1 след за телом. При Fr < 5 наблюдается генерация гравитационных внутренних волн, характеризующихся горизонтальной и вертикальной плоскостями симметрии, проходящими через ось Z [3, 911].

Рассмотрим рассчитанное установившееся трехмерное поле течения при Fr = 0.5, Re = 50, А = 981.6 (рис. 1). В силу горизонтальной плоскости Y–Z симметрии поля векторов скоростей этого течения анализируем только верхнее полупространство течения.

Линии трения (предельные линии тока) на поверхности диска [15, 16] при X > 0 показаны на рис. 1л–1н. Набегающий слева поток (рис. 1а) фокусируется в некоторой точке на передней стороне диска (ближе к ее верхней границе), а потом растекается по ней в разных направлениях (рис. 1л). На боковых сторонах диска поток стремится вниз к плоскости Y–Z (рис. 1м). В силу существования вертикальной плоскости симметрии X–Z поля векторов скоростей этого течения, в плоскости X–Z можно наблюдать линии тока (рис. 1з–1к). Более наглядна картина “синусоидальных” линий тока в СК1 (рис. 1з–1и). Здесь наблюдаются гребни и впадины волн. На рис. 1и расстояние между первым и вторым гребнем волны вдоль оси Z приблизительно равно λ, где ${\lambda } = U{{T}_{b}} = 2{\pi }U{\text{/}}N = 2{\pi Fr}d$ – длина внутренних волн в вертикальной плоскости X–Z.

В экспериментах, как правило, наблюдаются картины изолиний различных производных плотности. Например, популярная и информативная теневая картина “вертикальная щель – нож Фуко” [2] дает поле изолиний горизонтального градиента плотности Sz. Ниже для сравнения с экспериментом строятся изолинии Sz в вертикальной плоскости X–Z (рис. 1д–1е). Темные (Sz < 0) и светлые (Sz ≥ 0) полосы на рис. 1е визуализируют фазовые поверхности внутренних волн, а на их границах располагаются линии гребней и впадин, ясно видимые и на рис. 1и. В этом смысле рис. е и и подобны. В то же время поле изолиний Sz около тыльной стороны диска (рис. 1д) дает больше структурных элементов течения, чем картина линий тока на рис. 1з. Темные (S < 0) и светлые (S ≥ 0) полосы на картине изолиний S на рис. 1ж визуализируют полуволны впадин и гребней соответственно.

Поскольку в экспериментах тело обычно движется относительно покоящейся жидкости, то введем систему координат СК2 (x, y, z) (ось x – вертикальна, z – параллельна Z), которая равномерно двигается слева направо со скоростью U относительно СК1 (рис. 1а). В СК2 наблюдаются линии тока, для которых значения горизонтальных компонент векторов скорости, рассчитанные в СК1, уменьшаются на единицу – ср. з–и и к, а значения переменных S и p в СК2 те же, что и в СК1.

Картины мгновенных линий тока в СК2 визуализируют циркуляционные (вихревые) ячейки течения (рис. 1к). При t = 0 начало СК2 совпадает центром Q тыльной стороны диска. За время T = $[tf]{\text{/}}{{T}_{b}} = [t{\text{/}}(2{\text{Fr}}N)][N{\text{/}}2\pi ] = t{\text{/}}(4{\pi Fr})$ начало СК2 смещается в СК1 на расстояние s = $U[tf]{\text{/}}(0.5d) = t$ = 4πFrT.

Можно показать [9], что для установившегося течения в плоскости симметрии X–Z картина изолиний S (рис. 1, ж) подобна картине линий тока в СК2 (рис. 1к), где четко выделяются волнообразное движение жидкости около оси Z за равномерно движущимся справа налево тыловым торцом диска и ряд вытянутых циркуляционных ячеек над ним. Назовем часть жидкости между циркуляционной ячейкой номера M и осью Z на рис. 1к “базой” M этой ячейки. В циркуляционной ячейке 0, прилегающей к диску на рис. 1к, циркуляция жидкости идет по часовой стрелке, в следующей вытянутой “базе” 1, прилегающем справа к ячейке 0, жидкость описывает петлю, двигаясь против часовой стрелки; в ячейке 2, прилегающей справа к базе 1, циркуляция жидкости идет по часовой стрелке и т.д. Номера ячеек 0, 2–5 и базы 1 указаны на рис. 1ж. Таким образом, в ячейках 0, 2 и 4 на рис. 1к жидкость циркулирует по часовой стрелке, а в ячейках 3 и 5 (полуволнах впадин) – против часовой стрелки. На рис. 1ж и 1к четко прослеживаются базы 1, 3 и 5.

Для отображения пространственной вихревой структуры течения в каждом центре ячейки расчетной сетки определялась функция β. Если в центре ячейки существуют комплексно-сопряженные собственные значения ${{{\sigma }}_{{1{\text{,}}2}}} = {\alpha }\, \pm \,i{\gamma }$ тензора градиента скорости G, то β = γ > 0, иначе β = –1. Далее строится изоповерхность β = β0 > 0 (рис. 1б–1в). Если в некоторой фиксированной точке течения β > 0, то в декартовой СК x с началом в этой точке и двигающейся со скоростью этой точки можно записать обыкновенное дифференциальное уравнение ${\mathbf{v}} = {{d{\mathbf{x}}} \mathord{\left/ {\vphantom {{d{\mathbf{x}}} {dt}}} \right. \kern-0em} {dt}} \approx {\mathbf{Gx}}$, где v – скорость движения жидкости в СК x. Можно показать [16], что фазовой траекторией частицы жидкости в СК x является плоская спираль, по которой она двигается вокруг выбранной точки течения с угловой скоростью β. Пространственные вихревые структуры двухнитевого следа и цепочки шпилькообразных вихревых петель в следе за сферой, обтекаемой однородной вязкой жидкостью, полученные при помощи β-визуализации в [15, 16], хорошо согласуются с вихревыми структурами, полученными в экспериментах [18, 19] при помощи визуализации краской или эмульсией.

Картина изолиний β > 0 в вертикальной плоскости X–Z на рис. 1г подобна картине линий тока в СК2 на рис. 1к. Каждой полуволне на рис. 1ж, 1к можно поставить в соответствие полуволну на рис. 1г. Таким образом, на рис. 1г в плоскости X–Z можно выделить три полуволны впадин (соответствующие темным полосам 1, 3, 5 с S < 0 на рис. 1ж) и две полуволны гребней (соответствующие светлым полосам 2, 4).

Из рис. 1г следует переход от двухмерной вихревой структуры внутренних волн в плоскости X–Z к пространственной структуре волн на рис. 1б–1в, где показаны U-образные структуры первых полуволн впадин и гребней, часть второй полуволны впадин, а также V-образные структуры “осевых частей” первой (только на рис. 1в), второй и третьей полуволн гребней около оси Z, которые связаны между собой и с D1 при помощи горизонтальных вихревых нитей.

3. ЭВОЛЮЦИЯ ТЕЧЕНИЯ ЖИДКОСТИ ПРИ T ≤ 0.7

При старте диска в горизонтальном направлении в устойчиво стратифицированной вязкой жидкости с плотностью ${{\rho }_{{00}}}(X) = 1 - \frac{X}{{2{\text{A}}}}$, частички жидкости около диска выйдут из состояния покоя и начнут колебаться с частотой плавучести N в вертикальном направлении. Со временем эти колебания в ближнем следе затухнут, сформировав внутренние волны (рис. 1), распространяющиеся в СК2 справа налево вместе с телом и со скоростью тела. Ниже проводится детализация процесса формирования этих волн.

Для Fr = 4, Re = 50 в неподвижной системе координат СК1 в картине линий тока при T ≤ 0.001 наблюдается ламинарное обтекание диска; при T = 0.008 – отрыв потока около тыльных кромок диска и присоединение потока к тыльной стороне диска выше оси Z. Здесь работает основной механизм формирования течения вязкой жидкости 1k (генерация вихревого кольца (или полукольца) у поверхности обтекаемого тела) [16]. Цифра 1 в названии механизма 1k обозначает, что этот механизм функционирует в D1.

На рис. 2 при Т = 0.04, 0.08, 0.12, 0.14, 0.32 за диском наблюдаются D1, длиной L/d = 0.388, 0.565, 0.653, 0.671, 0.729. При Т > 0.32 L/d больше не растет. При этом толщина D1 у тыльной стороны диска равна 0.928d. Для Fr = 4, Re = 50 вихревая оболочка (область течения D2 на рис. 2г) и D1 осесимметричны.

Рис. 2.

Течение за диском при Fr = 4, Re = 50, А = 2776.4: (а–д) – изолинии β > 0 при φ = π/4 с шагами 0.2, 0.1, 0.0002, 0.0005, 0.002 при T = 0.02, 0.04, 0.06, 0.08, 0.19; (е–и) – изоповерхности β = 0.003 при T = 0.08, 0.16, 0.24, 0.32

Если q – вертикальная прямая, проходящая через место Q импульсного старта центра тыльной стороны диска, то в СК2 этот центр за время T переместился влево на расстояние $s = 4{\pi Fr}T$ от неподвижной прямой q. При 0.02 < Т < 0.04 между правым краем оболочки и прямой q на месте резкого сдвига жидкости в окрестностях плоскостей φ = π/4 и φ = 3π/4 начинают формироваться четыре вихревые структуры, симметричные относительно плоскостей X–Z и Y–Z (правая половина рис. 2б). Дальнейшее развитие этих структур при 0.06 ≤ Т ≤ 0.32 показано на рис. 2 как при помощи изолиний β > 0 в плоскости φ = π/4 (в правых половинах рис. 2в–2д), так и при помощи изоповерхностей β (рис. 2е–2и). При Т ≥ 0.08 эти четыре структуры уже похожи на вихревые нити. На рис. 2б–2ж при 0.04 ≤ Т ≤ 0.19 прямая q проходит через правые края этой первой четверки нитей. При 0.14 < T < 0.16 левые концы первой четверки нитей (около D1) индуцируют вторую четверку вихревых нитей (рис. 2д, 2ж). (Цифры 1 и 2 на рис. 2ж указывают на первую и вторую четверку вихрей соответственно).

При 0.22 < T < 0.24 правее каждой пары правых концов первой и второй четверки нитей индуцируется по одному головному вихрю (рис. 2з). При 0.3 < T < 0.5 каждая пара нитей с их головным вихрем превращаются в шпилькообразную вихревую петлю (рис. 2и и 3а). Таким образом, первая и вторая четверки вихревых нитей превращаются в первую и вторую пары вихревых петель; ножки второй пары петель стыкуются с оболочкой D2. Из рис. 4г следует, что вращение жидкости в сечении головной части первой вихревой петли вертикальной плоскостью X–Z (ячейка –1) при X > 0 идет по часовой стрелке. Следовательно, в вертикальном сечении ножки первого вихря при X > 0, Y > 0 и Z > 0 вращение жидкости идет против часовой стрелки (если смотреть на тыльную сторону диска).

Рис. 3.

Течение за диском при Fr = 4, Re = 50, А = 2776.4, T = 0.72: (а–б) – изоповерхности β = 0.003, 0.001; (в–е) – изолинии Sz с шагом 10–7 (в), β с шагом 0.001 (г) и S с шагом 10–6 (е) и мгновенные линии тока в СК2 (д) в вертикальной плоскости X–Z

Рис. 4.

Течение за диском при Fr = 4, Re = 50, А = 2776.4 в вертикальной плоскости X–Z: (а–г) – мгновенные линии тока в СК2 (I), изолинии Sz × 106 с шагами 10, 0.1, 0.1, 1 (II) и S ⋅ 106 с шагами 3, 0.1, 0.1, 5 (III) при T = 0.01, 0.19, 0.22, 0.24

Согласно [16] детальный механизм формирования вихрей при T < 0.5 можно символически записать как M1 = {2k–1k–3f–3f–3t/b–3r/l}, где 2k – формирование оболочки диска, цифрой 3 обозначается течение D3 вне оболочки (рис. 2г), 3f – генерация в D3 четырех вихревых нитей (сначала первой четверки, а потом второй), 3t/b/r/l – генерация в D3 головных частей вихревых петель, ориентированных вверх (t), вниз (b), вправо (r) и влево (l). Формирование этих вихревых петель (механизм M1) обусловлено сдвиговой и гравитационной неустойчивостями стратифицированной вязкой жидкости, инициируемыми движением диска.

Динамика картин мгновенных линий тока в СК2 и изолиний S и Sz в вертикальной плоскости X–Z при 0.01 ≤ T ≤ 0.24 на рис. 4 дополняет описанный выше процесс формирования вихревых петель (пространственный механизм M1). Начало движения диска формирует при X > 0 одну большую циркуляционную ячейку 0 в картине мгновенных линий тока в подвижной системе координат СК2, две ячейки в картине изолиний S (левую и правую) и три ячейки в картине изолиний Sz (рис. 4а). В картинах изолиний Sz (рис. 4, II) прямая q (черная вертикальная линия) проходит немного правее границы между центральной и правой ячейками, т.е. центральная ячейка визуализирует путь, пройденный телом в СК2. В картинах изолиний S (рис. 4, III) при 0.01 < T < 0.24 правая ячейка тоже сдвинулась вправо. По картинам изолиний S и Sz невозможно понять механизм эволюции течения, который наглядно показывает мгновенные линии тока в плоскости X–Z. Вызванное движением диска интенсивное вращение жидкости по часовой стрелке левее прямой q приводит к формированию циркуляционной ячейки 1 правее q (с вращением жидкости на рис. 4в против часовой стрелки), которая в свою очередь индуцирует (совместно с правыми концами первой четверки вихрей на рис. 2з) формирование циркуляционной ячейки –1 (головной части первой вихревой петли) около оси Z правее q (рис. 4г).

При 0.25 ≤ T ≤ 0.7 диск вместе с нулевой ячейкой сдвигается влево, освобождая место для приближения ячейки 1 к оси Z (рис. 3д), т.е. левее прямой q происходит формирование первой полуволны впадин (рис. 3е) из ячейки 1.

При 0.4 ≤ T ≤ 0.7 на границе циркуляционных ячеек 1 (с циркуляцией жидкости против часовой стрелки) и –1 (с циркуляцией жидкости по часовой стрелке) более тяжелая жидкость поднимается вдоль прямой q на равновесный уровень для более легкой жидкости. Поэтому отрицательные значения S около прямой q становятся положительными, и в картине изолиний S вдоль прямой q при 0.4 ≤ T ≤ 0.7 опускается некоторая новая светлая ячейка с S > 0 (рис. 3е и 5а–б (III)), которая делит правую ячейку с S < 0 на рис. 4г (III) на ячейку 1 (левее q) и ячейку –1 (правее q). Таким образом, ячейки 1 и –1 в картинах изолиний S формируются позднее, чем циркуляционные ячейки 1 и –1 в картинах мгновенных линий тока, и формирование первой и второй пары вихревых петель при 0.4 ≤ T ≤ 0.7 сопровождается формированием первой полуволны впадин.

Тот же самый пространственный механизм M1 работает и при Fr = 0.5 (рис. 5а–5б и 6а–6д). Правда, при одном и том же T путь $s = 4{\pi Fr}T$, который прошел диск в СК2, для Fr = 0.5 будет в 8 раз короче пути для Fr = 4. Можно было бы предположить, что картины течения для Fr = 0.5 получаются из картин для Fr = 4 путем простого сжатия их в 8 раз по горизонтали. Но в реальности это не совсем так. Например, сильное вращение жидкости в циркуляционной ячейке –1 (головной части первой петли) (рис. 5а–5б (II)) при 0.25 ≤ T ≤ 0.6 генерирует левее прямой q высокую базу 1 первой полуволны впадин, а циркуляционной ячейки 1 над базой 1 на рис. 5а–5б уже не видно. А при 0.01 < T < 0.1 для Fr = 0.5 первая четверка вихрей зарождается впритык к вихревой оболочке (рис. 6а–6б). Поэтому вторая четверка формируется при 0.14 < T < 0.16 уже около боковых сторон диска (рис. 6в).

Рис. 5.

Течение за диском при Fr = 0.5, Re = 50, А = 981.6 в плоскости X–Z: а–з – изолинии β с шагом 0.002 (I), мгновенные линии тока (II), изолинии S × 106 с шагами 2, 2, 6, 6, 6, 3, 3, 3 (III) при T = 0.35, 0.6, 0.7, 0.75, 0.8, 0.9, 1.2, 1.25

Рис. 6.

Течение за диском при Fr = 0.5, Re = 50, А = 981.6: (а) – изолинии β > 0 в плоскости φ = π/4 при T = 0.1 с шагом 0.1; (б–и) – изоповерхности β = 0.005 при T = 0.1, 0.2, 0.25, 0.7, 0.8, 0.9, 1.3, 1.5

При 0.24 < T < 0.25 для Fr = 0.5 также начинают формироваться головные части перед первой четверкой вихрей (рис. 6г). Головные части второй четверки вихрей для Fr = 0.5 не видны на рис. 6г. Таким образом, для 0.5 ≤ Fr ≤ 4 этапы формирования первой пары вихревых петель при T ≤ 0.7 практически слабо зависят от Fr, а определяются только моментом времени T.

4. ЭВОЛЮЦИЯ ТЕЧЕНИЯ ЖИДКОСТИ ПРИ T > 0.7

При 0.25 < T < 0.8 вокруг прямой q сформировалась головная часть первой петли (рис. 6г–6д), включающая в себя циркуляционные ячейки 1 и –1 на рис. 5а–5д (II). Будем считать ее первым вихревым полукольцом. При 0.6 < T ≤ 0.7 на границе ячеек 1 и –1 в окрестности прямой q наблюдается безвихревое движение жидкости вверх (рис. 5б–5в (II)). При 0.7 < T ≤ 0.9 за счет сдвиговой и гравитационной неустойчивостей вокруг прямой q (рис. 5г–5е) формируется деформированное вихревое кольцо 2 (рис. 6е–6ж), которое при 0.9 < T ≤ 1.2 сдвигается вниз (ближе к точке Q), и его левая половина превращается в первую полуволну гребней (рис. 6з) (циркуляционную ячейку 2 на рис. 5е–5ж). А правая половинка кольца 2 (соответствующая циркуляционной ячейке –2 в плоскости XZ на рис. 5ж) остается около прямой q. Скорость вращения жидкости в циркуляционной ячейке 2 на рис. 5ж много больше, чем в ячейке –2, поэтому правой половинки кольца 2 не видно на рис. 6з. При 1.2 < T ≤ 1.7 формируется и сдвигается к оси Z кольцо 3 (рис. 5ж–5з и 6з–6и), создавая вторую полуволну впадин (рис. 7а–7б). Таким образом, при T > 0.7 в течение каждого $\Delta T = 0.5$, над точкой Q наблюдается процесс зарождения вихревых колец (пространственный механизм M2 = {3k}), левые половины которых трансформируются в полуволны, заполняющие пространство между диском и точкой Q. При этом правые половины колец утончаются со временем под давлением новорожденных правых полуколец, давящих на них сверху. Универсальный механизм M2 работает и в случае покоящегося диска.

Рис. 7.

Течение за диском при Fr = 0.5, Re = 50, А = 981.6 в плоскости X–Z: (а–з) – изолинии S ⋅ 106 с шагами 2, 2, 2, 1, 1, 0.5, 0.5, 5 при T = 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3

Рисунок 5ж–5з (II) наглядно подробно демонстрируют, как реализуется гравитационная неустойчивость при 1.2 < T ≤ 1.25. Безвихревой поток около прямой q на границе между ячейками 2 и –2 при 1.1 < T ≤ 1.15 направлен вниз. При T = 1.2 скорость вдоль правой границы q2 циркуляционной ячейки 2, наклоненной под некоторым небольшим углом к вертикали, становится немного больше, чем у левой границы ячейки –2, из-за того, что ячейка 2 расположена ближе к двигающемуся в СК2 диску. Более легкая жидкость опускается вдоль прямой q2 на уровень более тяжелой жидкости. А силы плавучести стремятся вернуть более легкую жидкость вверх на ее уровень. В результате при T = 1.25 линия тока q2 становится волнообразной и около нее в вертикальной плоскости XZ зарождаются два небольших вихря (рис. 5з (II)), которые в пространстве являются двумя частями нового вихревого кольца 3 (рис. 6и) в сечении его плоскостью XZ. В результате при T = 1.25 на “пустом” месте генерируется кольцо 3.

При T > 0.4 в картине изолиний S в вертикальной плоскости X–Z на рис. 5 (III) и на рис. 7 периодически в течение $\Delta T = 0.5$ вдоль прямой q опускается некоторая новая ячейка, которая делит ячейку под ней на две ячейки: левую (новую полуволну) и правую. На рис. 7з–7ж при T = 3 хорошо видны ячейки (полуволны) 1–5 (левее q), более тонкие ячейки –1 – –5 (правее q) и наполовину разделенная ячейка (через которую проходит прямая q), из которой при T = 3.2 сформируются ячейки 6 и –6. Картины изолиний S при Fr = 0.5, Re = 50 для T = 3 и 11.46 (рис. 7з и 1ж) очень похожи, т.е. пространственная вихревая структура внутренних волн при T = 3 будет примерно такой же, как на рис. 1б–1в.

Обратимся к процессу формирования течения при Fr = 4, Re = 50 для T > 0.72. На рис. 3а, 3г видно начало формирования кольца 2 при T = 0.72. На рис. 3б (вверху и внизу) вихревая структура первой полуволны впадин похожа на сплюснутую вихревую петлю, ножки которой параллельны ножкам второй пары вихревых петель. При T > 0.995 путь, пройденный телом, равен s = 4π ∙ FrT > 50 = zmax, т.е. в СК1 прямая q вместе с головной частью первой петли выйдет за границы расчетной области. Поэтому в установившемся при T = 1.2 течении для Fr = 4 и Re = 50 наблюдаются только левые половины ножек первой пары вихревых петель (рис. 8а, справа). А почти все пространство между диском и правой границей расчетной области при T = 1.2 занимает полуволна впадин 1 (рис. 8б–8е). На картине изолиний возмущения давления p в плоскости X–Z (в центре рис. 8е) можно увидеть линию впадин, как и на картине изолиний Sz на рис. 8г. Установившееся при Fr = 4 и Re = 50 течение около тыльной стороны диска является квази-осесимметричным (рис. 8з).

Рис. 8.

Установившееся со временем течение за диском при Fr = 4.0, Re = 50, А = 2776.4, Т = 1.2 (L/d = 0.729, Сd = 1.89) в пространстве (а), в плоскостях X–Z (б, г–е) и φ = π/4 (в) и на диске (ж–з): (а) – изоповерхность β = 0.003; б – линии тока в СК2; (в–е) – изолинии β > 0 с шагом 0.002 (в), Sz с шагом 5 × 10–8 (г), S с шагом 10–6 (д) и p с шагом 10–3 (е); ж–з – линии трения на передней (ж) и тыльной (з) поверхностях диска в СК1

ЗАКЛЮЧЕНИЕ

В результате численного решения системы уравнений Навье–Стокса в приближении Буссинеска и визуализации пространственной вихревой структуры рассчитанного течения впервые детально рассмотрен процесс формирования трехмерных гравитационных внутренних волн над местом Q импульсного старта центра тыльной стороны диска с диаметром d и толщиной h = 0.76d в горизонтальном направлении вдоль оси симметрии диска Z справа налево в линейно стратифицированной по плотности вязкой жидкости при Fr = 0.5 и 4 для Re = 50.

Рассчитанные поля векторов скоростей и возмущения давления обладают горизонтальной и вертикальной плоскостями симметрии, проходящими через ось Z. Поэтому процесс формирования течения, обусловленный сдвиговой и гравитационной неустойчивостями, описывается в верхнем полупространстве следующим образом. Сначала между тыльным торцом диска и местом Q импульсного старта этого торца формируются две горизонтальные вихревые нити, которые потом трансформируются в ножки шпилькообразной вихревой петли, головная часть которой расположена правее Q. Далее над точкой Q периодически в течение промежутка времени, равного половине периода плавучести жидкости, формируются вихревые кольца, левые половины которых трансформируются в полуволны, заполняющие пространство между диском и точкой Q. При этом правые половины колец утончаются со временем под давлением новорожденных правых полуколец, давящих на них сверху.

Приведены детальное описание и результаты анализа динамики формирования пространственных вихревых структур в линейно стратифицированной вязкой сплошной среде, создаваемые движущимися в горизонтальном направлении объектами в форме диска.

Список литературы

  1. Lighthill J. Waves in Fluids. Cambridge: CUP, 1978. = Лайтхилл Дж. Волны в жидкостях. М.: Мир, 1981. 598 с.

  2. Миткин В.В., Чашечкин Ю.Д. Трансформация висящих разрывов в вихревые системы в стратифицированном течении за цилиндром // Изв. РАН. МЖГ. 2007. № 1. С. 15–28.

  3. Матюшин П.В. Эволюция течения стратифицированной вязкой жидкости при начале движения тела // Процессы в геосредах. 2016. № 4 (9). С. 333–343.

  4. Чашечкин Ю.Д., Воейков И.В. Вихревые системы за цилиндром в непрерывно стратифицированной жидкости // Изв. РАН. Физика атмосферы и океана. 1993. Т. 29. № 6. С. 821–830.

  5. Boyer D.L., Davies P.A., Fernando H.J.S., Zhang X. Linearly Stratified Flow Past a Horizontal Circular Cylinder // Philosophical Transactions of the Royal Society of London. Series A: Mathematical and Physical Sciences. 1989. V. 328. № 1601. P. 501–528.

  6. Lin Q., Lindberg W.R., Boyer D.L., Fernando H.J.S. Stratified flow past a sphere // J. Fluid Mech. 1992. V. 240. P. 315–354.

  7. Chomaz J.M., Bonneton P., Hopfinger E.J. The structure of the near wake of a sphere moving horizontally in a stratified fluid // J. Fluid Mech. 1993. V. 254. P. 1–21.

  8. Boyer D.L., Davies P.A. Laboratory studies of orographic effects in rotating and stratified flows // Ann. Rev. Fluid Mech. 2000. V. 32. P. 165–202.

  9. Gushchin V.A., Matyushin P.V. Mathematical Modeling of the Incompressible Fluid Flows // AIP Conf. Proc. 2014. V. 1631. P. 122–134.

  10. Гущин В.А., Матюшин П.В. Моделирование и исследование течений стратифицированной жидкости около тел конечных размеров // Журн. вычисл. матем. и матем. физики. 2016. Т. 56. № 6. С. 1049–1063.

  11. Матюшин П.В. Классификация режимов течений стратифицированной вязкой жидкости около диска // Процессы в геосредах. 2017. № 4 (13). С. 678–687.

  12. Hanazaki H. A numerical study of three-dimensional stratified flow past a sphere // J. Fluid Mech. 1988. V. 192. P. 393–419.

  13. Байдулов В.Г., Матюшин П.В., Чашечкин Ю.Д. Эволюция течения, индуцированного диффузией на сфере, погруженной в непрерывно стратифицированную жидкость // Изв. РАН. МЖГ. 2007. № 2. С. 130–143.

  14. Белоцерковский О.М., Гущин В.А., Коньшин В.Н. Метод расщепления для исследования течений стратифицированной жидкости со свободной поверхностью // ЖВМ и МФ. 1987. Т. 27. № 4. С. 594–609.

  15. Матюшин П.В. Численное моделирование пространственных отрывных течений однородной несжимаемой вязкой жидкости около сферы. Дис. … канд. физ.-мат. наук: 05.13.18. М., 2003. 194 с.

  16. Гущин В.А., Матюшин П.В. Механизмы формирования вихрей в следе за сферой при 200 < Re < 380 // Изв. РАН. МЖГ. 2006. № 5. С. 135–151.

  17. Bobinski T., Goujon-Durand S., Wesfreid J.E. Instabilities in the wake of a circular disk // Phys. Rev. 2014. V. E 89. P. 053021.

  18. Magarvey R.H., Bishop R.L. Transition ranges for three-dimensional wakes // Can. J. Phys. 1961. V. 39. P. 1418–1422.

  19. Sakamoto H., Haniu H. A study on vortex shedding from spheres in a uniform flow // Trans. ASME: J. Fluids Engng. 1990. V. 112. P. 386–392.

Дополнительные материалы отсутствуют.