Стационарные режимы и параметризация экмановского трения в кармановской модели течения вязкой жидкости, возбуждаемого внешней вихревой объемной силой

Полный текст

Введение

Проблема влияния внешнего и внутреннего трения на течения в тонких слоях вязкой несжимаемой жидкости представляет значительный интерес в связи с многочисленными приложениями к задачам гидродинамики и геофизики [Гринспен, 1975]. В качестве примера можно привести пограничные слои в атмосфере, которые являются источником развития движений воздуха, формирующих погоду. Аналогичные по своей постановке задачи возникают и в таких научных направлениях, как астрофизика и магнитогидродинамика [Незлин и др., 1990; Гурбатов и др., 1983].

Для описания крупномасштабных движений в настоящее время общепринятой является концепция линейного трения, восходящая к Чандрасекару [Педлоски, 1984; Вайнштейн и др., 1989; Горькаый и др., 1994]. Это приближение справедливо для течений на мелкой воде при малых числах Рейнольдса. Подробное лабораторное и теоретическое исследование подобного течения вязкой жидкости было проведено Ф.В. Должанским и коллективом сотрудников ИФА им. А.М. Обухова РАН [Kostrykin et al., 2014; Должанский, 2011, 1999; Кострыкин и др., 2011; Пономарев и др., 2009]. При этом было изучено влияние внешнего трения на параметры и устойчивость течений в тонких слоях вращающейся жидкости. В работах Ф.В. Должанского с соавторами на основе анализа вертикальной структуры течения вязкой жидкости в линейном приближении было получено выражение для коэффициента линейного трения и указаны пределы применимости этого приближения в зависимости от чисел Рейнольдса и Россби [Должанский, 1999; Должанский и др., 1990].

При больших числах Рейнольдса и Россби приближение линейного трения является недостаточным. На это было указано в работах [Kostrykin et al., 2014; Кострыкин и др., 2011; Пономарев и др., 2009] в связи с анализом данных лабораторных экспериментов с вихревыми течениями во вращающемся сосуде, который показал различное влияние придонного трения на вихри с циклоническим и антициклоническим вращением.

Течение жидкости в тонком слое над вращающимся дном может рассматриваться как квазидвумерное [Гринспен, 1975]. Такое течение однородной несжимаемой жидкости может быть представлено в виде суммы двумерно-вихревого и двумерно-потенциального течений. При этом скорость в горизонтальной плоскости U имеет вид

$U=k\times \nabla \Psi +\nabla \Phi ,$ (1)

где  $\nabla =i{\partial }_x+j{\partial }_y,\Delta ={\partial }_{xx}+{\partial }_{yy},i,j,k$ – единичные орты декартовой системы координат, $\omega =\Delta \Psi ,\delta =\Delta \Phi$ – вертикальная компонента относительной завихренности и дивергенция горизонтальной скорости соответственно. Вертикальная скорость определяется из уравнения неразрывности: $W_z=-\Delta \Phi$.

Уравнение для горизонтальной скорости однородной вязкой жидкости в слое толщиной $h$, расположенном на вращающемся с частотой $\Omega$ основании, записывается в виде

$\frac{dU}{dt}+2\Omega k\times U+W{U}_z=-\frac{1}{{\rho }_0}\nabla p+\nu U_{zz}+\frac{F_{out}}{{\rho }_0},$ (2)

$U\left(\mathrm{0,}t\right)=\mathrm{0,}{U}_z\left(h,t\right)=\mathrm{0,}W\left(\mathrm{0,}t\right)=W\left(h,t\right)=\mathrm{0,}$ (3)

где $\Omega$ – угловая скорость вращения основания, $p$ – давление, ${\rho }_0$ – постоянная плотность, $\nu$ – коэффициент кинематической вязкости, $F_{out}=k\times \nabla q\left(x,y,t\right)$ – внешняя сила, создающая вихрь с вертикальной осью, $q$ – функция тока внешней силы. Определим следующие параметры, которые мы будем использовать в дальнейшем: толщину экмановского погранслоя – ${\delta }_E={\left(\nu /\Omega \right)}^{1/2}$ и число Экмана (обратное число Рейнольдса) – $E={\mathrm{Re}}^{-1}={\left({\delta }_E/h\right)}^2$.

