Известия РАН. Механика жидкости и газа, 2023, № 4, стр. 117-130

1. ПОСТАНОВКА ЗАДАЧИ

Диск с диаметром d и толщиной h = 0.76d вносится в покоящуюся линейно стратифицированную по плотности несжимаемую вязкую жидкость и сразу же начинает равномерное движение со скоростью U вдоль горизонтальной оси симметрии Z диска справа налево. Пусть точка Q старта центра тыльной стороны диска будет началом покоящейся декартовой системы координат СК2 (z, x, y,), где ось x – вертикальна, ось z совпадает с осью Z. Пусть начало двигающейся декартовой системы координат СК1 (Z, X, Y) с вертикальной осью X связано с геометрическим центром О диска. Для решения поставленной задачи моделируется вспомогательная задача обтекания диска равномерным потоком жидкости со скоростью U в СК1. Направление потока жидкости совпадает с положительным направлением оси Z.

Пусть T_b – период плавучести жидкости, T – обезразмеренное на T_b время, прошедшее с момента внесения диска в жидкость. Плотность жидкости ${{\rho }}(Z,X,Y) = 1{\kern 1pt} - {\kern 1pt} 0.5 \cdot {X \mathord{\left/ {\vphantom {X {\text{A}}}} \right. \kern-0em} {\text{A}}} + S(Z,X,Y)$ обезразмерена на плотность ρ₀ на уровне центра диска, а координаты Z, X, Y – на d/2, где ${\text{A}} = \Lambda {\text{/}}d$ – отношение масштабов, N_b= 2π/T_b и $\Lambda = {\text{g/}}N_{{\text{b}}}^{2}$ – частота и масштаб плавучести жидкости, g – ускорение свободного падения, S(Z, X, Y) – обезразмеренное на ρ₀ возмущение солености, которое при T = 0 равно нулю.

Для математического моделирования вспомогательной задачи обтекания диска в СК1 решается следующая обезразмеренная система уравнений Навье–Стокса в приближении Буссинеска [12]:

Она записана в цилиндрической системе координат (Z, R, φ): Z = Z, $X = R\cos \varphi $, $Y = R\sin \varphi $, где ${\mathbf{v}} = \left( {{{{\text{v}}}_{Z}},{{{\text{v}}}_{X}},{{{\text{v}}}_{Y}}} \right)$ – вектор скорости, нормированный по U; p – возмущение давления, обезразмеренное на ${{\rho }_{0}}U_{{}}^{2}$; t – время (обезразмеренное на $f = d{\text{/}}\left( {2U} \right) = 1{\text{/}}\left( {2{\text{Fr}} \cdot {{N}_{{\text{b}}}}} \right)$); Re = Ud/ν – число Рейнольдса, Fr = UT_b/(2πd) – внутреннее число Фруда, Sc = ν/κ = 709.22 – число Шмидта; ν = = 0.01 см²/c – коэффициент кинематической вязкости воды, κ = 4.1 × 10^–5 см²/c – коэффициент диффузии соли; ∇ и Δ – операторы Гамильтона и Лапласа.

Для решения этой задачи используется численный метод расщепления по физическим факторам МЕРАНЖ [13], который успешно применялся для моделирования течений несжимаемой вязкой жидкости около сфер, цилиндров и дисков [4, 6–8, 11, 14–16], а также для течений со свободной поверхностью [13]. Детали численного метода, цилиндрической расчетной сетки, граничных условий и результатов тестирования программного комплекса математического моделирования и визуализации 3D-течений стратифицированной вязкой жидкости около диска опубликованы в [11, 14]. Приведенная в [14] классификация течений стратифицированной вязкой жидкости около диска толщиной h = 0.76d при 0.05 < Fr < 100 и 50 < Re < 500 хорошо согласуется c трехмерными расчетами [15, 16] и c экспериментами [17, 18].

Внешняя правая граница расчетной области удалена от центра диска на 25d. При Fr ≥ 10 длина внутренней волны λ = 2πFrd ≥ 20πd ≈ 62.83d, т.е. на выбранной расчетной сетке уместится менее чем 40% длины первой волны. Более того, амплитуда этой волны будет мала. Поэтому характер течения около диска будет эквивалентен течению однородной вязкой жидкости [6]. При Fr = = 4 ‒ λ = 8πd ≈ 25.13d, т.е. только одна волна будет доступна для наблюдения [11]. Поэтому для рассмотренной расчетной области именно при Fr ≤ 4 возможно исследовать генерацию внутренних волн, которые занимают все пространство расчетной сетки левее точки Q [11] (рис. 1а). Расчеты проводились на вычислительных ресурсах Межведомственного суперкомпьютерного центра Российской академии наук (МСЦ РАН).

