Instability of solid-body rotation of heton type

Full Text

1. ВВЕДЕНИЕ

Одна из центральных проблем геофизической гидродинамики – проблема гидродинамической неустойчивости атмосферных и океанических сдвиговых течений. С развитием неустойчивости связывают образование вихрей различных масштабов – от масштаба отдельного облака до масштаба циклона. Проблеме гидродинамической неустойчивости посвящены многочисленные обзоры и монографии [Pedlosky, 1987; Vallis 2013; Dolzhansky, 2013; Badin, 2013; Sokolovskiy and Verron, 2014; Kalashnik et. al., 2021, 2022]. Широко известны и монографии на русском языке [Дикий, 1976; Шакина, 1990], монография и статьи [Булатов и др., 2015; 2017].

В настоящей работе рассмотрена одна из малоизученных задач теории устойчивости – задача о неустойчивости осесимметричных вихрей с твердотельным вращением. Хорошо известно, что в однородной по плотности жидкости такие баротропные вихри устойчивы. Как показано в настоящей работе, ситуация принципиально меняется для бароклинных вихрей с различным направлением вращения на верхней и нижней границе слоя атмосферы. Такие вихри являются неустойчивыми. В работе получено аналитическое выражение для инкремента нарастания неустойчивых возмущений, исследована его зависимость от внешних параметров. Показано, что для неустойчивости в модели учитывающей экмановское трение, горизонтальное волновое число возмущения должно быть меньше некоторого критического значения. Этот факт приводит к существованию кривой нейтральной устойчивости, разделяющей волновые числа растущих и нейтральных возмущений.

Следует отметить, что бароклинные вихри с различным направлением вращения на двух уровнях в литературе принято называть хетонами. Этот термин был впервые предложен в работе [Hogg and Stommel, 1985], для подчеркивания возможности бароклинных вихрей переносить тепло. Обширный список публикаций, посвященных динамике хетонов, представлен в обзоре [Гряник и др., 2003] и монографии [Sokolovskiy and Verron, 2014].

2. ДВУХУРОВЕННАЯ ДИСКРЕТНАЯ КВАЗИГЕОСТРОФИЧЕСКАЯ МОДЕЛЬ И ПОСТАНОВКА ЗАДАЧ УСТОЙЧИВОСТИ ДЛЯ ОСЕСИММЕТРИЧНЫХ ТЕЧЕНИЙ

Для изучения устойчивости вихрей используем дискретный вариант поверхностной геострофической модели (Surface quasigeostrophic (SQG) model), описывающей движения слоя стратифицированной вращающейся жидкости толщины H с нулевой потенциальной завихренностью [Blumen, 1978; Held et al., 1995]. Толщина слоя H принимается в качестве вертикального масштаба. В безразмерных переменных уравнения SQG модели включают уравнение Лапласа для функции тока, к которому присоединяются уравнения переноса плавучести на границах. Используется дискретная аппроксимация оператора Лапласа, включающая разбиение вертикального отрезка интегрирования на четыре отрезка толщины h = 1 / 4 и привлечение уравнений переноса плавучести на границах. Система уравнений модели имеет вид

$$\begin{gathered} q_{1t}+\left[{\psi }_1,q_1\right]=\mathrm{0,} \\ {q}_{2t}+\left[{\psi }_2,q_2\right]=\mathrm{0,} \\ {q}_{\mathrm{1,2}}=\Delta {\psi }_{\mathrm{1,2}}\pm ({\lambda }^2/2)({\psi }_2-{\psi }_1). \end{gathered}$$ (1)

Здесь ψ₁, ψ₂ – значения функции тока на верхней и нижней границах слоя, Δ – двумерный оператор Лапласа, $\left[a,b\right]=a_x{b}_y-a_y{b}_x$ – двумерный якобиан, параметр ${\lambda }^2=1/2h^2=8$. Частные производные обозначены нижними буквенными индексами. В качестве горизонтального масштаба в безразмерных уравнениях модели принят бароклинный радиус деформации Россби D = NH / f, где N – частота Брента, f – параметр Кориолиса. Масштаб времени, функции тока и скорости соответственно $T_{\ast }=D/U_{\ast }$, ${\Psi }_{\ast }=U_{\ast }D$ и $U_{\ast }={\Theta }_{\ast }/N$, где ${\Theta }_{\ast }$ – характерное значение перепада плавучести на границах. Безразмерные горизонтальные компоненты скорости u, v связаны с функцией тока ψ соотношениями u = –ψ_y, v = –ψ_x. Плавучесть $\vartheta$ находится из дискретной аппроксимации соотношения $\vartheta ={\psi }_z$.

Уравнения (1) совпадают с уравнениями классической двухуровенной модели Филлипса [Pedlosky, 1987; Phillips, 1954], однако отличаются от нее физической интерпретацией переменных. Распределения PV q₁, – q₂ на нижнем и верхнем уровнях в модели Филлипса, теперь имеют смысл удвоенных распределений плавучести на верхней и нижней границах слоя атмосферы. Подробный вывод уравнений (1) представлен в работах [Kalashnik, 2020; Kalashnik et al., 2021, 2022].

Из уравнений (1) следует закон сохранения полной энергии

$E_t=\mathrm{0,} E=(1/2)\iint \left({(\nabla {\psi }_1)}^2+{(\nabla {\psi }_2)}^2+4{({\psi }_2-{\psi }_1)}^2\right)dx$, (2)