Внешняя вихревая сила будет создавать первичное течение в виде вихря с вертикальной осью. В то же время граничные условия (3) требуют учета вторичного течения в виде вихря с горизонтальной осью. Связь между двумя вихрями определяется течением в пограничных вязких слоях. Как показывают лабораторные и численные исследования, при этом могут реализовываться различные режимы течений, характеризующиеся разным направлением вращения слоев жидкости [Kostrykin et al., 2014; Кострыкин и др., 2011; Кострыкин, 2018].

Описание этих режимов в рамках чисто двумерной задачи может быть получено, если нам известна связь давления с параметрами поверхностного или осредненного по глубине течения. Такое представление давления даст нам параметризацию исходной задачи, т.е. сокращение первоначально заданных степеней свободы, и позволит определить коэффициент экмановского трения. Если предположить, что числа Экмана и Россби малы, а вертикальный профиль скорости соответствует градиентному экмановскому течению, то после осреднения по вертикали можно получить значение для коэффициента линейного трения [Должанский и др., 1900]

$\lambda =\nu /\left(h{\delta }_E\right)=E^{1/2}\Omega .$ (4)

Выражение для коэффициента экмановского трения при произвольных значениях числа Экмана получено в работе [Козлов и др., 1992]. Кроме того, оставаясь в рамках экмановского профиля скорости в работах [Пермяков и др., 2018; Калашник и др., 2014], получены нелинейные поправки к уравнению квазигеострофического вихря, связанные с трением жидкости о дно. Также в работах [Pedlosky, 2008; Benthuysen et al., 2012] с помощью разложения по малому параметру Россби получены поправки для толщины и скорости накачки экмановского слоя для геострофического течения.

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

Для определения давления через компоненты среднего течения необходимо первоначальное упрощение исходной задачи. Таким упрощением является обращение к модифицированному течению Кармана [Кострыкин и др., 2011], содержащему зависимость течения только от двух основных параметров – чисел Экмана и Россби. Кроме того, в данной модели предполагается осесимметричное течение и линейная зависимость горизонтальной скорости от радиальной координаты, что с достаточно хорошей точностью выполняется в центральной части квазидвумерных вихрей согласно некоторыми теоретическим представлениям, а также данным лабораторных и численных экспериментов [Parfenyev et al., 2012; Орлов и др., 2018].

Анализу различных стационарных режимов течения в рамках этой задачи были посвящены работы [Kostrykin et al., 2014; Кострыкин и др., 2011; Кострыкин, 2018]. Настоящая статья предлагает общий подход к вычислению нелинейного коэффициента трения хотя бы для некоторых режимов вынужденного стационарного течения. Такой подход позволяет объяснить существование разных стационарных режимов течения и допускает обобщение на более широкий класс задач.

Модель Кармановского течения с вихревой объемной вынуждающей силой

В классической постановке модель Кармана используется для описания течения жидкости над вращающимся основанием для слоя бесконечной глубины или течения между двумя вращающимися дисками [Holodniok et al., 1981; Zandbergen, 1980; Чефранов, 2016; Hewitt et al., 2009]. В нашем случае мы используем модифицированную постановку, в которой, во-первых, рассматривается слой жидкости со свободной границей и, во-вторых, течение возбуждается постоянной вихревой силой. Данная постановка оказывается полезной при анализе результатов лабораторных экспериментов по генерации течений магнитогидродинамическим методом, а также их затуханию [Kostrykin et al., 2014; Кострыкин и др., 2011]. Важное преимущество кармановской модели по сравнению с большинством других модельных течений состоит в том, что ее решение являются точным решением 3-мерных уравнений Навье–Стокса и при этом не требуется больших вычислительных затрат для ее численной реализации.

Обратимся к модифицированной задаче Кармана для трехмерного течения, развивающегося под действием вихревой силы в слое глубиной $h$, расположенном над вращающейся с угловой скоростью $\Omega$ бесконечной плоскостью. В цилиндрической системе координат $\left(r\mathrm{,}\phi \mathrm{,}z\right)$ предполагается отсутствие зависимости скорости от азимутального угла, а также линейная (квадратичная) зависимость от радиальной координаты для горизонтальных компонент скорости (давления и функции тока внешней вихревой силы). Уравнения движения записываются для комплекснозначной функции: $K\left(z\right)=G\left(z\right)+iF\left(z\right)$, зависящей только от вертикальной координаты, где $2G$ и $2F$ – вертикальная компонента завихренности и дивергенция горизонтальной скорости, нормированные на характерную угловую скорость вращения ${\Omega }_r$. Также время и вертикальная координата нормируются на характерный временной и пространственный масштабы $1/{\Omega }_r$ и $h$ соответственно.