В СК2 значения безразмерных горизонтальных компонент v_Z векторов скорости, рассчитанные в СК1, уменьшаются на единицу (v_z = v_Z – 1), а значения переменных v_x, v_y, S и p в СК2 те же, что и в СК1. Таким образом, в СК2 скорость набегающего на диск потока равна нулю, а безразмерная скорость диска равна –1, т.е. диск равномерно двигается справа налево в покоящейся жидкости. Центр тыльной стороны диска за время $T = \left[ {tf} \right]{\text{/}}{{T}_{b}} = \left[ {t{\text{/}}\left( {2{\text{Fr}}{{N}_{{\text{b}}}}} \right)} \right]\left[ {{{N}_{{\text{b}}}}{\text{/}}2\pi } \right] = t{\text{/}}\left( {4{{\pi Fr}}} \right)$ смещается в покоящейся СК2 на расстояние $L = U\left[ {tf} \right]{\text{/}}\left( {0.5d} \right) = t = 4{{\pi Fr}}T$ влево от точки Q (рис. 1е).

2. ВИЗУАЛИЗАЦИЯ ВИХРЕВОЙ СТРУКТУРЫ ТЕЧЕНИЯ ЖИДКОСТИ

Из-за стратификации жидкости, медленной скорости диска (Re = 50) и симметричности диска относительно плоскости x–z рассчитанные 3D-поля векторов скорости в любой момент времени характеризуются горизонтальной y–z и вертикальной x–z плоскостями симметрии. Следовательно, имеют смысл картины мгновенных линий тока в плоскостях x–z (рис. 1е–1ж) и y–z (рис. 1з). На рис. 1е для СК2 визуализируются вихревые ячейки 1, 2, 3, соответствующие внутренним полуволнам 1, 2, 3, сгенерированным при T = 1.5 для Fr = 0.3. На рис. 1ж–1з для СК1 показана рециркуляционная зона R, сформированная в результате отрыва набегающего на диск потока от тыльной стороны диска.

Известно, что построение 3D мгновенных линий тока с целью визуализации 3D вихревой структуры течения – трудное и неблагодарное дело. Можно легко потерять многие вихри. Поэтому для визуализации этой 3D-структуры естественно было бы использовать изоповерхности модуля завихренности ω = |ω| = |rot v| (рис. 1а). Рассмотрим сначала картины изолиний фитой компоненты завихренности ω_φ в вертикальной плоскости x–z (рис. 1в). Из определения ω_φ следует, что в ячейках 1, 3 и –2 на рис. 1в, где ω_φ > 0, жидкость вращается против часовой стрелки, что подтверждается рисунком 1е. В ячейках 0, 2, –3, –1 на рис. 1в, где ω_φ < 0, жидкость вращается по часовой стрелке (см. рис. 1е). Изоповерхность ω = 0.01 на рис. 1а визуализирует внутренние полуволны 1–3 (области течения со значительным вращением сплошной среды), но не способна показать структуру более слабых вихрей около оси z [4]. Для визуализации всей вихревой структуры 3D течения несжимаемой жидкости в статье [2] было предложено использовать изоповерхность λ₂ < 0 (рис. 2б). Здесь λ₂ – второе собственное значение симметричного тензора ${{{\mathbf{F}}}^{2}} + {{{\mathbf{B}}}^{2}}$, где ${\mathbf{F}} + {\mathbf{B}} = {\mathbf{G}}$ – тензор градиента скорости c элементами ${{{\text{v}}}_{{i{\text{,}}j}}} = {{\partial {{{\text{v}}}_{i}}} \mathord{\left/ {\vphantom {{\partial {{{\text{v}}}_{i}}} {\partial {{x}_{j}}}}} \right. \kern-0em} {\partial {{x}_{j}}}}$, ${{F}_{{i,j}}} = 0.5\left( {{{{\text{v}}}_{{i,j}}} + {{{\text{v}}}_{{j,i}}}} \right)$, B_{i, j} = = $0.5({{{\text{v}}}_{{i,j}}} - {{v}_{{j,i}}})$. В [2] также упоминается визуализация 3D-течений при помощи изоповерхности мнимой части β двух комплексно-сопряженных собственных значений тензора G (рис. 1б, 1и). Интересно отметить, что в [2] эта мнимая часть никак не обозначается. Поэтому в разных статьях ее обозначают при помощи разных символов: Im(σ_1,2) [6, 15], β [8, 11, 14, 16], λ_ci [19], λ_cx [20] и т.п. Сначала использовалась λ₂-визуализация [4], а потом β-визуализация [6, 8, 11, 14–16], которая имеет следующий ясный физический смысл. Выберем в СК1 произвольную точку М со скоростью жидкости v_M в ней. Перейдем в систему отсчета СК3(v_M), двигающуюся относительно СК1 со скоростью v_M. Тогда скорость в точке M в СК3(v_M) станет равной нулю, а для скоростей в малой окрестности этой точки будет справедливо обыкновенное дифференциальное уравнение (ОДУ) ${\mathbf{v}} = {{d{\mathbf{x}}} \mathord{\left/ {\vphantom {{d{\mathbf{x}}} {dt}}} \right. \kern-0em} {dt}} \approx {\mathbf{Gx}}$. Известно, что наличие двух комплексно-сопряженных собственных значений ${{\sigma }} = {{\alpha }} + i{{\beta }}$ и ${{\bar {\sigma }}} = {{\alpha }} - i{{\beta }}$ тензора G в точке М дает следующее решение этого ОДУ [21]:

Здесь величины С и С₃ – комплексная и действительная константы соответственно, σ₃ – третье (действительное) собственное значение тензора G, h₃ – третий (действительный) собственный вектор, ${\mathbf{h}} = {{{\mathbf{h}}}_{1}} - i{{{\mathbf{h}}}_{2}}$ и ${\mathbf{\bar {h}}}$ – два комплексно-сопряженных собственных вектора тензора G, h₁ и h₂ – два линейно независимых действительных вектора, которые являются базисом в плоскости P. Пусть комплексное число ζ = ξ₁ + iξ₂ = Сexp(σt), тогда мы имеем: x_I = ξ₁ ∙ h₁ + ξ₂ ∙ h₂. Отобразим аффинно плоскость P на вспомогательную плоскость P* комплексного переменного ζ, чтобы вектор h₁ перешел в единицу, а вектор h₂ – в i. Тогда вектору x_I = ξ₁ ∙ h₁ + ξ₂ ∙ h₂ будет соответствовать комплексное число ζ = ξ₁ + iξ₂. В силу этого отображения траектория в плоскости P перейдет в траекторию в плоскости P*, описываемую уравнением ζ = Сexp(σt). В полярных координатах (r, θ) ζ = rexp(iθ), т.е. ξ₁ = rcosθ, ξ₂ = rsinθ. Пусть С = С_realexp(iD), тогда r = С_realexp(αt), θ = βt + D. Таким образом, если β ≠ 0, то можно говорить о вихревом характере движения жидкости в малой окрестности точки М в СК3(v_M): молекулы жидкости движутся по овалам (при α = 0) или овальным спиралям (при α ≠ 0) и β – это усредненная по времени угловая скорость вращения молекул жидкости около точки М в СК3(v_M) [21]. Для удобства построения одноцветной изоповерхности β = β₀ > 0 (рис. 1и), демонстрирующей вихревую структуру течения жидкости, в тех ячейках расчетной сетки, где β = 0, функция β переопределяется: β = –0.01 [6, 8, 11, 14–16].

В [22] была введена двухцветная β+-визуализация, которая при помощи знака ω_φ окрашивает полуволны впадин 1, 3 и гребней 2 на рис. 1б разными цветами, что очень удобно для исследования вихревой структуры внутренних волн (рис. 3–5, IV). Функция β+ определена при β > 0: β+ = sign(ω_φ)β, где функция sign(ω_φ) = 1 при ω_φ ≥ 0, sign(ω_φ) = –1 при ω_φ < 0. Аналогичным образом можно ввести двухцветную ω+-визуализацию: ω+ = sign(ω_φ)ω (рис. 1а).

Рис. 3.

Течение около диска при Fr = 0.3, Re = 50, T_b = 2π c: а–в – мгновенные линии тока в СК2 (I) и изолинии 10³β+ с шагами 500, 0.1, 1 (II) в вертикальной плоскости x–z, изолинии β+ с шагами 0.5, 0.1, 0.1 в плоскости φ = π/4 (III), изоповерхности 10³β = 0.5, 5 и 10³β+ = ±0.5 (IV) при T = 0.1, 0.2, 0.24.

Рис. 4.

Течение около диска при Fr = 0.3, Re = 50, T_b = 2π c: а–в – мгновенные линии тока в СК2 (I) и изолинии 10²β+ с шагами 0.2, 5, 5 (II) в плоскости x–z, изолинии β+ с шагами 0.1, 0.08, 0.04 в плоскости φ = π/4 (III), изоповерхности β+ = ±0.0052, β = 0.15, β+ = ±0.05 (IV) при T = 0.28, 0.45, 0.68.

Рис. 5.

Течение около диска при Fr = 0.3, Re = 50, T_b = 2π c: а–в – мгновенные линии тока в СК2 (I) и изолинии β+ с шагом 10^–2 (II) в плоскости x–z, изолинии β+ с шагом 0.08 в плоскости φ = π/4 (III), изоповерхности (β–) = = 5 × 10^–3, β+ = ±10^–2, ±10^–2 (IV) при T = 0.8, 0.86, 1.

Далее было замечено, что на бóльшей площади рис. 1в ω_φ < 0. Поэтому, если на рис. 1в убрать вихревые структуры, для которых ω_φ ≥ 0, то структура течения станет в два раза проще, но при этом удаленные структуры можно легко мысленно восстановить. Поэтому в [22] была введена одноцветная (β–)-визуализация течения жидкости. Функция (β–) определена при β > 0 и ω_φ < 0: (β–) = β (рис. 2а). В 3D-случае при сравнении изоповерхностей β+ на рис. 1б и изоповерхности β– на рис. 2а видно, что на рис. 2а все полуволны впадин, кроме первой, бесследно исчезают. На рис. 1б черной линией отмечена примерная граница между областями (0 + S) и 1, где S – это головная часть боковой вихревой петли (рис. 2е) [11]. Удаление области (S+1) оставляет грубый след на оболочке следа 0 (сравните рис. 4в, IV и 5а, IV).