$x=(x,y)$. Сохраняются также квадратичные интегралы

${\left(V_{\mathrm{1,2}}\right)}_t=\mathrm{0,} V_{\mathrm{1,2}}=\iint q_{\mathrm{1,2}}^{2}dx$, (3)

отражающие законы сохранения энстрофии в непрерывной модели.

Отметим, что для характерных значений параметров тропосферы средних широт H = 10 км, N = 10^–2 (сек)^–1, f = 10^–2 (сек)^–4, радиус деформации D = NH / f = 1000 км. Характерному значению скорости $U_{\ast }=10 м/сек$ отвечает перепад плавучести ${\Theta }_{\ast }={10}^{-1}м{\left(сек\right)}^{-2}$ (перепад потенциальной температуры порядка) и временной масштаб $T_{\ast }=D/U_{\ast }=27.8 час$ (порядка суток).

Введением баротропного $\sigma =(1/2)({\psi }_1+{\psi }_2)$ и бароклинного компонентов $\theta =(1/2)({\psi }_2-{\psi }_1)$ функции тока, система (1) сводится к системе

$$\begin{gathered} {(\Delta \theta -{\lambda }^2\theta )}_t+\left[\sigma ,\Delta \theta -{\lambda }^2\theta \right]+\left[\theta ,\Delta \sigma \right]=\mathrm{0,} \\ {\left(\Delta \sigma \right)}_t+\left[\sigma ,\Delta \sigma \right]+\left[\theta ,\Delta \theta \right]=\mathrm{0,} \end{gathered}$$ (4)

удобной для анализа задач устойчивости.

Будем рассматривать систему (4) в цилиндрических координатах r, α, когда x = rcosα, y = rsinα. При этом якобиан и двумерный оператор Лапласа определяются выражениями $\left[\sigma ,\theta \right]=r^{-1}\left({\sigma }_r{\theta }_{\alpha }-{\sigma }_{\alpha }{\theta }_r\right)$, $\Delta \sigma ={\sigma }_{rr}+r^{-1}{\sigma }_r+r^{-2}{\sigma }_{\alpha \alpha }$. Не зависящая от угла часть оператора Лапласа определяется как $\Delta \sigma ={\sigma }_{rr}+r^{-1}{\sigma }_r$. Точное решение системы (4) $\theta =\overline{\theta }\left(r\right)$, $\sigma =\overline{\sigma }\left(r\right)$ отвечает осесимметричному течению с баротропным и бароклинным составляющими азимутальной компоненты скорости $\overline{V}=-{\overline{\sigma }}_r$, $\tilde{V}={\overline{\theta }}_r$. Полная азимутальная компонента скорости течения на двух уровнях определяется выражением $V_{\mathrm{1,2}}=\overline{V}\left(r\right)\mp \tilde{V}\left(r\right)$. Линеаризованная форма системы (4) в цилиндрических координатах имеет вид