В нестационарном режиме $K$ удовлетворяет уравнению, следующему из (2)

$K_t=E{K}_{zz}-H{K}_z+i\left(K^2-K_0^2\right),$ (5)

$K_0^2=-P+iQ,$

где $E=1/\mathrm{Re}=\nu /\left({\Omega }_r{h}^2\right)$ – параметр Экмана, $P=-1/\left({\rho }_0{\Omega }_r^2\right)\Delta p$ – константа, характеризующая давление, $Q=1/\left({\rho }_0{\Omega }_r^2\right)\Delta q$ – константа, характеризующая внешнюю вихревую силу, $H=W/\left({\Omega }_{r}h\right)$ – безразмерная вертикальная скорость, удовлетворяющая условию

$H_z=-2F.$ (6)

Граничные условия имеют вид:

$\begin{array}{c}K\left(z=\mathrm{0,}t\right)=S\equiv \frac{\Omega }{{\Omega }_r},H\left(z=\mathrm{0,}t\right)=\mathrm{0,}\\ K_z\left(z=\mathrm{1,}t\right)=\mathrm{0,}H\left(z=\mathrm{1,}t\right)=0.\end{array}$ (7)

Также для нахождения эволюции по времени решения (5)–(7) задаются соответствующие начальные профили $G\left(z\mathrm{,0}\right)$ и $F\left(z\mathrm{,0}\right)$.

Анализируя уравнения (5)–(7), можно заметить, что их решение зависит только от двух из трех исходных параметров задачи $\left(E,Q,S\right)$, в качестве которых можно, например, выбрать пару $\left(Q/E^2,S/E\right)$. Это становится очевидным, если рассмотреть замену переменных

$\begin{array}{c}\tilde{K}=K/E,\tilde{H}=H/E,\tilde{S}=S/E,\tilde{P}=\\ =P/E^2,\tilde{Q}=Q/E^2.\end{array}$

Здесь $\tilde{K}$ явно не зависит от числа Экмана, а только через $\tilde{Q}\mathrm{,}\tilde{S}$.

Анализ численного решения

Для численного решения системы (5)–(7) используется метод, основанный на схеме “кабаре“ для описания адвективных слагаемых и полунеявный метод для описания диффузионных слагаемых [Головизнин и др., 2013]. Решается нестационарная задача с начальным состоянием вида

$G\left(z\mathrm{,0}\right)=S+A\sin \left(\pi z/2\right),F\left(z\mathrm{,0}\right)=0.$

Для достаточно больших моментов времени ($t>50$) численное решение выходит на стационарный режим и предполагается, что такое установившееся решение является стационарным решением системы (5)–(7).

Число узлов сетки по вертикали выбирается достаточно большим ($N=100$), так, чтобы численное решение было устойчивым и не зависело от разрешения сетки. Как показали расчеты [Kostrykin et al., 2014], зависимость стационарного решения задачи от параметра начального распределения $A$ наблюдается только при малых отрицательных значениях параметра $Q$, и при этом эта зависимость существенно не влияет на полученные далее выводы.

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

Как мы видим из рис. 1а, на котором изображено стационарное давление в зависимости от параметра Россби, и в нашем случае стационарное течение имеет два режима. При положительных и больших отрицательных числах Россби – $R{o}_q=signQ{\left|Q\right|}^{1/2}$ давление велико и сильно зависит от этого числа. При малых отрицательных $R{o}_q$ давление близко к нулю. Далее будем называть такие режимы режимами Бэтчелора и Стюартсона соответственно, поскольку похожие стационарные режимы были указаны Бэтчелором [Batchelor, 1951] и Стюартсоном [Stewartson, 1953] для задачи о кармановском течении между двумя дисками.

 