Рис. 1г–1д демонстрируют физический смысл β для двух точек М₁ и М₂ со скоростями в СК1 (v_Z, v_R, v_φ) = (0.973, 0.005, 0) (при x > d/2) и (0, –0.01, 0), соответственно, отмеченных черными крестами. Рис. 1д разделен черной прямой x = d/2 на верхнюю и нижнюю части. Картины линий тока в верхней и нижней частях рис. 1д показаны в СК3 (0.973, 0.005, 0) и СК3(0, –0.01, 0), соответственно, и демонстрируют вращение жидкости в окрестностях точек М₁ и М₂. Картины линий тока в верхней и нижней частях рис. 1д похожи на картины линий тока в СК2 (рис. 1е) и СК1 (рис. 1ж) соответственно. Так, на рис. 1д и 1ж при 0 ≤ x ≤ d/2 мы видим зону R. А на рис. 1д и 1е при x > d/2 мы видим полуволны впадин и гребней. Поэтому для приблизительного описания изменения кинематики внутренних волн и отрывных течений около оси z в вертикальной плоскости можно использовать мгновенные линии тока в СК2 и СК1 соответственно. Так, на рис. 1, 3–6 для каждой картины изолиний β+ в плоскости x–z приводится соответствующая ей картина мгновенных линий тока в СК2.

Рис. 6.

Течение около диска при Fr = 0.3, Re = 50, T_b = 2π c: а–д – мгновенные линии тока в СК2 (I) и изолинии 10³β+ с шагом 2 (II) в плоскости x–z, изоповерхности β+ = ±10^–3 (III) при T = 1.25, 1.3, 1.4, 1.73, 1.9.

По аналогии с β+ и (β–)-визуализациями можно определить λ₂+ и (λ₂–)-визуализации (рис. 2б). Функция λ₂+ определена при λ₂ < 0: λ₂+ = sign(ω_φ)|λ₂|. Функция (λ₂–) определена при λ₂ < 0 и ω_φ < 0: (λ₂–) = λ₂. Сравнение рис. 2а и 2б для T = 1.5 и Fr = 0.8 показывает, что топологии вихревых структур, полученные при помощи (β–) и (λ₂–)-визуализаций, практически совпадают, но (β–)-визуализация более гармонично представляет вихревую структуру течения жидкости. Так, на рис. 2а нити f₁ соединяются с кольцом R, а на рис. 2б не соединяются. На рис. 2а вихревое кольцо –1 цельное, а на рис. 2б – разорвано на три части.

3. МЕХАНИЗМ ФОРМИРОВАНИЯ ВНУТРЕННИХ ВОЛН ПРИ Fr ≤ 4

При T = 0 диск вносится в покоящуюся стратифицированную вязкую жидкость и сразу же начинает равномерное движение вдоль горизонтальной оси z симметрии диска справа налево. На потревоженную жидкость, отклонившуюся от своего начального положения, начинают действовать силы плавучести (последний член уравнения (1.2)). С другой стороны, в результате прилипания молекул жидкости к поверхности диска, перед диском и по бокам диска на жидкость действуют силы, направленные справа налево. Эти силы приводят к генерации вихревого тора с осью симметрии z, заполняющего все пространство (рис. 3а, I при T = 0.1), а силы плавучести при этом генерируют гравитационные внутренние полуволны над Q (рис. 3в, I при T = 0.24). На рис. 3–5, IV показана динамика 3D вихревой структуры течения жидкости при 0 < T < 1, обусловленная гравитационной и сдвиговой неустойчивостями. Гравитационная неустойчивость, в свою очередь, обусловлена силами плавучести жидкости. Каждую 3D-структуру на рис. 3–5 дополняют изолинии β+ в ее сечениях при φ = 0 и φ = π/4 и мгновенные линии тока в СК2 в вертикальной плоскости x–z (при φ = 0). На рис. 3а, IV одна четверть изоповерхности β = 5 × 10^–4 вырезана, чтобы показать поверхность диска.

При T ≤ 0.1 и Fr ≤ 4 вне осесимметричной вихревой оболочки 0 диска у тыловой острой кромки диска при x > 0 появляются две небольшие вихревые структуры (рис. 3а, III–IV и рис. 2а–2г, 2е, 6а–6б из [11]) – зачатки двух вихревых нитей f₁ (рис. 3б–3в, III–IV). Картина изолиний β+ в сечении φ = π/4 демонстрирует динамику роста нитей f₁ и рециркуляционной зоны R. При 0.2 ≤ ≤ T ≤ 0.24 гравитационная и сдвиговая неустойчивости приводят к тому, что форма мгновенных линий тока над точкой Q в СК2 при φ = 0 из прямолинейной становится волнообразной (рис. 3в, I). Почему это происходит, подробно описано в [11]. В результате над Q формируются два каскада вихрей: 1 (слева) и –1 (справа) (рис. 3в, II). При T = 0.28 для Fr = 0.3 сердцевины каскадов 1 и –1 образуют кольцо 1 (рис. 4а, IV).