$$\begin{gathered} {\left(\Delta {\theta }^{\text{'}}-{\lambda }^2{\theta }^{\text{'}}\right)}_t+r^{-1}{\left(\overline{\sigma }\right)}_r{\left(\Delta {\theta }^{\text{'}}-{\lambda }^2{\theta }^{\text{'}}\right)}_{\alpha }- \\ -r^{-1}{\left({\sigma }^{\text{'}}\right)}_{\alpha }{\left(\Delta \overline{\theta }-{\lambda }^2\overline{\theta }\right)}_r+r^{-1}{\left(\overline{\theta }\right)}_r{\left(\Delta {\sigma }^{\text{'}}\right)}_{\alpha }=0, \\ {\left(\Delta {\sigma }^{\text{'}}\right)}_t+r^{-1}{\left(\overline{\sigma }\right)}_r{\left(\Delta {\sigma }^{\text{'}}\right)}_{\alpha }-r^{-1}{\left({\sigma }^{\text{'}}\right)}_{\alpha }{\left(\Delta \overline{\sigma }\right)}_r+ \\ +r^{-1}{\left(\overline{\theta }\right)}_r{\left(\Delta {\theta }^{\text{'}}\right)}_{\alpha }-r^{-1}{\left({\theta }^{\text{'}}\right)}_{\alpha }{\left(\Delta \overline{\theta }\right)}_r=0, \end{gathered}$$ (5)

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

3. НЕУСТОЙЧИВОСТЬ БАРОКЛИННОГО ВИХРЯ С ТВЕРДОТЕЛЬНЫМ ВРАЩЕНИЕМ

Далее будем рассматривать случай $\overline{\theta }\left(r\right)=S{r}^2/2=0$, $\overline{\sigma }\left(r\right)=B{r}^2/2$, отвечающий ситуации с твердотельным вращением в неограниченной области 0 < r < ∞. Учитывая, что в этом случае $\Delta \overline{\theta }=2S$, $\Delta \overline{\theta }-{\lambda }^2\overline{\theta }=2S-{\lambda }^{2}S{r}^2/2$, ${(\Delta \overline{\theta }-{\lambda }^2\overline{\theta })}_r=-{\lambda }^{2}Sr$, ${\overline{\sigma }}_r=Br$, $\Delta \overline{\sigma }=2B$, из (5) получим систему

${\left(\Delta {\sigma }^{\text{'}}\right)}_t+B{\left(\Delta {\sigma }^{\text{'}}\right)}_{\alpha }+S{\left(\Delta {\theta }^{\text{'}}\right)}_{\alpha }=0$,

${(\Delta {\theta }^{\text{'}}-{\lambda }^2{\theta }^{\text{'}})}_t+B{(\Delta {\theta }^{\text{'}}-{\lambda }^2{\theta }^{\text{'}})}_{\alpha }+S{({\lambda }^2{\sigma }^{\text{'}}+\Delta {\sigma }^{\text{'}})}_{\alpha }=0$. (6)

Первое уравнение (6) можно сократить на оператор Лапласа. Возможность такого сокращения следует из общей теории гармонических функций – гармоническая функция достигает максимальных и минимальных значений на границе области. Поэтому если она равна нулю на границе, то она равна нулю тождественно. С учетом этого сокращения уравнения (6) примут вид

${{\sigma }^{\text{'}}}_t+B{{\sigma }^{\text{'}}}_{\alpha }+S{{\theta }^{\text{'}}}_{\alpha }=0$,

${(\Delta {\theta }^{\text{'}}-{\lambda }^2{\theta }^{\text{'}})}_t+B{(\Delta {\theta }^{\text{'}}-{\lambda }^2{\theta }^{\text{'}})}_{\alpha }+S{({\lambda }^2{\sigma }^{\text{'}}+\Delta {\sigma }^{\text{'}})}_{\alpha }=0$. (7)

Рассмотрим вначале случай баротропного течения, когда S = 0. При этом из (6) получим независимые и распадающиеся уравнения

${{\sigma }^{\text{'}}}_t+B{{\sigma }^{\text{'}}}_{\alpha }=0$,

${(\Delta {\theta }^{\text{'}}-{\lambda }^2{\theta }^{\text{'}})}_t+B{(\Delta {\theta }^{\text{'}}-{\lambda }^2{\theta }^{\text{'}})}_{\alpha }=0$. (8)

Будем искать экспоненциально нарастающие во времени и периодические по азимутальной координате решения вида ${\theta }^{\text{'}}=\cos n\alpha e^{\gamma t}\tilde{\theta }$, ${\sigma }^{\text{'}}=\cos n\alpha e^{\gamma t}\tilde{\sigma }$, где γ – инкремент нарастания. Из уравнений (8) для амплитуд, отмеченных волной, получим уравнения

$\gamma \tilde{\sigma }+Bn\tilde{\sigma }=0$,

$\gamma ({\Delta }_2\tilde{\theta }-{\lambda }^2\tilde{\theta })+Bn({\Delta }_2\tilde{\theta }-{\lambda }^2\tilde{\theta })=0$. (9)

Здесь обозначено ${\Delta }_2\tilde{\theta }={\tilde{\theta }}_{rr}+r^{-1}{\tilde{\theta }}_r-r^{-2}{n}^2\tilde{\theta }$ – двумерный оператор Лапласа для периодических по углу возмущений. Для инкремента нарастания из (9) получим выражение γ = –Bn. Поскольку γ < 0 все возмущения баротропного течения затухают. Здесь мы получили общеизвестный результат.

Рассмотрим теперь основной случай бароклинного течения, когда B = 0. При этом азимутальная компонента скорости $V_{\mathrm{1,2}}=\mp \tilde{V}\left(r\right)=\mp Sr$, т.е. мы имеем течение хетонного типа (с различным направлением вращения на двух уровнях). Для этого течения система (7) сводится к системе

${{\sigma }^{\text{'}}}_t+S{{\theta }^{\text{'}}}_{\alpha }=0$,

${(\Delta {\theta }^{\text{'}}-{\lambda }^2{\theta }^{\text{'}})}_t+S{({\lambda }^2{\sigma }^{\text{'}}+\Delta {\sigma }^{\text{'}})}_{\alpha }=0$, (10)

которую можно свести к одному уравнению второго порядка по времени относительно .

${(\Delta {\theta }^{\text{'}}-{\lambda }^2{\theta }^{\text{'}})}_{tt}-S^2{({\lambda }^2{\theta }^{\text{'}}+\Delta {\theta }^{\text{'}})}_{\alpha \alpha }=0$. (11)

Для отыскания решений полагаем ${\theta }^{\text{'}}=\cos n\alpha e^{\gamma t}{\tilde{\theta }}^{\text{'}}$. Тогда

${\gamma }^2({\Delta }_2\widetilde{{\theta }^{\text{'}}}-{\lambda }^2\widetilde{{\theta }^{\text{'}}})+S^2{n}^2({\lambda }^2\widetilde{{\theta }^{\text{'}}}+{\Delta }_2\widetilde{{\theta }^{\text{'}}})=0$, (12)

где, как и ранее, двумерный оператор Лапласа ${\Delta }_2\tilde{\theta }={\tilde{\theta }}_{rr}+r^{-1}{\tilde{\theta }}_r-r^{-2}{n}^2\tilde{\theta }$.

Будем искать решения уравнения (12) такие что ${\Delta }_2{\tilde{\theta }}^{\text{'}}=-\mu {\tilde{\theta }}^{\text{'}}$, или подробнее ${\widetilde{{\theta }^{\text{'}}}}_{rr}+r^{-1}{\widetilde{{\theta }^{\text{'}}}}_r-r^{-2}{n}^2\widetilde{{\theta }^{\text{'}}}=-\mu \widetilde{{\theta }^{\text{'}}}$. Замена переменной $x=\sqrt{\mu }r$ и $y\left(x\right)={\theta }^{\text{'}}\left(r\right)={\theta }^{\text{'}}\left(\frac{x}{\sqrt{\mu }}\right)$ приводит последнее уравнение к виду

$\frac{1}{x}\frac{d}{dx}\left(x\frac{dy}{dx}\right)+(1-\frac{n^2}{x^2})y=0$.

Решение этого уравнения представляется в терминах функций Бесселя y(x) = AJ_n(x). Для переменной $\widetilde{{\theta }^{\text{'}}}$ таким образом получим решение $\widetilde{{\theta }^{\text{'}}}=C{J}_n\left(\sqrt{\mu }r\right)$, где C – произвольная константа, определяемая по начальным данным. Полное представление для возмущения при C = 1, таким образом, имеет вид ${\theta }^{\text{'}}=\cos n\alpha e^{\gamma t}{J}_n\left(\sqrt{\mu }r\right)$. Величина $\sqrt{\mu }$ имеет смысл горизонтального (по радиусу) волнового числа.

С учетом приведенного решения, для инкремента нарастания возмущений из (12) следует уравнение ${\gamma }^2(-\mu -{\lambda }^2)+(S^2{\lambda }^2{n}^2-S^2{n}^2\mu )=0$, откуда

${\gamma }^2=\frac{S^2{n}^2({\lambda }^2-\mu )}{{\lambda }^2+\mu }$. (13)

Согласно (13), при μ < λ² существует положительное значение γ, отвечающее растущему возмущению. Таким образом, существует неустойчивость относительно длинноволновых (по радиусу) возмущений. При фиксированном значении n наибольшую скорость роста имеют возмущения с μ = 0. С ростом углового волнового числа n инкремент нарастания неограниченно растет. Нейтральной моде (с нулевым инкрементом) отвечает значение μ = λ².

Изолинии бароклинного компонента полной функции тока (среднее течение плюс возмущение) на верхнем уровне для двух моментов времени t = 2 (верх) и t = 3 (низ) представлены на рис. 1. Этот компонент определяется выражением $\theta =\overline{\theta }\left(r\right)+\epsilon {\theta }^{\text{'}}=S{r}^2/2+\epsilon {\theta }^{\text{'}}$, ${\theta }^{\text{'}}=\cos n\alpha e^{\gamma t}{J}_n\left(\sqrt{\mu }r\right)$, где ε – амплитудный параметр. Приняты значения S = 2, n = 2, μ = 2, ε = 0.001. Как видно, при потере устойчивости осесимметричного течения формируется система замкнутых вихревых ячеек.

 

Рис. 1. Изолинии бароклинного компонента полной функции тока на верхнем уровне для двух моментов времени. По вертикальной и горизонтальной оси отложены безразмерные координаты y и x

 

Используя уравнение (11), можно доказать неустойчивость течения с твердотельным вращением в ограниченной области 0 < r < r₀. В этой области, азимутальная и радиальная компоненты скорости в двух слоях выражаются через функцию тока соотношениями

$u_{\theta }=\frac{\partial {\psi }_{\mathrm{1,2}}}{\partial r}, u_r=-\frac{1}{r}\frac{\partial {\psi }_{\mathrm{1,2}}}{\partial \theta }$.

С учетом соотношений $\sigma =(1/2)({\psi }_1+{\psi }_2)$ и $\theta =(1/2)({\psi }_2-{\psi }_1)$, имеем ψ₁ = σ – θ, ψ₂ = σ + θ. Отсюда следует, что на твердой границе вихря условие отсутствия радиальной компоненты скорости для периодических по углу линейных возмущений сводится к условию $\theta \left|{}_{r=r_0}\right.=0$. Для экспоненциально растущих возмущений, при этом, автоматически выполняется условие $\sigma \left|{}_{r=r_0}\right.=0$, т.е. условие отсутствия радиальной компоненты скорости. Это прямо следует из первого уравнения (10). Как и ранее, линейные возмущения для переменной θ удовлетворяют уравнению (11). Полагая здесь ${\theta }^{\text{'}}=\cos n\alpha e^{\gamma t}{\tilde{\theta }}^{\text{'}}$, получим уравнение с краевым условием

$$\begin{gathered} {\gamma }^2({\Delta }_2\widetilde{{\theta }^{\text{'}}}-{\lambda }^2\widetilde{{\theta }^{\text{'}}})+S^2{n}^2({\lambda }^2\widetilde{{\theta }^{\text{'}}}+{\Delta }_2\widetilde{{\theta }^{\text{'}}})=\mathrm{0,} \\ \widetilde{{\theta }^{\text{'}}}\left|{}_{r=r_0}\right.=0. \end{gathered}$$

С учетом волновых решений этого уравнения и краевого условия, полное возмущение θ' представляется в виде ${\theta }^{\text{'}}=Ñ\cos n\alpha e^{\gamma t}{J}_n\left(\sqrt{\mu }r\right)$. Из краевого условия $J_n\left(\sqrt{\mu }{r}_0\right)=0$ при этом получим ${\sqrt{\mu }}_k{r}_0={\lambda }_k$, откуда ${\sqrt{\mu }}_k={\lambda }_k/r_0$, где λ_k – нули функции Бесселя порядка n. Величина $\sqrt{\mu }$, имеющая смысл горизонтального (по радиусу) волнового числа, теперь принимает дискретный набор значений. Инкремент нарастания определяется прежним выражением (13), где теперь μ = μ_k. Неустойчивость существует при условии μ_k < λ², которое наиболее легко удовлетворяется для медленно меняющихся по радиусу возмущений.

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

4. НЕУСТОЙЧИВОСТЬ ТВЕРДОТЕЛЬНОГО ВРАЩЕНИЯ В МОДЕЛИ С ЭКМАНОВСКИМ ТРЕНИЕМ

Рассмотрим теперь бароклинное течение в присутствие экмановского трения. Как и течение Колмогорова, это течение является точным решением уравнений динамики только в присутствии внешней силы. Уравнения двухуровенной модели с внешней силой и экмановским трением имеют вид

$q_{1t}+\left[{\psi }_1,q_1\right]=-R\Delta {\psi }_1+F_1$,

$q_{2t}+\left[{\psi }_2,q_2\right]=-R\Delta {\psi }_2+F_2$,

$q_{\mathrm{1,2}}=\Delta {\psi }_{\mathrm{1,2}}\pm ({\lambda }^2/2)({\psi }_2-{\psi }_1)$. (14)

Здесь F₁, F₂ – компоненты внешней силы, R – коэффициент трения, определенный выражением $R={\left(\frac{\nu f}{2{U_{\ast }}^2}\right)}^{1/2}\frac{D}{H}$, где ν – коэффициент вертикальной турбулентной вязкости. Подробный вывод уравнений в модели с трением представлен в монографии [Pedlosky, 1987]. Отметим, что коэффициент трения может быть представлен в виде $R=E^{1/2}/2Ro$, где $E={(h_E/H)}^2$ – число Экмана, $h_E={\left(2\nu /f\right)}^{1/2}$ – толщина экмановского пограничного слоя, $Ro=U_{\ast }/fD$ – число Россби для принятого масштаба скорости. Отметим также, что для приведенных значений параметров и h_E = 0.5 км, коэффициент трения R = 4. При U_* = 10 м/с, число Россби Ro = 0.1.

Введением баротропного $\sigma =(1/2)({\psi }_1+{\psi }_2)$ и бароклинного компонентов $\theta =(1/2)({\psi }_2-{\psi }_1)$ функции тока, система (14) теперь сводится к системе

${(\Delta \theta -{\lambda }^2\theta )}_t+\left[\sigma ,\Delta \theta -{\lambda }^2\theta \right]+\left[\theta ,\Delta \sigma \right]=-R\Delta \theta +\tilde{F}$,

${\left(\Delta \sigma \right)}_t+\left[\sigma ,\Delta \sigma \right]+\left[\theta ,\Delta \theta \right]=-R\Delta \sigma +\overline{F}$, (15)

где $\overline{F}=(1/2)(F_1+F_2)$ и $\tilde{F}=(1/2)(F_2-F_1)$. Далее будем рассматривать случай $\overline{F}=0$ или F₁ = –F₂. При этом $\tilde{F}=F_2$. Бароклинное твердотельное вращение $\overline{\theta }\left(r\right)=S{r}^2/2$, $\overline{\sigma }=0$ является точным решением уравнений (15), если $R\Delta \overline{\theta }=\tilde{F}$ или, с учетом $\Delta \overline{\theta }={\overline{\theta }}_{rr}+r^{-1}{\overline{\theta }}_r=2S$, при внешней силе $\tilde{F}=2SR$. Будем считать последнее равенство выполненным.

Полагаем $\theta =\overline{\theta }\left(r\right)+{\theta }^{\text{'}}$, $\sigma =\overline{\sigma }\left(r\right)+{\sigma }^{\text{'}}$. Линеаризованная система для возмущений в случае $\overline{F}=0$ примет вид

${{\sigma }^{\text{'}}}_t+S{{\theta }^{\text{'}}}_{\alpha }=-R{\sigma }^{\text{'}}$,

${(\Delta {\theta }^{\text{'}}-{\lambda }^2{\theta }^{\text{'}})}_t+S{({\lambda }^2{\sigma }^{\text{'}}+\Delta {\sigma }^{\text{'}})}_{\alpha }=-R\Delta {\theta }^{\text{'}}$. (16)

Система (16) получается из системы (10) добавлением в правые части диссипативных слагаемых.

Обозначая оператор $L=\partial /\partial t+R$, первое уравнение (16) можно записать как $L{\sigma }^{\text{'}}=-S{{\theta }^{\text{'}}}_{\alpha }$. Применяя оператор L ко второму уравнению, получим одно уравнение второго порядка $L{(\Delta {\theta }^{\text{'}}-{\lambda }^2{\theta }^{\text{'}})}_t-S^2{({\lambda }^2{\theta }^{\text{'}}+\Delta {\theta }^{\text{'}})}_{\alpha \alpha }=-RL\Delta {\theta }^{\text{'}}$, или, после арифметических преобразований, уравнение

$$\begin{gathered} {(\Delta {\theta }^{\text{'}}-{\lambda }^2{\theta }^{\text{'}})}_{tt}-S^2{({\lambda }^2{\theta }^{\text{'}}+\Delta {\theta }^{\text{'}})}_{\alpha \alpha }+ \\ +2R{\left(\Delta {\theta }^{\text{'}}\right)}_t-R{\lambda }^2{{\theta }^{\text{'}}}_t+R^2\Delta {\theta }^{\text{'}}=0. \end{gathered}$$ (17)

Аналогично предыдущему, мы полагаем ${\theta }^{\text{'}}=\cos n\alpha e^{\gamma t}\tilde{\theta }\left(r\right)$. Тогда с учетом $\Delta {\theta }^{\text{'}}=\cos n\alpha e^{\gamma t}{\Delta }_2\tilde{\theta }$, ${\Delta }_2\tilde{\theta }={\tilde{\theta }}_{rr}+r^{-1}{\tilde{\theta }}_r-r^{-2}{n}^2\tilde{\theta }$, для $\tilde{\theta }\left(r\right)$ получим уравнение

$\begin{array}{c}{\gamma }^2({\Delta }_2\tilde{\theta }-{\lambda }^2\tilde{\theta })+S^2{n}^2({\lambda }^2\tilde{\theta }+{\Delta }_2\tilde{\theta })+\\ +R\gamma (2{\Delta }_2\tilde{\theta }-{\lambda }^2\tilde{\theta })+R^2{\Delta }_2\tilde{\theta }=0.\end{array}$ (18)

Как и ранее, мы рассматриваем возмущения, для которых ${\Delta }_2\tilde{\theta }=-\mu \tilde{\theta }$ или подробнее ${\tilde{\theta }}_{rr}+r^{-1}{\tilde{\theta }}_r-r^{-2}{n}^2\tilde{\theta }+\mu \tilde{\theta }=0$. При этом $\tilde{\theta }=C{J}_n\left(\sqrt{\mu }r\right)$, где J_n(z) — функция Бесселя порядка n и C – произвольная постоянная интегрирования. Для инкремента нарастания возмущений из (18) получим квадратное уравнение

$\begin{array}{c}{\gamma }^2(\mu +{\lambda }^2)+\gamma R(2\mu +{\lambda }^2)-\\ -S^2{n}^2({\lambda }^2-\mu )+R^2\mu =0.\end{array}$ (19)

Возникновению неустойчивости отвечает появления корня этого уравнения с положительной действительной частью. Легко показать, что неустойчивость существует, если для корней уравнения γ₁, γ₂ выполнено условие γ₁γ₂ < 0. С учетом (18) это условие можно записать в виде ${\gamma }_1{\gamma }_2=-S^2{n}^2({\lambda }^2-\mu )+R^2\mu <0$. Уравнение γ₁γ₂ = 0, таким образом, определяет кривую нейтральной устойчивости, связывающей радиальное и азимутальное волновые числа m, n. Это уравнение можно записать в виде

$\mu =\frac{S^2{n}^2{\lambda }^2}{S^2{n}^2+R^2}$. (20)

График нейтральной кривой (20) для значений параметров ${\lambda }^2=1/2h^2=8$, S = 2, R = 1, представлен на рис. 2. Существование нейтральной кривой является основной особенностью, вносимой экмановским трением. Нейтральная кривая появляется только в присутствии трения. При R = 0 из (20) следует соотношение μ = λ², которое дает параметры нейтральной моды (13) в модели без трения.

 

Рис. 2. Кривая нейтральной устойчивости на плоскости волновых чисел m = μ и n для значений параметров λ² = 1/2h² = 8, S = 2, R = 1

 

Важно подчеркнуть, что из (20) нельзя однозначно определить параметры наиболее опасного возмущения с максимальной скоростью роста. Видимо, для этого надо рассматривать полную модель с учетом экмановского и внутреннего трения.

Отметим также, что условие возникновения неустойчивости γ₁γ₂ < 0 можно записать в виде $R^2>{\mu }^{-1}\left(S^2{n}^2\right({\lambda }^2-\mu \left)\right)$. С учетом $R=E^{1/2}/2Ro$ это условие дает ограничение сверху на число Россби.

5. АБСОЛЮТНАЯ НЕУСТОЙЧИВОСТЬ БАРОКЛИННОГО ТЕЧЕНИЯ ОТНОСИТЕЛЬНО ДЛИННОВОЛНОВЫХ ВОЗМУЩЕНИЙ

Из уравнений двухуровенной квазигеострофической модели (4) без трения можно получить уравнения, описывающие динамику длинноволновых возмущений. Эти уравнения получаются, если в первом уравнении (4) формально использовать аппроксимацию Δθ << λ²θ. При такой аппроксимации первое уравнение (4) сводится к уравнению ${\theta }_t+\left[\sigma ,\theta \right]=0$, а второе уравнение не меняется. Длинноволновой вариант уравнений двухуровенной модели, таким образом, сводится к системе уравнений

${\theta }_t+\left[\sigma ,\theta \right]=0, {\left(\Delta \sigma \right)}_t+\left[\sigma ,\Delta \sigma \right]+\left[\theta ,\Delta \theta \right]=0$. (21)

Строгий асимптотический вывод уравнений (21) из уравнений непрерывной SQG модели представлен в недавней работе [Kalashnik, 2021]. Полученные в этой работе уравнения, использующие длинноволновые асимптотики решений уравнения Лапласа, отличаются от уравнений (20) лишь численным значением коэффициентов при уравнениях.

Важная особенность системы (21) состоит в том, что из нее следует три квадратичных закона сохранения $d{I}_k/dt=\mathrm{0,} k=\mathrm{1,2,3}$.

$$\begin{gathered} I_1=<{\theta }^2>, I_2=<(\nabla \sigma ,\nabla \theta )>, \\ {I}_3=<{\left|\nabla \sigma \right|}^2-{\left|\nabla \theta \right|}^2>. \end{gathered}$$ (22)

Здесь для краткости угловыми скобками обозначена операция интегрирования по двумерной области $<a>=\iint adxdy$. Первый закон прямо следует из первого уравнения системы (21). Непосредственно из этого закона следует ограниченность во времени решения для θ. Умножение уравнений (21) соответственно на ω = Δσ и θ, сложение и интегрирование дает сохранение $I_2=-<\omega \theta >=<(\nabla \sigma ,\nabla \theta )>$. В справедливости третьего закона можно убедиться непосредственно – дифференцируя выражение для I₃ по времени и преобразуя подынтегральные выражения с использованием уравнений (21). Следуя работе [Kalashnik, 2021], можно показать, что законы сохранения (22) есть асимптотические формы законов сохранения (2), (3) полной двухуровенной модели. Обратим внимание, что в отличие от I₁, законы сохранения I₂, I₃ не являются положительно определенными. Это свидетельствует о внутренней неустойчивой динамике системы.

Как и ранее будем рассматривать систему (1) в цилиндрических координатах r, α, когда x = rcosα, y = rsinα. Аналогично предыдущему, при этом якобиан и оператор Лапласа определяются выражениями $\left[\sigma ,\theta \right]=r^{-1}\left({\sigma }_r{\theta }_{\alpha }-{\sigma }_{\alpha }{\theta }_r\right), \Delta \sigma ={\sigma }_{rr}+r^{-1}{\sigma }_r+r^{-2}{\sigma }_{\alpha \alpha }$. Точное решение системы (20) $\theta =\overline{\theta }\left(r\right)$, $\sigma =\overline{\sigma }\left(r\right)$, отвечает осесимметричному течению с баротропным и бароклинным составляющими азимутальной компоненты скорости $\overline{V}=-{\overline{\sigma }}_r$, $\tilde{V}={\overline{\theta }}_r$. Бароклинному течению отвечает значение $\overline{\sigma }\left(r\right)=0$. Линеаризованная система для возмущений бароклинного течения примет вид

${{\theta }^{\text{'}}}_t-r^{-1}{\overline{\theta }}_r{{\sigma }^{\text{'}}}_{\alpha }=0$,

${\left(\Delta {\sigma }^{\text{'}}\right)}_t+r^{-1}\left({\overline{\theta }}_r{\left(\Delta {\theta }^{\text{'}}\right)}_{\alpha }-{{\theta }^{\text{'}}}_{\alpha }{\left(\Delta \overline{\theta }\right)}_r\right)$. (23)

Дифференцируя второе уравнение (23) по времени, получим уравнение ${\left(\Delta {\sigma }^{\text{'}}\right)}_{tt}+r^{-1}\left({\overline{\theta }}_r{\left(\Delta {{\theta }^{\text{'}}}_t\right)}_{\alpha }-r^{-1}{\overline{\theta }}_r{\left(\Delta \overline{\theta }\right)}_r{{\sigma }^{\text{'}}}_{\alpha \alpha }\right)=0$. Теперь учтем, что из первого уравнения (23) следует ${{\theta }^{\text{'}}}_t=r^{-1}{\overline{\theta }}_r{{\sigma }^{\text{'}}}_{\alpha }$. Подставляя это выражение, окончательно получим одно уравнение для переменной σ'.

${\left(\Delta {\sigma }^{\text{'}}\right)}_{tt}+r^{-1}{\overline{\theta }}_r\left((\Delta {\left(r^{-1}{\overline{\theta }}_r{{\sigma }^{\text{'}}}_{\alpha }\right)}_{\alpha }-{\left(\Delta \overline{\theta }\right)}_r{{\sigma }^{\text{'}}}_{\alpha \alpha }\right)$. (24)

Данное уравнение удобно преобразовать, фиксируя зависимость от угловой координаты, т.е. полагая ${\sigma }^{\text{'}}=\sigma (r,t)\cos \left(n\alpha \right)$. Для амплитуды σ(r, t) получим уравнение

${\left({\Delta }_2\sigma \right)}_{\tau \tau }+n^2{r}^{-1}{\overline{\theta }}_r\left(\sigma {\left({\Delta }_1\overline{\theta }\right)}_r-\left({\Delta }_2\right(r^{-1}{\overline{\theta }}_r\sigma )\right)=0$, (25)

где обозначено ${\Delta }_2\sigma ={\sigma }_{rr}+r^{-1}{\sigma }_r-n^2{r}^{-2}\sigma$, ${\Delta }_1\sigma ={\sigma }_{rr}+r^{-1}{\sigma }_r$.

Рассмотрим в начале частный случай твердотельного вращения $\overline{\theta }\left(r\right)=S{r}^2/2$. Учитывая, что в этом случае $r^{-1}{\overline{\theta }}_r=S$, ${\Delta }_1\overline{\theta }=2S$, из (25) получим уравнение ${\left({\Delta }_2\sigma \right)}_{\tau \tau }-n^2{S}^2{\Delta }_2\sigma =0$, которое после сокращения на оператор Лапласа Δ₂, примет вид ${\sigma }_{\tau \tau }-n^2{S}^2\sigma =0$. Это простейшее уравнение имеет растущее решение σ = Ce^γ,t, где квадрат инкремента нарастания γ² = S²n². Легко видеть, что данное значение квадрата инкремента есть предел общего выражения (13) при μ → 0. Длинноволновой вариант модели, таким образом, приводит к асимптотически точному решению. Отметим, что амплитуда C в линейном приближении не определена.

Используя уравнение (25) можно доказать абсолютную неустойчивость бароклинных течений, т.е. неустойчивость относительно произвольных длинноволновых возмущений и для произвольных профилей скорости. Для доказательства умножим уравнение (25) на rσ и проинтегрируем результат по частям. После простых преобразований получим соотношение

$\begin{array}{c}{\gamma }^2\int (r{\sigma }_r^2+r^{-1}{n}^2{\sigma }^2)dr=\\ =\int (n^{2}r{\left(K\sigma \right)}_r^2+r^{-2}{n}^4{K}^2{\sigma }^2)dr+n^2\int rK{\left({\Delta }_1\overline{\theta }\right)}_r{\sigma }^{2}dr,\end{array}$ (26)

где обозначено $K=r^{-1}{\overline{\theta }}_r$. Обозначим также $A=n^{2}rK{\left({\Delta }_1\overline{\theta }\right)}_r{\sigma }^2=n^2{\overline{\theta }}_r{\left({\Delta }_1\overline{\theta }\right)}_r{\sigma }^2$. Рассмотрим теперь случай течений, слабо отклоняющихся от твердотельного вращения и таких, что ${\overline{\theta }}_r=rF\left(\epsilon r\right)$, где ε – малый параметр. Можно показать, что в этом случае ${\overline{\theta }}_r{\left({\Delta }_1\overline{\theta }\right)}_r=O\left(\epsilon \right)<<1$ и вторым слагаемым в правой части (26) можно пренебречь. При этом из (26) получим соотношение для квадрата инкремента

${\gamma }^2=\frac{{\displaystyle \int (n^{2}r{\left(K\sigma \right)}_r^2+r^{-2}{n}^4{K}^2{\sigma }^2)dr}}{{\displaystyle \int (r{\sigma }_r^2+r^{-1}{n}^2{\sigma }^2)dr}}$. (27)

Отсюда следует что γ² > 0, т.е. всегда существует неустойчивость. В частном случае $F\left(\epsilon r\right)=e^{-\epsilon r}$ мы имеем неустойчивость модифицированного бароклинного течения с твердотельным вращением в центре и затухании на периферии $\tilde{V}={\overline{\theta }}_r=rF\left(\epsilon r\right)=r{e}^{-\epsilon r}$. Напомним, что для возникновения неустойчивости течений однородной жидкости необходимо выполнение специальных условий, связанных с наличием точек перегиба на профиле скорости. Как показано выше, неустойчивость осесимметричных бароклинных течений существует всегда. Этот результат можно доказать и используя общее уравнение (17) для длинноволновых возмущений течения в двухуровенной модели с трением. Для таких возмущений можно считать в уравнении (17) Δθ' = 0. При этом (17) сводится к уравнению

${\left({\theta }^{\text{'}}\right)}_{tt}+S^2{\left({\theta }^{\text{'}}\right)}_{\alpha \alpha }+R{{\theta }^{\text{'}}}_t=0$. (28)

Отыскивая решения (28) вида ${\theta }^{\text{'}}={\theta }^{\text{'}}(r,t)\cos \left(n\alpha \right)$ получим уравнение ${{\theta }^{\text{'}}}_{tt}-n^2{S}^2{\theta }^{\text{'}}+R{{\theta }^{\text{'}}}_t=0$. Если теперь θ' = Ae^γt, то инкремент нарастания  находится из квадратного уравнения ${\gamma }^2+R\gamma -n^2{S}^2=0$, откуда

$\gamma =-(R/2)\pm \sqrt{{(R/2)}^2+n^2{S}^2}$. (29)

Как видно, всегда есть положительный корень γ, отвечающий неустойчивости.

Отметим, что аналогичный результат для плоских бароклинных течений был получен в нашей недавней работе [Kalashnik et. al., 2021].

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

Как известно, в однородной по плотности жидкости баротропное течение с твердотельным вращением устойчиво. В рамках двухуровенной квазигеострофической модели в настоящей работе показано, что это не так в случае бароклинного течения с различным направлением вращения на двух уровнях (течения хетонного типа). За счет разности скоростей по вертикали, такое течение всегда неустойчиво. В работе развита линейная теория неустойчивости таких течений как в модели без трения, так и модели с экмановским трением. Показано, что для неустойчивости в модели с трением горизонтальное волновое число возмущения должно быть меньше некоторого критического значения. Установлено, что неустойчивость относительно длинноволновых возмущений в модели без трения носит абсолютный характер, т.е. существует всегда.

Автор благодарит анонимных рецензентов за полезные замечания.

Работа выполнена при поддержке РНФ (проект № 23-17-00273).

About the authors

M. V. Kalashnik

Obukhov Institute of Atmospheric Physics of the Russian Academy of Sciences; Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences; Typhoon Research and Production Association

Author for correspondence.
Email: kalashnik-obn@mail.ru
Russian Federation, Moscow; Moscow; Obninsk

References

Supplementary files