Рис. 1. Результаты численной модели. (a) Зависимость давления от параметра Q при разных числах Экмана. E=1/100 – квадраты, E=1/200 – крестики (S=1, A=0). (б) Диаграмма режимов стационарного течения. Проведены изолинии P/E²=-10 и P/E²=-20 на плоскости параметров (Q/E², S/E)( A=0).

 

Диаграмма режимов стационарного течения на плоскости параметров $\left(Q/E^2,S/E\right)$ представлена на рис. 1б. Фактически там изображены изолинии нормированного давления $P/E^2$. Область малого отрицательного давления, соответствующая режиму Стюартсона, существует при достаточно больших значениях $S/E$ и находится левее оси $Q=0$. Вне этой области наблюдается режим Бэтчелора. Как показано в работе [Кострыкин и др., 2011], при $Q=0$ режим Стюартсона наблюдается в численных расчетах при $S/E\ge 46$.

Из рис. 2, 3 следует, что профили стационарной скорости сильно зависят от числа Экмана. В режиме Стюартсона при уменьшении Е вязкий слой отрывается от границы, что приводит к образованию локального вихря с горизонтальной осью, интенсивность которого может быть сравнима по величине с интенсивностью вихря с вертикальной осью. Очевидно, глубины, где $H=0$, могут являться областями с сильной диссипацией. Расположение этих областей определяет структуру течения, от которой зависит коэффициент трения для поверхностного или среднего по глубине течения.

 

Рис. 2. Профили стационарного численного решения при разных значениях числа Экмана (E=1/100 – сплошная линия, E=1/200 – штриховая, E=1/400 – штрих-пунктирная). Q=1.0, S=1. Вертикальная тонкая пунктирная линия – нулевое значение.

 

Рис. 3. То же, что на рис. 2, только при Q=-1.0, S=1.

 

Режим Бэтчелора

Теоретическое описание стационарного состояния

В режиме Бэтчелора центральное течение имеет простейшую структуру в виде интенсивного вихря с вертикальной осью и слабого вторичного течения. При малых числах Экмана пограничный слой занимает малую часть слоя жидкости и расположен около дна. Решение (5)–(7) представимо в виде

$K\left(z\right)=K_0+k\left(z\right),k\left(z\right)=\left(S-K_0\right)e^{\lambda z},$ (8)

где

$K_0^2=-P+iQ,$ (9)

$-\mathrm{Re}\lambda \gg 1.$

а $k$ удовлетворяет уравнению

$E{k}_{zz}-H{k}_z+i\left(2K_0+k\right)k=0.$ (10)

Его решение ищется при заданном давлении, которое затем определяется из условия $H\left(1\right)=0$.

Из (8), (10) следует, что

${\lambda }_{\mathrm{1,2}}=H/2E\pm {\left[H^2/4E^2-i\left(2K_0+k\right)/E\right]}^{1/2}.$ (11)

Предположим, что ${\lambda }_{\mathrm{1,2}}\left(z\right)\approx {\lambda }_{\mathrm{1,2}}\left(0\right)$. Вблизи нижней границы можно считать $H\approx \mathrm{0,}k\approx S-K_0$, тогда затухающему с высотой стационарному профилю согласно (11) соответствует

$\lambda =-{\left(-i\frac{S+K_0}{E}\right)}^{1/2}.$ (12)

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

$H=-2\left(F_{0}z-\mathrm{Im}\frac{S-K_0}{\lambda }\left(1-e^{\lambda z}\right)\right).$ (13)

Условие $H\left(1\right)=0$ можно записать в виде

$F_0-\left(S-G_0\right)\mathrm{Im}\frac{1}{\lambda }\approx \mathrm{0,}$ (14)

откуда с помощью (12) можно найти связь между $G_0$, $F_0$ и $P$

$\begin{array}{c}{E}^2{\left(S+G_0\right)}^2{\left(S-G_0\right)}^4=\\ =4F_0^3\left(F_0^2+{\left(S+G_0\right)}^2\right)\times \\ \times \left(F_0^3+F_0{\left(S+G_0\right)}^2+E{\left(S-G_0\right)}^2\right).\end{array}$ (15)

Из (9) следуют другие два уравнения, замыкающие систему

$2F_0{G}_0=Q,$

$P=-G_0^2+F_0^2.$ (16)