При Fr ≤ 1 и Fr ≥ 2 детали процесса формирования вихрей немного отличаются. Если при 0.1 < < T ≤ 0.28, Fr ≥ 2 и x > 0 оболочка 0 следа практически не меняется со временем, а две первые нити f₁ и две вторые (боковые) нити f_S формируются вне оболочки 0 (рис. 2е) [11], то при 0.1 < T ≤ 0.28, Fr ≤ 1 и x > 0 нити f₁ и f_S размещаются в самой оболочке 0. Поэтому при Fr ≤ 1 имеет смысл разделить оболочку 0 на внутреннюю 0₁ и внешнюю 0₂ части при помощи изоповерхности β₀. Модуль β₀ приведен на рис. 3–5, III. Оболочка 0₁ слабо меняется со временем и состоит из двух колец r_I, r_II, которые генерируются двумя острыми кромками диска (рис. 3а, IV). На рис. 3б, III в сечении оболочки 0₂ при φ = π/4 можно выделить нити f₁, заштрихованные в клеточку, зону R у тыльной стороны диска и остальную часть 0₂ над боковой стороной диска. При x > 0 и 0.2 ≤ T ≤ 1 две вихревые нити f₁ (ножки вихревой петли –1) являются частью каскада –1, а петля –1 – сердцевиной каскада –1. Генерация двух симметричных петель –1 при x > 0 и x < 0 сопровождается зарождением боковых петель S, как боковых частей двух симметричных каскадов 1. При x > 0 и 0.24 ≤ T ≤ 0.32 для Fr = 0.3 наблюдаются четыре боковые нити f_S (рис. 3в, 4а, IV), для Fr = 4 – две нити f_S (рис. 2е). Появление петель S связано с присутствием диска и не упоминается в обсуждаемом здесь МФВВ, так как этот МФВВ посвящен непосредственно генерации внутренних волн.

Таким образом, при 0.1 ≤ T ≤ 0.32 формируются две вихревые петли –1 с ножками f₁ и боковые петли S. При T = 0.28 для Fr = 0.3 сердцевины каскадов 1 и –1 образуют кольцо 1 (рис. 4а, IV), а для Fr = 4 полукольцо 1 похоже на арку (как на рис. 3б в статье [11]), радиус которой примерно в пять раз больше радиуса полукольца –1. Поэтому в формулировке МФВВ надо говорить о формировании не одного вихревого кольца 1, а двух каскадов вихрей 1 и –1, сердцевинами которых будут полукольца 1 и –1.

При 0.28 < T ≤ 0.45 и Fr = 0.3 диаметр D(1) сечения полукольца 1 плоскостью x–z уменьшается в три раза, а далее исчезает совсем. D(1) = 0 говорит о том, что полукольцо 1 при x > 0 разделилось на две половинки, симметричные относительно плоскости x–z. При x > 0 и 0.28 ≤ T ≤ 0.8 диск сдвигается влево, а каскад –1 остается у точки Q, поэтому ножки f₁ петли –1 вытягивается в горизонтальном направлении (рис. 4б–4в, 5а, IV); полукольцо –1 сплющивается в вертикальном направлении и превращается в полукруг –1. На рис. 4б и 4в, IV показаны только половины изоповерхностей β для того, чтобы показать диск. При 0.5 ≤ T ≤ 0.8 над Q из-за гравитационной и сдвиговой неустойчивостей генерируются каскады вихрей 2 и –2. Каскад вихрей 2 со временем станет полуволной гребней 2. Вихревые нити f₂, обведенные черными кругами на рис. 5а, IV для Fr = 0.3 или представленные в форме полукольцевой перемычки между нитями f₁ на рис. 3а из [11] для Fr = 4, – это нижние части каскада вихрей 2. Таким образом, при Fr = 4 полуволны гребней формируются из четных нитей около Q, а не из полуколец над Q как при Fr = 0.3.

Пусть β_max(k) – значение локального максимума функции β в сечении полукольца k плоскостью x–z. При T = 0.8 и Fr = 0.3 сердцевина каскада 2 приближается к оси z, β_max(2) = 0.025, β_max(–1) = = 0.142, т.е. полукруг –1 обладает наибольшим β_max в сечении φ = 0 вне оболочки диска (рис. 5а, II). При x > 0 и T = 0.86 каскад 2 садится на левый край полукруга –1, β_max(2) увеличивается до 0.088, значения β перераспределяются между левым и правым краями полукруга –1: β_max(–1, левый) = = 0.12, β_max(–1, правый) = 0.068 (рис. 5б, II). Теперь в СК2 через нити f₁ прокачивается жидкость, поступившая как из полукруга –1, так и из каскада 2, что в дальнейшем при T = 1 приводит к трансформации полукруга –1 в кольцо –1 (рис. 5в, II, β_max(2) = 0.22, β_max(–1) = 0.038). На рис. 5б и 5в, IV изоповерхности β показаны при –π/2 ≤ φ ≤ π/10 для того, чтобы показать диск. Таким образом, в результате слияния при 0.86 ≤ T ≤ 1 сердцевины каскада 2 и левой части полукруга –1 каскад 2 стал обладать наибольшим β_max в сечении φ = 0 вне оболочки диска.

Итак, при T ≤ 1 над Q периодически в течение каждой ∆T = 0.5 формируются каскады вихрей k и –k, где k = 1, 2. Сердцевины каскадов –1 и 1 состоят из вихревой петли –1 (с ножками f₁ и полукольцом –1) и из полукольца 1 с боковыми петлями S соответственно. Сначала каскад 2 состоит из нитей f₂ (около Q) и полукольца 2, которое потом садится на левую часть полукольца –1. В результате формируются полуволна впадин 1 (над нитями f₁) и полуволна гребней 2. Фундаментом вихревой структуры первой волны при T = 1 является горизонтальное кольцо –1, сформировавшееся из полукольца –1.

Описанный выше процесс формирования вихрей для T ≤ 1 схематически представлен на рис. 2в–2г. Его можно обобщить на любой другой интервал времени (n–1) < T ≤ n, где n = 2, 3, 4, … . Действительно, при T > 0 над Q периодически в течение каждого интервала времени 0.5(k–1) < < T  ≤ 0.5k формируются левый k и правый –k каскады вихрей. Они состоят из нитей f_k около Q и деформированных полуколец k и –k, где k = 1, 2, 3, … . Для каждого нечетного k у оси z формируется вихревая петля –k. Она состоит из нитей f_k и нижнего полукольца –k, на которое потом садится четный каскад (k + 1). Далее нижнее полукольцо –k становится деформированным кольцом. Таким образом, в течение каждого интервала времени (n–1) < T ≤ n формируется внутренняя волна n, состоящая из полуволн впадин k = (2n–1) и гребня (k + 1) = 2n. Левая часть кольца –k становится осевой частью гребня (k + 1). В результате осевые части гребней оказываются связанными друг с другом нечетными нитями в цепочку. МФВВ для Fr ≤ 4 и ∆T = 1 можно сформулировать короче. Нечетные каскады k и –k, нити f_k + нижний вихрь –k; четные каскады (k + 1) и –(k + 1), нити (k + 1), каскад (k + 1) садится на нижний вихрь –k (рис. 2в–2д).

На рис. 2ж приведен график зависимости коэффициента C_d лобового сопротивления диска от времени t при Fr = 0.3. Волнообразная форма графика при T ≤ 1 обусловлена тем, что при T ≤ 1 точка Q, над которой работает МФВВ, находится около диска. Реализация МФВВ приводит к возмущению полей векторов скорости и давления у поверхности диска, которые и определяют C_d.

Если при T ≤ 4 и Fr = 0.3 реализации МФВВ для разных временных интервалов ∆T = 1 немного отличаются друг от друга, то при T > 3 и Fr = 0.3 похожи друг на друга. Ниже подробно описываются реализации МФВВ при 1 < T ≤ 3, Fr = 0.3, T_b = 2π c и Re = 50, опираясь на рис. 1, 6–7. Описание особенностей МФВВ при T > 3 приведено без помощи рисунков.

Рис. 7.

Течение около диска при Fr = 0.3, Re = 50, T_b = 2π c: а–г – изолинии 10³β+ с шагом 2 (I) в плоскости x–z, изоповерхности (β–) = 10^–3 (II) при T = 2.05, 2.6, 2.8, 2.95.

4. ДЕТАЛИЗАЦИЯ ПРОЦЕССА ГЕНЕРАЦИИ ВОЛН ПРИ Fr = 0.3 И T > 1