Решения уравнения (15) представляют собой корни многочлена 8-й степени от переменных $G_0\mathrm{,}{F}_0$ и лежат на двух ветвях (нижней и верхней), проходящих через точку $\left(G_0=S,F_0=0\right)$ (рис. 4). Можно определить, какие ветви решения уравнения (15) реализуются в численной модели, если построить стационарное решение задачи (5)–(7) на плоскости $\left(G_0,F_0\right)$. При построении результатов численного решения предполагается, что $G_0\approx G\left(1\right),F_0\approx F\left(1\right)$, поскольку $\left|k\left(1\right)\right|\ll \left|K_0\right|$ согласно (8), (9). Из рис. 4 мы видим, что в численном решении в режиме Бэтчелора реализуется верхняя ветвь аналитического решения.

 

Рис. 4. Изолинии, соответствующие решению уравнения (15). Значки – точки параметрической кривой (G(1,Q), F(1,Q)), полученной по данным численной модели при -2≤Q≤2 и разных начальных профилях завихренности (A=-2 – кружки и A=0 – плюсы). E=1/50, S=1.

 

Можно оценить величину максимальной вертикальной скорости, генерируемой погранслоем (скорость экмановской подкачки). Подставляя (12), (14) в (13) и пренебрегая первым слагаемым и экспоненциальным множителем во втором, получим

$\begin{array}{c}{H}_{\max }\approx 2\mathrm{Im}\left(\left(S-G_0\right)/\lambda \right)=\\ ={\left(2E\right)}^{1/2}{\left|S+G_0\right|}^{-1/2}\left(G_0-S\right)sign\left(S+G_0\right).\end{array}$ (17)

При малых числах Россби $G_0-S\ll S$

$H_{\max }\approx =E^{1/2}{S}^{-1/2}\left(G_0-S\right).$ (18)

совпадает со скоростью экмановской накачки в линейной теории. При больших числах Россби $\left|G_0\right|\gg S$ выражение (17) примет вид

$H_{\max }\approx {\left(2E\right)}^{1/2}{\left|G_0\right|}^{1/2}sign\left(G_0\right).$ (19)

Из (17) можно сделать вывод, что отношение амплитуды вторичной циркуляции к амплитуде вертикальной завихренности – $|H_{\max }/(G_0-S\left)\right|\le {\left(2E\right)}^{1/2}|S+G_0{|}^{-1/2}$. Поскольку в режиме Бэтчелора $S+G_0$ согласно численным расчетам не может быть малой величиной, то отношение $|H_{\max }/(G_0-S\left)\right|\ll 1$. Следовательно, вторичная циркуляция в этом режиме всегда существенно слабее квазидвухмерного вихря.

Параметризация экмановского трения

Рассмотрим нестационарное решение задачи (5)–(7). Проинтегрировав уравнения (4) по глубине и учитывая граничные условия, получим

${⟨K⟩}_t=-E{K}_z\left(0\right)-⟨2FK⟩+i\left(⟨K^2⟩-K_0^2\right),$ (20)

Здесь $⟨\cdot ⟩$ означает оператор осреднения по глубине.

Предположим, что в режиме Бэтчелора нестационарное решение имеет такой же профиль, как и стационарное решение, но с коэффициентами, зависящими от времени

$K\left(z,t\right)=K_0\left(t\right)+\left(S-K_0\left(t\right)\right)e^{\lambda z},$ (21)

где $\lambda$ определяется по формуле (12).

Учитывая, что $\left|\frac{S-K_0}{K_0\lambda }\right|\ll 1$, а значит $⟨K⟩\approx K_0$, из (20) следует

$K_{0t}=\lambda E\left(K_0-S\right)-2F_0{K}_0.$ (22)

Первое слагаемое в правой части уравнения (22) отвечает за придонное трение, второе – за перенос вторичным течением.

Если рассмотреть решение уравнения (22) вблизи стационарного состояния ${\overline{K}}_0,K_0={\overline{K}}_0+k,\left|k\right|\ll \left|{\overline{K}}_0\right|,k=g+if$, то, пренебрегая квадратичными слагаемыми, получим

$k_t=-\mu k-2{\overline{K}}_{0}f+C,$ (23)

где