Рассмотрим теперь более детально МФВВ при 1 < T ≤ 2 (рис. 6, 2д, 1). Так как при T = 1 и Fr = 4 точка Q достигла внешней границы расчетной сетки, то дальнейшее исследование МФВВ при Fr = 4 стало невозможным. При 1 < T ≤ 1.3 и Fr = 0.3 над точкой Q генерируются каскады вихрей 3 и –3 (рис. 6а–6б). При 1 < T ≤ 1.3 β_max(2) увеличивается до 0.33. При T > 1.3 β_max(2) больше не меняется, что говорит о том, что при T = 1.3 гребень 2 и вихревая структура левее его уже примерно сформировались. На рис. 6, III и 7, II изоповерхности β показаны при π/2 ≤ φ ≤ π для того, чтобы показать диск и для лучшего понимания 3D вихревой структуры течения жидкости. При T = 1.25 на рис. 6а, III широкий “амфитеатр” полуволны 2 стоит перед меньшей по размерам “сценой” полукольца –2, боковые части которого выделены черной окружностью. Эти боковые части индуцируют появление нитей f₃ над точкой Q при T = 1.3. Нити f₃, выделенные черным прямоугольником на рис. 6б–6в, III, являются частью формирующегося каскада –3. При T = 1.4 полукольцо –2 отрывается от полуволны 2, β_max(3) = 0.021, β_max(–3) = 0.013, размеры нитей f₃ увеличиваются (рис. 6в, III). Нити f₃ расположены выше нитей f₁. Картина течения при T = 1.5, представленная на рис. 1а–1з и 2а–2б, уже была подробно описана в подразделе 2 настоящей статьи. Здесь к этому описанию остается добавить устоявшуюся внутреннюю структуру (“скелет”) первой волны на рис. 1и, где боковые нити f_S соединены с кольцом r_I, а нити f₁ – с кольцом r_II. Здесь так же видно, что при T = 1.5 сердцевина полуволны 2 – это деформированное кольцо 2 (рис. 1к). На рис. 6г, III при T = 1.73 нити f₃ соединены с кольцом 2. Таким образом, в стратифицированной жидкости каждое вихревое кольцо связано с четырьмя нитями. В то время как в однородной жидкости вихревое кольцо связано только с двумя нитями [6].

При 1.5 < T ≤ 1.73 над точкой Q генерируются каскады вихрей 4 и –4 (рис. 6г). При 1.3 < T ≤ 1.73 β_max(–2) уменьшается до 0.003. При 1.4 < T ≤ 1.73 β_max(3) и β_max(–3) увеличиваются до 0.027 и 0.052, соответственно (рис. 6г, II), полукольцо –3 превращается в полукруг –3, из нитей f₃ и полукруга –3 формируется вихревая петля –3 (рис. 6г, III); D(3) уменьшается примерно в полтора раза (аналогично уменьшению D(1) при 0.28 < T ≤ 0.45), но при этом объем полуволны впадин 3 значительно увеличивается. При 1.73 < T ≤ 1.9 β_max(3) уменьшается до 0.023, β_max(–3) увеличивается до 0.057 (рис. 6д, II). При T = 1.9 β_max(4) = 0.0124, β_max(–4) = 0.0056.

При 1.9 < T ≤ 2.05 развиваются две полуволны 4 и –4. Вертикальная полуволна 4 садится на левую часть вытянутого горизонтального полукруга –3, который затем трансформируется в деформированное кольцо –3, выделенное на рис. 7а, II двумя черными кругами с черной ломаной кривой между ними; β_max(–3) = 0.0046, β_max(4) и β_max(–4) увеличиваются до 0.08 и 0.0092 соответственно. Таким образом, при 1 < T ≤ 2 формируются полуволна впадин 3 над нитями f₃ и полуволна гребней 4. Фундаментом вихревой структуры второй волны при T = 2 является горизонтальное кольцо –3. При 2.05 < T ≤ 2.95 диаметры деформированных колец –1 и –3 постепенно увеличиваются (рис. 7, II).

При 2.05 < T ≤ 2.5 над точкой Q генерируются каскады вихрей 5 и –5. Если каскады –1 и –3 во время своего зарождения при T = 0.28 и 1.3, соответственно, имели одну ярко выраженную сердцевину, то на рис. 7б–7г при 2.6 ≤ T ≤ 2.95 наблюдаются сразу три сердцевины каскада –5 в форме полуколец. При 2.1 ≤ T ≤ 2.5 боковые части полукольца –4 формируют нити f₅ (аналогично формированию нитей f₃), которые стыкуются с нижним полукольцом каскада –5, которое при T = 2.6 соприкасается с верхней частью левого внутреннего куска –3* правой части кольца –3 (рис. 7б). В силу того, что кусок –3* расположен ближе к оси z, то β_max(–3) = 0.011 > 0.005 = β_max(–5). Поэтому при 2.6 < T ≤ 2.8 кусок –3* поглощает полукольцо –5, т.е. –3* := (–3*) + (–5). В результате нити f₅ соединяются с куском –3* (рис. 7в). При 2.5 < T ≤ 2.95 над точкой Q генерируются каскады вихрей 6 и –6. При 2.8 < T ≤ 2.95 кусок –3* отрывается от правой части кольца –3 и трансформируется в сильно деформированное кольцо –3* ≡ –5, боковая часть которого выделена черным прямоугольником на рис. 7г, II. Далее вертикальный каскад 6 садится на левую часть кольца –5, а полуволна 1, состоящая ранее из двух половинок, снова становится полноценным полукольцом, β_max(1) = 0.008 (рис. 7г, I). Таким образом, при 2 < T ≤ 3 формируются полуволна впадин 5 над нитями f₅ и полуволна гребней 6. Фундаментом вихревой структуры третьей волны является кольцо –5. При T ≥ 2.8 β_max(4) = 0.14 < 0.33 = β_max(2), т.е. при T = 2.8 гребень 4 и вихревая структура левее его уже примерно сформировались. Поэтому ширина полуволн 2-4 на рис. 7в–7г, измеряемая вдоль прямой, перпендикулярной им, приблизительно равна πFrd ≈ 0.94d. Это хорошо согласуется с линейной теорией внутренних волн.