$$\begin{gathered} \mu =2{\overline{F}}_0-\overline{\lambda }E\frac{S+3{\overline{K}}_0}{2\left(S+{\overline{K}}_0\right)}, \\ C=-\overline{\lambda }E\left(S-{\overline{K}}_0\right)-2{\overline{F}}_0{\overline{K}}_0, \\ \overline{\lambda }=-{\left(\frac{-i\left(S+{\overline{K}}_0\right)}{E}\right)}^{1/2}. \end{gathered}$$ (24)

Поскольку стационарному состоянию системы отвечает решение $k=0$, то из (23) следует, что $c=0$. Если умножить (23) на $k^{\text{*}}$ и сложить его с комплексно-сопряженным уравнением, то получим

${\left|k\right|}_t^2=-2\mathrm{Re}\mu {\left|k\right|}^2-4f{E}^{-1}\left({\overline{K}}_0{k}^{*}\right).$ (25)

Далее, предполагая, что $k=k_0{e}^{-{\lambda }_{E}t}$, где ${\lambda }_E$ – действительное число, а $k\left(0\right)=k_0$ – начальное значение отклонения решения от стационарного, из (25) следует, что

${\lambda }_E=\mathrm{Re}\mu +\frac{2f_0{E}^{-1}\left({\overline{K}}_0{k}_0^{*}\right)}{|k_0{|}^2}.$ (26)

Заметим, что вторым слагаемым в правой части (26) можно пренебречь, поскольку всегда можно выбрать начальное условие $k_0$ таким образом, чтобы второе слагаемое было существенно меньше первого, например, при $f_0/g_0\ll 1$.

Таким образом, ${\lambda }_E=\mathrm{Re}\mu$ является декрементом затухания завихренности, или, иначе говоря, коэффициентом экмановского трения.

Согласно (12) можно считать, что 

$\mathrm{Re}\overline{\lambda }\approx -{\left(2E\right)}^{-1/2}{\left|S+{\overline{G}}_0\right|}^{1/2}.$

$\begin{array}{c}{\lambda }_E\approx {\left(2E\right)}^{1/2}{\left|S+{\overline{G}}_0\right|}^{-1/2}\times \\ \times sign\left(S+{\overline{G}}_0\right)\left(\frac{7}{4}{\overline{G}}_0-\frac{3}{4}S\right).\end{array}$ (27)

Из формулы (27) следует существование запретной зоны для значений параметра ${\overline{G}}_0$ в режиме Бэтчелора; $-S<{\overline{G}}_0<3/7S$, что подтверждается результатами численных расчетов с этой моделью. В соответствии с построенной теорией, декремент затухания для дивергенции также будет выражаться формулой (27). Далее для удобства обозначения отбросим черту над переменными $G_0\mathrm{,}{F}_0$.

Рассмотрим предельные случаи для формулы (27) при S >0. При $\left|Q\right|\ll \mathrm{1,}\left|G_0-S\right|\ll S,\left|F_0\right|\ll S$.

$F_0=\frac{E^{1/2}}{2S^{1/2}}\left(G_0-S\right),{\lambda }_E\approx {\left(ES\right)}^{1/2}\left(1+\frac{3}{2}\frac{G_0-S}{S}\right)$

или в размерном виде

${\lambda }_E=E^{1/2}\Omega \left(1+\frac{3}{2}Ro\right),$ (28)

где $Ro=G_0-1$ число Россби для стационарного течения.

Для $\left|Q\right|\gg 1$ согласно асимптотическим оценкам

$F_0\approx {(E/2)}^{1/2}{\left|G_0\right|}^{1/2},{\lambda }_E=\frac{7}{2}{\left(E/2\right)}^{1/2}{\left|G_0\right|}^{1/2},$

или в размерном виде

${\lambda }_E\approx \frac{7}{2^{3/2}}{E}^{1/2}\Omega |Ro{|}^{1/2}.$ (29)

Коэффициент экмановского трения (28) при малых числах Россби отличается от коэффициента, полученного с помощью линейной теории (4), только при $Ro\ne 0$. В первом случае за счет добавки в правой части формулы (28) течение будет затухать к стационарному циклоническому состоянию несколько быстрее, чем к антициклоническому. Следует отметить, что максимальное отличие в скорости диссипации для циклонов и антициклонов, вероятней всего, будет наблюдаться при числах Россби порядка –1. Однако в этом случае течение находится в режиме Стюартсона, для которого построенная теория неприменима.

Сравнение декремента затухания завихренности, полученного по формуле (27) с аналогичным коэффициентом, полученным по результатам расчетов численной модели, приведено на рис. 5. Мы видим довольно хорошее согласие для этого параметра при всех значениях $Q$, соответствующих режиму Бэтчелора, как при малых числах Россби, так и при умеренных.

 

Рис. 5. Зависимость декремента затухания по данным численной (сплошная линия) и теоретической (кружки) моделей от параметра Q: (a) -2≤Q≤2, (б) -0.03≤Q≤0.1. E=1/100, S=1.

 

Режим Стюартсона

В режиме Бэтчелора амплитуда вихря с вертикальной осью определяется равновесием между сторонней силой и диссипацией в пограничном слое. Вместе с тем генерируемая в погранслое вертикальная скорость уравновешивается слабым вторичным течением в центральной части. В режиме Стюартсона второе равновесие нарушается и возникает спиральный вихрь, когда вихри с вертикальной и горизонтальной осями уравновешивают друг друга при малом давлении. Похожий режим был указан Стюартсоном на основе линейного решения задачи Кармана с двумя дисками, вращающимися с одинаковой угловой скоростью в противоположных направлениях [Stewartson, 1953]. Только в последнем случае жидкость вне пограничных слоев находится в покое, а в нашем случае, согласно проведенным численным экспериментам, в стационарном режиме присутствует значительная вторичная циркуляция в центральном слое жидкости и в обоих случаях давление мало. Отметим, что данный режим реализуется как для вынужденного течения, так и в случае свободно затухающего течения при сильно антициклонических начальных условиях [Kostrykin et al., 2014].

Рассмотрим стационарное течение в центральном слое. Полное решение уравнения (5) для невязкого течения представимо в виде $K=U+V$, где $U=G^u+i{F}^u$ – свободное течение и $V=G^v+i{F}^v$ – течение, зависящее от величины внешней силы.

Запишем систему уравнений для компоненты течения, не зависящей от форсинга

$-H^u{U}_z+i\left(U^2+P\right)=\mathrm{0,}{H}_z^u=-2\mathrm{Im}U.$

Ее решение имеет вид

$\begin{array}{c}{H}^u=A+B\cos \left(2m\left(z-z_0\right)\right),\\ U=m{H}^u-\frac{i}{2}{H}_z^u=m\left(A+B{e}^{2i\left(z-z_0\right)}\right),\\ P=m^2\left(B^2-A^2\right)\end{array}$ (30)

Оно зависит от четырех произвольных постоянных $A\mathrm{,}B\mathrm{,}m\mathrm{,}{z}_0$. При этом давление в этом режиме близко к нулю, и, кроме того, имеются два условия на вертикальную скорость на границах

$H^u\left(z_{\mathrm{1,2}}\right)\approx 0.$ (31)

Как следует из (30), $P\approx 0$ при $B\approx A$, и, следовательно,

$H^u=2A{\cos }^{2}m\left(z-z_0\right),U=mA\left(1+e^{2i\left(z-z_0\right)}\right).$ (32)

Если считать погранслои достаточно тонкими, то $z_1\approx \mathrm{0,}{z}_2\approx 1$, и условие (31) с учетом (32) запишется в виде

 $m{z}_0=\pi /2+\pi k,m\left(1-z_0\right)=\pi /2+\pi l,$

или 

$m=n\pi ,z_0=\frac{1}{2n}+\frac{k}{n},$ (33)

где $k\mathrm{,}n\ne 0$ — произвольные целые числа, $l=n-k-1$.

Далее по данным численной модели можно построить отношение минимального значения завихренности к максимальному значению вертикальной скорости во всем слое, а также положение экстремумов функций $G\left(z\right)\mathrm{,}H\left(z\right)$ в зависимости от параметра $Q$ (рис. 6).

 

Рис. 6. Зависимость характеристик стационарного решения в режиме Стюартсона от параметров Экмана и величины форсинга Q: (a) координата максимума вертикальной скорости, (б) координата минимума завихренности, (в) максимальное значение H(z), (г) отношение минимального значения G(z) к максимальному значению H(z)·E=1/100 – треугольники, E=1/200 – квадраты, E=1/400 – кружки.

 