При T = 2.95 вихревые кольца –1 и –3 вместе с полукольцами –5 и –6 создают вихревое “гнездо” (рис. 7г). В нем потом генерируются каскады вихрей 7 и –7 с несколькими сердцевинами, далее у оси z появляется вихревая петля –7, затем формируются вихревые каскады 8 и –8 и каскад 8 садится на головную часть петли –7. При T = 4 и π/2 ≤ φ ≤ π боковая часть кольца –5 становится дополнительной нитью, связывающей осевые части гребней –5 и –7 наряду с нитью f₇. Для числа Fr = 0.3 все МФВВ при T > 4 похожи на МФВВ при 3 < T ≤ 4. В результате осевые части гребней оказываются связанными друг с другом в цепочку нечетными нитями как “позвонки” “скелета” вихревой структуры течения (рис. 7г, II). На каждый такой “позвонок” прикреплены две полуволны (впадин и гребня). По мере уменьшения Fr горизонтальные размеры между центрами “позвонков”, равные 2πFrd, уменьшаются.

ЗАКЛЮЧЕНИЕ

Проведено математическое моделирование равномерного движения диска с диаметром d и толщиной h = 0.76d в горизонтальном направлении вдоль его оси симметрии Z в покоящейся линейно стратифицированной по плотности несжимаемой вязкой жидкости при Re = 50, T_b = 2π с и в широком диапазоне Fr. Диск генерирует пространственные (3D) гравитационные внутренние волны, занимающие все пространство между диском и местом Q его старта. Существенно дополнен опубликованный в [11] МФВ, работающий при формировании гравитационных внутренних волн, полученный на основе анализа динамики 3D вихревой структуры течения жидкости, визуализируемой при помощи изоповерхности мнимой части β > 0 комплексно-сопряженных собственных значений тензора градиента скорости. Для нахождения β в каждом узле расчетной сетки при помощи численного метода МЕРАНЖ [13] решалась система уравнений Навье–Стокса в приближении Буссинеска, записанная в цилиндрической системе координат. Для большей наглядности полуволны гребней и впадин были окрашены разными цветами в зависимости от знака угловой компоненты завихренности ω_φ на них. Такая визуализация получила наименование β+ [22].

В [11] основное внимание уделяется периодическому процессу зарождения деформированных вихревых колец над точкой Q, происходящему в силу гравитационной и сдвиговой неустойчивостей. Тогда левое полукольцо трансформируется в полуволну впадин или гребней, а правое полукольцо исчезает со временем. В настоящей работе подчеркивается важная роль нечетных правых полуколец в МФВ и сформулирован следующий универсальный МФВВ за диском для Re = 50, T_b = 2π c и Fr ≤ 4 в верхнем полупространстве над осью z. При T > 0 периодически в течение каждой ∆T = 0.5 над точкой Q формируются левый k и правый –k каскады вихрей, состоящие из нитей f_k около Q и деформированных полуколец k и –k, где k = 1, 2, 3, … (см. рис. 3б–4а и 5а–5б для k = 1 и 2 соответственно). Для каждого нечетного k у оси z формируется вихревая петля –k, состоящая из нитей f_k и нижнего полукольца –k, на которое потом садится четный каскад (k + 1). Далее головная часть петли –k становится деформированным кольцом. Таким образом, в течение каждой ∆T = 1 формируется новая внутренняя волна, состоящая из полуволн впадин k и гребня (k + 1). Левая часть кольца –k становится осевой частью гребня (k + 1). В результате осевые части гребней оказываются связанными друг с другом в цепочку нечетными нитями. Чтобы увидеть эту вихревую цепочку, нужно визуализировать только ту часть изоповерхности β = β₀, на которой ω_φ < 0. Такая визуализация получила наименование (β–) [22]. Для разных интервалов ∆T = 1 МФВВ имеет свои особенности, подробно описанные в настоящей работе. Существует небольшая зависимость МФВВ от Fr. Например, при 0.8 < Fr ≤ 4 полуволны гребней формируются из четных нитей около Q, а не из полуколец над Q как при Fr ≤ 0.5. При Fr = 0.8 и T ≥ 1 все МФВВ для разных временных интервалов ∆T = 1 похожи на МФВВ при Fr = = 0.3 и T ≥ 3.

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

Расчеты проводились на вычислительных ресурсах Межведомственного суперкомпьютерного центра Российской академии наук (МСЦ РАН).