Из рис. 6 следует, что $z_0\approx m\approx -\pi$. Этим значениям параметров, согласно (33), соответствуют $k=n=-1$. При этом, как показывают расчеты, значения $z_0\mathrm{,}m$ слабо зависят от чисел Экмана и Россби, более существенным образом от них зависит параметр $H_{max}$. Поскольку $\left|G_{\min }/H_{\max }\right|\approx \pi$, то вертикальная циркуляция составляет существенную долю от горизонтальной в данном режиме, в отличие от стационарного течения в режиме Бэтчелора. Из результатов численных экспериментов следует, что в режиме Стюартсона при достаточно малых числах Экмана и при умеренно отрицательных числах Россби в большой части слоя жидкости между точками $z_1$ и $z_2$ доминирует течение вида (32). В граничных точках $z_i$ невязкое решение сшивается с вязким решением в пограничных слоях и появляется зависимость параметров невязкого течения от параметра Экмана и величины форсинга. Таким образом, получается стационарное решение (5)–(7), на основе которого в линейном приближении можно получить коэффициент экмановского трения аналогично процедуре, приведенной выше для режима Бэтчелора. Однако разработка данной параметризации требует дальнейшего теоретического исследования, и поэтому далее мы приведем только результаты численного моделирования.

На рис. 7 приведен декремент затухания завихренности на поверхности в зависимости от величины форсинга для всего рассматриваемого диапазона параметра $\left|Q\right|\le 2$ при $E=1/100$ и $E=1/200$. Согласно рис. 1а при умеренно отрицательных значениях параметра $Q$ наблюдается режим Стюартсона для стационарного течения. И как следует из рис. 7, именно в данном диапазоне параметра $Q$ декремент затухания завихренности имеет резкий максимум.

 

Рис. 7. Декремент затухания завихренности на поверхности по данным численной модели при разных значениях параметра Q· E=1/100 (звездочки), E=1/200 (квадраты).

 

Некоторые выводы и обобщения

В данной работе показано соответствие стационарных режимов течений задаче о ветровой циркуляции [Кострыкин, 2018] и задаче о течениях, возбуждаемых вихревой объемной силой. Проведено детальное численное исследование стационарных режимов и декремента затухания течения кармановской модели в широком диапазоне значений чисел Россби и при малых числах Экмана. На основе характерного поведения стационарного давления построена диаграмма режимов на плоскости параметров задачи. В отличие от более ранних работ авторов [Kostrykin et al., 2014; Кострыкин и др., 2011] численное исследование стационарных режимов проводилось в более широком диапазоне параметров течения.

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

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

В рамках задачи Кармана аналогичным образом может быть рассмотрена задача об установлении стационарного режима. Для заданного момента времени величина $dK/dt$ рассматривается как входящая в $K_0$ (в $P$ и $Q$). Уравнение стационарного режима при этом становится уравнением эволюции. В случае неоднородного распределения вертикальной завихренности в горизонтальной плоскости, определяя локально по пространству и времени вертикальную структуру, мы получаем систему уравнений для описания поверхностного (среднего) течения, представляющего собой сумму вихревой и потенциальной компонент. Дивергенция потенциальной компоненты скорости зависит от завихренности в фиксированной точке и в свою очередь определяет ее эволюцию. Очевидно, что потенциальная компонента приводит к сжатию (увеличению градиентов) завихренности там, где локально реализуется режим Стюартсона (в области слабых антициклонов). Именно там образуются внутренние (вертикальные) пограничные слои с дополнительной диссипацией. Течение в этих слоях создает нормальную к границе компоненту скорости, которая определяет коэффициент нелинейного трения, учитывающий внешнее и внутреннее трение. С помощью указанного подхода, вероятно, можно улучшить истолкование результатов лабораторных экспериментов с многовихревыми течениями.

Об авторах

С. В. Кострыкин

Институт вычислительной математики им. Г.И. Марчука РАН

Автор, ответственный за переписку.
Email: kostr@mail.ru
Россия, 119333, Москва, ул. Губкина, 8

И. Г. Якушкин

Институт физик и атмосферы им. А.М. Обухова РАН

Email: kostr@mail.ru
Россия, 119017, Москва, Пыжевский пер., 3

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

Дополнительные файлы