Журнал вычислительной математики и математической физики, 2021, T. 61, № 1, стр. 136-149

Спектральный анализ оптимальных возмущений стратифицированного турбулентного течения Куэтта

Г. В. Засько^1, *, Ю. М. Нечепуренко^1, 2, **

¹ Отделение Московского центра фундаментальной и прикладной математики в ИВМ им. Г.И. Марчука РАН
119333 Москва, ул. Губкина, 8, Россия

² ИВМ им. Г.И. Марчука РАН
119333 Москва, ул. Губкина, 8, Россия

^* E-mail: zasko.gr@bk.ru
^** E-mail: yumnech@yandex.ru

Поступила в редакцию 10.12.2019
После доработки 27.07.2020
Принята к публикации 18.09.2020

DOI: 10.31857/S0044466921010105

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

Аннотация

Рассматриваются собственные моды и оптимальные возмущения уравнений стратифицированного турбулентного течения Куэтта, осредненных по горизонтальным пространственным переменным и линеаризованных относительно стационарного состояния. Установлено, что спектр таких уравнений симметричен относительно вещественной оси и лежит строго в левой полуплоскости, т.е. все собственные моды устойчивые, а главная часть оптимального возмущения представляет собой линейную комбинацию большого числа мод, отвечающих собственным значениям с наибольшими вещественными частями. При этом число наиболее значимых мод в этой линейной комбинации растет с ростом числа Рейнольдса. Библ. 20. Фиг. 5. Табл. 2.

Ключевые слова: стратифицированное турбулентное течение Куэтта, мелкомасштабная турбулентность, крупномасштабные структуры, собственные моды, максимальная амплификация, оптимальные возмущения.

1. ВВЕДЕНИЕ

В геофизических пограничных слоях на фоне мелкомасштабной турбулентности часто наблюдаются крупномасштабные структуры. Они вносят существенный вклад в обмены импульсом, теплом и влагой между свободной атмосферой и подстилающей поверхностью (см. [1]). Примером таких структур являются слоистые структуры в поле температуры, наблюдаемые в природе и при прямом (DNS) или вихреразрешающем (LES) численном моделировании устойчиво-стратифицированной атмосферной турбулентности (см. [2]–[4]).

В [3], [5] было проведено DNS-моделирование стратифицированного турбулентного течения Куэтта. Это модельное течение близко к турбулентным течениям в пограничных слоях атмосферы и океана. Установлено, что при больших числах Рейнольдса, в широком диапазоне статической устойчивости, характеризуемой различными значениями числа Ричардсона, наряду с хаотической турбулентностью течение Куэтта содержит крупные структуры, которые могут быть выделены из результатов DNS-моделирования путем разложения мгновенных полей в ряды Фурье по горизонтальным переменным и отбора крупномасштабных гармоник. В случаях, близких к нейтральной стратификации, эти структуры представляют собой крупномасштабные вихри приблизительно круглой формы в поперечном сечении канала, а в случае устойчивой стратификации – крупномасштабные наклонные слои в поле температуры в продольном сечении канала.

Иногда возникновение крупномасштабных структур удается объяснить гидродинамической неустойчивостью осредненного течения (см. [6], [7]). Однако в случае течения Куэтта среднее течение устойчиво. В [8]–[10] образование крупномасштабных структур связывается с возникновением и развитием в мелкомасштабном турбулентном потоке оптимальных возмущений. Оптимальные возмущения вычислялись с помощью технологии, предложенной в [11], [12], на основе уравнений турбулентного течения, осредненных по горизонтальным пространственным переменным и линеаризованных относительно стационарного состояния. Качественное и количественное сравнение крупномасштабных структур с соответствующими им по волновым числам оптимальными возмущениями показало совпадение их пространственных масштабов и конфигураций.

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

2. ОСРЕДНЕННЫЕ УРАВНЕНИЯ ТУРБУЛЕНТНОГО ТЕЧЕНИЯ

Рассмотрим в декартовых координатах $x$ (продольная), $y$ (вертикальная), $z$ (поперечная) турбулентное течение вязкой несжимаемой жидкости в поле силы тяжести в бесконечном трехмерном канале

$\{ - \infty < x < + \infty ,\; - {\kern 1pt} h < y < h,\; - {\kern 1pt} \infty < z < + \infty \} $

полувысоты $h$, верхняя стенка которого движется со скоростью $({{U}_{0}}{\text{/}}2,0,0)$, нижняя – со скоростью $( - {{U}_{0}}{\text{/}}2,0,0)$. На стенках поддерживаются постоянные значения температуры ${{T}_{2}}$ и ${{T}_{1}} < {{T}_{2}}$ соответственно. Определим числа Рейнольдса, Ричардсона и Прандтля как

${\text{Re}} = {{U}_{0}}h{\text{/}}\nu ,\quad {\text{Ri}} = g({{T}_{2}} - {{T}_{1}}){{U}_{0}}{\text{/}}{{T}_{1}}{{h}^{2}},\quad {\text{Pr}} = \nu {\text{/}}\mu ,$

где $g$ – ускорение свободного падения, $\nu $ и $\mu $ – кинематическая вязкость и теплопроводность соответственно.

Следуя [8]–[10], будем считать, что эволюция крупномасштабных составляющих течения в нормированных переменных определяется следующей системой уравнений:

$\begin{gathered} \frac{{\partial U}}{{\partial t}} + \frac{{\partial {{U}^{2}}}}{{\partial x}} + \frac{{\partial UV}}{{\partial y}} + \frac{{\partial UW}}{{\partial z}} + \frac{{\partial P}}{{\partial x}} - \left( {\bar {\nu }\frac{{{{\partial }^{2}}U}}{{\partial {{x}^{2}}}} + \frac{\partial }{{\partial y}}\bar {\nu }\frac{{\partial U}}{{\partial y}} + \bar {\nu }\frac{{{{\partial }^{2}}U}}{{\partial {{z}^{2}}}}} \right) = 0, \\ \frac{{\partial V}}{{\partial t}} + \frac{{\partial UV}}{{\partial x}} + \frac{{\partial {{V}^{2}}}}{{\partial y}} + \frac{{\partial VW}}{{\partial z}} + \frac{{\partial P}}{{\partial y}} - \left( {\bar {\nu }\frac{{{{\partial }^{2}}V}}{{\partial {{x}^{2}}}} + \frac{\partial }{{\partial y}}\bar {\nu }\frac{{\partial V}}{{\partial y}} + \bar {\nu }\frac{{{{\partial }^{2}}V}}{{\partial {{z}^{2}}}}} \right) - \operatorname{Ri} T = 0, \\ \end{gathered} $

(2.1)

$\frac{{\partial W}}{{\partial t}} + \frac{{\partial UW}}{{\partial x}} + \frac{{\partial VW}}{{\partial y}} + \frac{{\partial {{W}^{2}}}}{{\partial z}} + \frac{{\partial P}}{{\partial z}} - \left( {\bar {\nu }\frac{{{{\partial }^{2}}W}}{{\partial {{x}^{2}}}} + \frac{\partial }{{\partial y}}\bar {\nu }\frac{{\partial W}}{{\partial y}} + \bar {\nu }\frac{{{{\partial }^{2}}W}}{{\partial {{z}^{2}}}}} \right) = 0,$

$\begin{gathered} \frac{{\partial T}}{{\partial t}} + \frac{{\partial UT}}{{\partial x}} + \frac{{\partial VT}}{{\partial y}} + \frac{{\partial WT}}{{\partial z}} - \left( {\bar {\mu }\frac{{{{\partial }^{2}}T}}{{\partial {{x}^{2}}}} + \frac{\partial }{{\partial y}}\bar {\mu }\frac{{\partial T}}{{\partial y}} + \bar {\mu }\frac{{{{\partial }^{2}}T}}{{\partial {{z}^{2}}}}} \right) = 0, \\ \frac{{\partial U}}{{\partial x}} + \frac{{\partial V}}{{\partial y}} + \frac{{\partial W}}{{\partial z}} = 0, \\ \end{gathered} $

где взаимодействие с мелкомасштабной турбулентностью параметризовано с помощью коэффициентов турбулентной вязкости $\bar {\nu }$ и теплопроводности $\bar {\mu }$, зависящих только от вертикальной координаты $y$. Здесь $x$, $y$, $z$ – нормированные декартовы координаты; $U$, $V$, $W$, $T$, $P$ – нормированные компоненты вектора скорости, температура и удельное давление соответственно. Верхняя и нижняя стенки канала движутся со скоростями $1{\text{/}}2$ и $ - 1{\text{/}}2$ в направлении $x$, и на них поддерживаются постоянные значения температуры $2$ и $1$ соответственно.

Система уравнений (2.1) имеет стационарное решение вида

(2.2)

$U = \bar {U}(y),\quad V = 0,\quad W = 0,\quad T = \bar {T}(y),\quad P = \bar {P}(y)$

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

(2.3)

$\bar {\nu }\frac{{d\bar {U}}}{{dy}} = {\text{const}},\quad \bar {\mu }\frac{{d\bar {T}}}{{dy}} = {\text{const}},\quad \frac{{d\bar {P}}}{{dy}} = \operatorname{Ri} \bar {T}(y).$

Это течение мы далее будем называть основным.

Представим произвольное решение системы (2.2) в окрестности основного течения в виде

(2.4)

$\begin{gathered} U = \bar {U} + \delta u{\kern 1pt} {\text{'}} + o(\delta ),\quad V = \delta v{\kern 1pt} {\text{'}} + o(\delta ),\quad W = \delta w{\kern 1pt} {\text{'}} + o(\delta ), \\ T = \bar {T} + \delta T{\kern 1pt} {\text{'}} + o(\delta ),\quad P = \bar {P} + \delta p{\kern 1pt} {\text{'}} + o(\delta ), \\ \end{gathered} $

где $\delta $ – малый параметр. Требуя, чтобы для любого сколь угодно малого $\delta $ (2.4) было решением системы (2.2), получаем следующие уравнения распространения малых возмущений:

$\begin{gathered} \frac{{\partial u{\kern 1pt} '}}{{\partial t}} + \bar {U}\frac{{\partial u{\kern 1pt} '}}{{\partial x}} + \frac{{d\bar {U}}}{{dy}}v{\kern 1pt} {\text{'}} + \frac{{\partial p{\kern 1pt} '}}{{\partial x}} - \left( {\bar {\nu }\frac{{{{\partial }^{2}}u{\kern 1pt} '}}{{\partial {{x}^{2}}}} + \frac{\partial }{{\partial y}}\bar {\nu }\frac{{\partial u{\kern 1pt} '}}{{\partial y}} + \bar {\nu }\frac{{{{\partial }^{2}}u{\kern 1pt} '}}{{\partial {{z}^{2}}}}} \right) = 0, \\ \frac{{\partial v{\kern 1pt} '}}{{\partial t}} + \bar {U}\frac{{\partial v{\kern 1pt} '}}{{\partial x}} + \frac{{\partial p{\kern 1pt} '}}{{\partial y}} - \left( {\bar {\nu }\frac{{{{\partial }^{2}}v{\kern 1pt} '}}{{\partial {{x}^{2}}}} + \frac{\partial }{{\partial y}}\bar {\nu }\frac{{\partial v{\kern 1pt} '}}{{\partial y}} + \bar {\nu }\frac{{{{\partial }^{2}}v{\kern 1pt} '}}{{\partial {{z}^{2}}}}} \right) - \operatorname{Ri} T{\kern 1pt} ' = 0, \\ \end{gathered} $

(2.5)

$\frac{{\partial w{\kern 1pt} '}}{{\partial t}} + \bar {U}\frac{{\partial w{\kern 1pt} '}}{{\partial x}} + \frac{{\partial p{\kern 1pt} '}}{{\partial z}} - \left( {\bar {\nu }\frac{{{{\partial }^{2}}w{\kern 1pt} '}}{{\partial {{x}^{2}}}} + \frac{\partial }{{\partial y}}\bar {\nu }\frac{{\partial w{\kern 1pt} '}}{{\partial y}} + \bar {\nu }\frac{{{{\partial }^{2}}w{\kern 1pt} '}}{{\partial {{z}^{2}}}}} \right) = 0,$

$\begin{gathered} \frac{{\partial T{\kern 1pt} '}}{{\partial t}} + \bar {U}\frac{{\partial T{\kern 1pt} '}}{{\partial x}} + \frac{{d\bar {T}}}{{dy}}v{\kern 1pt} {\text{'}} - \left( {\bar {\mu }\frac{{{{\partial }^{2}}T{\kern 1pt} '}}{{\partial {{x}^{2}}}} + \frac{\partial }{{\partial y}}\bar {\mu }\frac{{\partial T{\kern 1pt} '}}{{\partial y}} + \bar {\mu }\frac{{{{\partial }^{2}}T{\kern 1pt} '}}{{\partial {{z}^{2}}}}} \right) = 0, \\ \frac{{\partial u{\kern 1pt} '}}{{\partial x}} + \frac{{\partial v{\kern 1pt} '}}{{\partial y}} + \frac{{\partial w{\kern 1pt} '}}{{\partial z}} = 0. \\ \end{gathered} $

Уравнения (2.5) рассматриваются с нулевыми граничными условиями для $u{\kern 1pt} '$, $v{\kern 1pt} '$, $w{\kern 1pt} '$, $T{\kern 1pt} '$ при $y = \pm 1$.

Нас будут интересовать периодические по $x$ и $z$ возмущения вида

(2.6)

$(u{\kern 1pt} ',v{\kern 1pt} ',w{\kern 1pt} ',T{\kern 1pt} ',p{\kern 1pt} ') = \operatorname{Real} \{ ({{u}_{{\alpha \gamma }}},{{v}_{{\alpha \gamma }}},{{w}_{{\alpha \gamma }}},{{T}_{{\alpha \gamma }}},{{p}_{{\alpha \gamma }}}){{e}^{{{\text{i}}\alpha x + {\text{i}}\gamma z}}}\} ,$

где $\alpha $ и $\gamma $ – вещественные продольное и поперечное волновые числа соответственно, а ${{f}_{{\alpha \gamma }}}$ – комплексные амплитуды, зависящие только от $y$ и $t$. Амплитуды таких возмущений удовлетворяют системе уравнений

$\begin{gathered} \frac{{\partial {{u}_{{\alpha \gamma }}}}}{{\partial t}} + {\text{i}}\alpha \bar {U}{{u}_{{\alpha \gamma }}} + \frac{{d\bar {U}}}{{dy}}{{v}_{{\alpha \gamma }}} + {\text{i}}\alpha {{p}_{{\alpha \gamma }}} + \left( {{{\alpha }^{2}}\bar {\nu }{{u}_{{\alpha \gamma }}} - \frac{\partial }{{\partial y}}\bar {\nu }\frac{{\partial {{u}_{{a\gamma }}}}}{{\partial y}} + {{\gamma }^{2}}\bar {\nu }{{u}_{{\alpha \gamma }}}} \right) = 0, \\ \frac{{\partial {{{v}}_{{\alpha \gamma }}}}}{{\partial t}} + {\text{i}}\alpha \bar {U}{{{v}}_{{\alpha \gamma }}} + \frac{{\partial {{p}_{{\alpha \gamma }}}}}{{\partial y}} + \left( {{{\alpha }^{2}}\bar {\nu }{{{v}}_{{\alpha \gamma }}} - \frac{\partial }{{\partial y}}\bar {\nu }\frac{{\partial {{v}_{{\alpha \gamma }}}}}{{\partial y}} + {{\gamma }^{2}}\bar {\nu }{{{v}}_{{\alpha \gamma }}}} \right) - \operatorname{Ri} {{T}_{{\alpha \gamma }}} = 0, \\ \end{gathered} $

(2.7)

$\frac{{\partial {{w}_{{\alpha \gamma }}}}}{{\partial t}} + {\text{i}}\alpha \bar {U}{{w}_{{\alpha \gamma }}} + {\text{i}}\gamma {{p}_{{\alpha \gamma }}} + \left( {{{\alpha }^{2}}\bar {\nu }{{w}_{{\alpha \gamma }}} - \frac{\partial }{{\partial y}}\bar {\nu }\frac{{\partial {{w}_{{\alpha \gamma }}}}}{{\partial y}} + {{\gamma }^{2}}\bar {\nu }{{w}_{{\alpha \gamma }}}} \right) = 0,$

$\begin{gathered} \frac{{\partial {{T}_{{\alpha \gamma }}}}}{{\partial t}} + {\text{i}}\alpha \bar {U}{{T}_{{\alpha \gamma }}} + \frac{{d\bar {T}}}{{dy}}{{{v}}_{{\alpha \gamma }}} + \left( {{{\alpha }^{2}}\bar {\mu }{{T}_{{\alpha \gamma }}} - \frac{\partial }{{\partial y}}\bar {\mu }\frac{{\partial {{T}_{{\alpha \gamma }}}}}{{\partial y}} + {{\gamma }^{2}}\bar {\mu }{{T}_{{\alpha \gamma }}}} \right) = 0, \\ {\text{i}}\alpha {{u}_{{\alpha \gamma }}} + \frac{{\partial {{{v}}_{{\alpha \gamma }}}}}{{\partial y}} + {\text{i}}\gamma {{w}_{{\alpha \gamma }}} = 0. \\ \end{gathered} $

3. СОБСТВЕННЫЕ МОДЫ И ОПТИМАЛЬНЫЕ ВОЗМУЩЕНИЯ

Решение вида

(3.1)

$({{u}_{{\alpha \gamma }}},{{{v}}_{{\alpha \gamma }}},{{w}_{{\alpha \gamma }}},{{T}_{{\alpha \gamma }}},{{p}_{{\alpha \gamma }}}) = ({{\tilde {u}}_{{\alpha \gamma }}},{{\tilde {v}}_{{\alpha \gamma }}},{{\tilde {w}}_{{\alpha \gamma }}},{{\tilde {T}}_{{\alpha \gamma }}},{{\tilde {p}}_{{\alpha \gamma }}}){{e}^{{\lambda t}}}$

системы (2.7), где $\lambda $ – комплексное число, а ${{\tilde {f}}_{{\alpha \gamma }}}$ – комплексные амплитуды, зависящие только от $y$, называется ее собственной модой, отвечающей собственному значению $\lambda $. Подставляя (3.1) в (2.7), получаем следующую проблему собственных значений:

$\begin{gathered} \lambda {{{\tilde {u}}}_{{\alpha \gamma }}} + {\text{i}}\alpha \bar {U}{{{\tilde {u}}}_{{\alpha \gamma }}} + \frac{{d\bar {U}}}{{dy}}{{{\tilde {v}}}_{{\alpha \gamma }}} + {\text{i}}\alpha {{{\tilde {p}}}_{{\alpha \gamma }}} + \left( {{{\alpha }^{2}}\bar {\nu }{{{\tilde {u}}}_{{\alpha \gamma }}} - \frac{d}{{dy}}\bar {\nu }\frac{{d{{{\tilde {u}}}_{{\alpha \gamma }}}}}{{dy}} + {{\gamma }^{2}}\bar {\nu }{{{\tilde {u}}}_{{\alpha \gamma }}}} \right) = 0, \\ \lambda {{{\tilde {v}}}_{{\alpha \gamma }}} + {\text{i}}\alpha \bar {U}{{{\tilde {v}}}_{{\alpha \gamma }}} + \frac{{d{{{\tilde {p}}}_{{\alpha \gamma }}}}}{{dy}} + \left( {{{\alpha }^{2}}\bar {\nu }{{{\tilde {v}}}_{{\alpha \gamma }}} - \frac{d}{{dy}}\bar {\nu }\frac{{d{{{\tilde {v}}}_{{\alpha \gamma }}}}}{{dy}} + {{\gamma }^{2}}\bar {\nu }{{{\tilde {v}}}_{{\alpha \gamma }}}} \right) - \operatorname{Ri} {{{\tilde {T}}}_{{\alpha \gamma }}} = 0, \\ \end{gathered} $

(3.2)

$\lambda {{\tilde {w}}_{{\alpha \gamma }}} + {\text{i}}\alpha \bar {U}{{\tilde {w}}_{{\alpha \gamma }}} + {\text{i}}\gamma {{\tilde {p}}_{{\alpha \gamma }}} + \left( {{{\alpha }^{2}}\bar {\nu }{{{\tilde {w}}}_{{\alpha \gamma }}} - \frac{d}{{dy}}\bar {\nu }\frac{{d{{{\tilde {w}}}_{{\alpha \gamma }}}}}{{dy}} + {{\gamma }^{2}}\bar {\nu }{{{\tilde {w}}}_{{\alpha \gamma }}}} \right) = 0,$

$\begin{gathered} \lambda {{{\tilde {T}}}_{{\alpha \gamma }}} + {\text{i}}\alpha \bar {U}{{{\tilde {T}}}_{{\alpha \gamma }}} + \frac{{d\bar {T}}}{{dy}}{{{\tilde {v}}}_{{\alpha \gamma }}} + \left( {{{\alpha }^{2}}\bar {\mu }{{{\tilde {T}}}_{{\alpha \gamma }}} - \frac{d}{{dy}}\bar {\mu }\frac{{d{{{\tilde {T}}}_{{\alpha \gamma }}}}}{{dy}} + {{\gamma }^{2}}\bar {\mu }{{{\tilde {T}}}_{{\alpha \gamma }}}} \right) = 0, \\ {\text{i}}\alpha {{{\tilde {u}}}_{{\alpha \gamma }}} + \frac{{d{{{\tilde {v}}}_{{\alpha \gamma }}}}}{{dy}} + {\text{i}}\gamma {{{\tilde {w}}}_{{\alpha \gamma }}} = 0. \\ \end{gathered} $

Теорема 1. Пусть $\bar {U}(y)$, $\bar {T}(y)$, $\bar {\nu }(y)$, $\bar {\mu }(y)$ удовлетворяют соотношениям (2.3), причем профили скорости и температуры являются нечетными функциями. Тогда спектр проблемы собственных значений (3.2) симметричен относительно вещественной оси, и если

$(\lambda ,{{\tilde {u}}_{{\alpha \gamma }}}(y),{{\tilde {v}}_{{\alpha \gamma }}}(y),{{\tilde {w}}_{{\alpha \gamma }}}(y),{{\tilde {T}}_{{\alpha \gamma }}}(y),{{\tilde {p}}_{{\alpha \gamma }}}(y))$

является решением проблемы (3.2), то

$(\lambda {\kern 1pt} *,{{\tilde {u}}_{{\alpha \gamma }}}( - y){\kern 1pt} *,{{\tilde {v}}_{{\alpha \gamma }}}( - y){\kern 1pt} *,{{\tilde {w}}_{{\alpha \gamma }}}( - y){\kern 1pt} *,{{\tilde {T}}_{{\alpha \gamma }}}( - y){\kern 1pt} *, - {{\tilde {p}}_{{\alpha \gamma }}}( - y){\kern 1pt} *)$

также является решением проблемы (3.2), где $ * $ означает комплексное сопряжение.

Отметим, что приведенная теорема доказывается комплексным сопряжением уравнений (3.2) и заменой переменных $y \to - y$. Она является непосредственным обобщением соответствующего утверждения, справедливого для классического течения Куэтта (см. [13]), т.е. ламинарного течения с ${\text{Ri}} = 0$ и $\bar {\nu } = 1{\text{/Re}} = {\text{const}}$. Классическое течение Куэтта имеет линейный профиль продольной компоненты скорости $\bar {U}(y) = y{\text{/}}2$. Известно, что оно устойчиво при любом числе Рейнольдса, т.е. спектр, соответствующий проблеме (3.2), лежит строго в левой полуплоскости [14]. Доказательством последнего утверждения для проблемы (3.2) в условиях теоремы 1 мы не располагаем, но, как показывают приведенные ниже результаты расчетов, оно по-видимому верно.

Обозначим через $r_{{max}}^{{\alpha \gamma }}$ максимальную вещественную часть собственных значений проблемы (3.2) при фиксированных значениях волновых чисел $\alpha $, $\gamma $, а через

${{r}_{{max}}} = \mathop {max}\limits_{\alpha ,\gamma } r_{{max}}^{{\alpha \gamma }}$

глобальную максимальную вещественную часть соответственно, где максимум берется по всем вещественным $\alpha $ и $\gamma $. Собственную моду вида (3.1), на которой достигается ${{r}_{{max}}}$, будем далее называть глобальной ведущей.

Определим среднюю плотность полной энергии возмущения вида (2.6) в момент времени $t$ как

(3.3)

${{\mathcal{E}}_{t}} = \frac{1}{2}\int\limits_{ - 1}^1 \,\left( {{{{\left| {{{u}_{{\alpha \gamma }}}} \right|}}^{2}} + {{{\left| {{{{v}}_{{\alpha \gamma }}}} \right|}}^{2}} + {{{\left| {{{w}_{{\alpha \gamma }}}} \right|}}^{2}} + \frac{{{\text{Ri}}}}{{d\bar {T}{\text{/}}dy}}{{{\left| {{{T}_{{\alpha \gamma }}}} \right|}}^{2}}} \right)dy.$

Максимально возможное увеличение

${{\Gamma }^{{\alpha \gamma }}}(t) = max\frac{{{{\mathcal{E}}_{t}}}}{{{{\mathcal{E}}_{0}}}}$

средней плотности полной энергии возмущения, где максимум берется по всем решениям системы (2.7), будем называть максимальной амплификацией средней плотности полной энергии при фиксированных значениях $\alpha $, $\gamma $ и $t$. Введем обозначения

$\Gamma _{{max}}^{{\alpha \gamma }} = \mathop {max}\limits_{t \geqslant 0} {{\Gamma }^{{\alpha \gamma }}}(t),\quad {{\Gamma }_{{max}}} = \mathop {max}\limits_{\alpha ,\gamma } \Gamma _{{max}}^{{\alpha \gamma }},\quad ({{\alpha }_{{{\text{opt}}}}},{{\gamma }_{{{\text{opt}}}}},{{t}_{{{\text{opt}}}}}) = \mathop {\arg \max }\limits_{\alpha ,\gamma ,t} {{\Gamma }^{{\alpha \gamma }}}(t)$

для максимальной амплификации средней плотности полной энергии при фиксированных волновых числах, глобальной максимальной амплификации и оптимальных значений волновых чисел и оптимального момента времени соответственно. Начальное возмущение, на котором достигается ${{\Gamma }_{{max}}}$, будем называть глобальным оптимальным возмущением.

Для классического течения Куэтта известно (см. [15]), что глобальная максимальная амплификация достигается при значении продольного волнового числа $\alpha $, близкого к нулю, а глобальное оптимальное возмущение представляет собой крупномасштабные чередующиеся по направлению вращения вихри в поперечном сечении канала.

4. ЧИСЛЕННАЯ МОДЕЛЬ

Следуя работам [8]–[10], пространственную аппроксимацию по $y$ системы уравнений (2.7) будем выполнять методом Галеркина-коллокаций, описанным в [16], [17]. В качестве узлов сетки для давления выберем корни ${{y}_{1}} < \cdots < {{y}_{n}}$ производной $L_{{n + 1}}^{'}$ многочлена Лежандра степени $n + 1$, а в качестве узлов сетки для компонент скорости и температуры – те же узлы вместе с точками ${{y}_{0}} = - 1$ и ${{y}_{{n + 1}}} = 1$ (узлы Гаусса–Лобатто). В качестве базисных функций для компонент скорости и температуры с учетом нулевых граничных условий будем использовать элементарные интерполяционные многочлены Лагранжа, представимые в виде

${{\psi }_{i}}(y) = \frac{{({{y}^{2}} - 1)L_{{n + 1}}^{'}(y)}}{{(n + 1)(n + 2)(y - {{y}_{i}}){{L}_{{n + 1}}}({{y}_{i}})}},$

а для давления – элементарные интерполяционные многочлены Лагранжа, представимые в виде

${{\varphi }_{i}}(y) = \frac{{(y_{i}^{2} - 1)L_{{n + 1}}^{'}(y)}}{{(n + 1)(n + 2)(y - {{y}_{i}}){{L}_{{n + 1}}}({{y}_{i}})}}.$

Таким образом, компоненты амплитуд возмущения будем аппроксимировать как

$g(y,t) \approx \sum\limits_{i = 1}^n \,{{g}_{i}}(t){{\psi }_{i}}(y),\quad {{p}_{{\alpha \gamma }}}(y,t) \approx \sum\limits_{i = 1}^n \,{{p}_{{\alpha \gamma ,i}}}(t){{\varphi }_{i}}(y),$

где $g$ означает ${{u}_{{\alpha \gamma }}}$, ${{{v}}_{{\alpha \gamma }}}$, ${{w}_{{\alpha \gamma }}}$ или ${{T}_{{\alpha \gamma }}}$. Коэффициенты ${{g}_{i}}(t)$, ${{p}_{i}}(t)$ при такой аппроксимации являются значениями аппроксимантов в узле ${{y}_{i}}$.

В качестве пробных функций для каждого из первых четырех уравнений в (2.7) будем использовать функции ${{\psi }_{i}}(y)$, а для уравнения неразрывности – ${{\varphi }_{i}}(y)$. Для расчета фигурирующих в слабой постановке скалярных произведений будем использовать квадратурную формулу с узлами и весами Гаусса–Лобатто (см. [16]):

$\int\limits_{ - 1}^1 \,f(y)dy \approx \sum\limits_{k = 0}^{n + 1} \,{{\kappa }_{i}}f({{y}_{i}}),\quad {{\kappa }_{i}} = \frac{2}{{(n + 1)(n + 2)L_{{n + 1}}^{2}({{y}_{i}})}},$

точную для многочленов от $y$ степени не выше $2n + 1$.

Обозначим через ${{K}_{0}}$ положительно определенную диагональную матрицу порядка $n + 2$ квадратурных коэффициентов, а через $K$ – ее подматрицу порядка $n$, отвечающую внутренним узлам. Введем также следующие диагональные матрицы порядка $n$: $U$, ${{U}_{y}}$, $N$, $M$ и ${{T}_{y}}$, составленные соответственно из значений профиля $\bar {U}$ и производной $d\bar {U}{\text{/}}dy$ профиля продольной компоненты скорости основного течения, значений коэффициентов турбулентной вязкости $\bar {\nu }$ и диффузии $\bar {\mu }$ и значений производной $d\bar {T}{\text{/}}dy$ профиля температуры основного течения в узлах ${{y}_{1}} < \ldots < {{y}_{n}}$, а также диагональные матрицы ${{N}_{0}}$, ${{M}_{0}}$ порядка $n + 2$, составленные из значений коэффициентов турбулентной вязкости и диффузии в узлах сетки ${{y}_{0}} < \ldots < {{y}_{{n + 1}}}$. Для вычисления значений в узлах ${{y}_{0}} < \ldots < {{y}_{{n + 1}}}$ производной функции, заданной в узлах ${{y}_{1}} < \ldots < {{y}_{n}}$ и удовлетворяющей нулевым граничным условиям, будем использовать матрицу дифференцирования $D$ размера $(n + 2) \times n$. Так же нам потребуется матрица проектирования $P$ размера $(n + 2) \times n$, восстанавливающая по значениям функции в узлах ${{y}_{1}} < \ldots < {{y}_{n}}$ ее значения в узлах ${{y}_{0}} < \ldots < {{y}_{{n + 1}}}$. Эффективные методы вычисления матриц $D$ и $P$ на основе интерполяционных многочленов Лагранжа описаны в [18].

Выполнив дискретизацию системы уравнений (2.7) методом галеркина-коллокаций, получим систему обыкновенных дифференциальных и алгебраических уравнений, которую можно привести к следующему виду:

(4.1)

$\frac{{d{\mathbf{v}}}}{{dt}} = J{\mathbf{v}} - Gp,\quad F{\mathbf{v}} = 0,$

где

${\mathbf{v}} = {{E}^{{1/2}}}{{({{u}^{{\text{T}}}},{{v}^{{\text{T}}}},{{w}^{{\text{T}}}},{{T}^{{\text{T}}}})}^{{\text{T}}}},\quad E = \frac{1}{2}{\text{diag}}(K,K,K,\operatorname{Ri} KT_{y}^{{ - 1}}),$

а $u$, $v$, $w$, $T$, $p$ суть $n$-компонентные векторные функции, компонентами которых являются соответственно значения амплитуд ${{u}_{{\alpha \gamma }}}$, ${{v}_{{\alpha \gamma }}}$, ${{w}_{{\alpha \gamma }}}$, ${{T}_{{\alpha \gamma }}}$, ${{p}_{{\alpha \gamma }}}$ в узлах ${{y}_{1}} < \ldots < {{y}_{n}}$. Отметим, что при выбранной нормировке дискретным аналогом функционала (3.3) средней плотности полной энергии будет $\left\| {\mathbf{v}} \right\|_{2}^{2}$.

Матрицы в (4.1) устроены следующим образом:

$J = \left[ {\begin{array}{*{20}{c}} {{{S}_{\nu }}}&{ - {{U}_{y}}}&0&0 \\ 0&{{{S}_{\nu }}}&0&R \\ 0&0&{{{S}_{\nu }}}&0 \\ 0&{ - R}&0&{{{S}_{\mu }}} \end{array}} \right],\quad G = \left[ {\begin{array}{*{20}{c}} {{\text{i}}\alpha I} \\ G \\ {{\text{i}}\gamma I} \\ 0 \end{array}} \right],\quad F = - G{\kern 1pt} *,$

где $I$ – единичная матрица порядка $n$, $G = - {{K}^{{ - 1/2}}}{{D}^{{\text{T}}}}{{K}_{0}}P{{K}^{{ - 1/2}}}$,

${{S}_{\nu }} = - {\text{i}}\alpha U - {{\alpha }^{2}}N + {{L}_{\nu }} - {{\gamma }^{2}}N,\quad {{S}_{\mu }} = T_{y}^{{ - 1/2}}( - {\text{i}}\alpha U - {{\alpha }^{2}}M + {{L}_{\mu }} - {{\gamma }^{2}}M)T_{y}^{{1/2}},$

${{L}_{\nu }} = - {{K}^{{ - 1/2}}}{{D}^{{\text{T}}}}{{K}_{0}}{{N}_{0}}D{{K}^{{ - 1/2}}},\quad {{L}_{\mu }} = - {{K}^{{ - 1/2}}}{{D}^{{\text{T}}}}{{K}_{0}}{{M}_{0}}D{{K}^{{ - 1/2}}}.$

Отметим, что матрицы $G$ и $F$ являются дискретными аналогами оператора $\partial {\text{/}}\partial y$ в градиенте давления и уравнении неразрывности соответственно, ${{L}_{\nu }}$ и ${{L}_{\mu }}$ – соответственно дискретные аналоги операторов

$\frac{\partial }{{\partial y}}\bar {\nu }\frac{\partial }{{\partial y}},\quad \frac{\partial }{{\partial y}}\bar {\mu }\frac{\partial }{{\partial y}},$

а $R = \sqrt {\operatorname{Ri} } T_{y}^{{1/2}}$.

Из второго уравнения системы (4.1) следует, что ее решение принадлежит ядру матрицы $F$. После замены переменных ${\mathbf{v}} = V{\mathbf{u}}$, где $V$ – прямоугольная матрица, столбцы которой образуют ортонормированный базис в ядре матрицы $F$, и умножения полученного уравнения слева на $V{\kern 1pt} *$, а также с учетом того, что $G = - F{\kern 1pt} *$, мы получим следующую систему обыкновенных дифференциальных уравнений:

(4.2)

$\frac{{d{\mathbf{u}}}}{{dt}} = H{\mathbf{u}},$

где $H = V{\kern 1pt} *JV$ – квадратная комплексная матрица порядка $3n + 1$ при $\alpha = \gamma = 0$ и порядка $3n$ в остальных случаях. Подробное обоснование такого типа редукций линейных дифференциально-алгебраических систем дано в [19].

Для вычисления вектора значений амплитуды $({{u}_{{\alpha \gamma }}},{{v}_{{\alpha \gamma }}},{{w}_{{\alpha \gamma }}},{{T}_{{\alpha \gamma }}})$ возмущения вида (2.6) в узлах расчетной сетки необходимо сделать обратную замену переменных

(4.3)

${{({{u}^{{\text{T}}}},{{v}^{{\text{T}}}},{{w}^{{\text{T}}}},{{T}^{{\text{T}}}})}^{{\text{T}}}} = {{E}^{{ - 1/2}}}V{\mathbf{u}}.$

5. ВЫЧИСЛЕНИЕ СОБСТВЕННЫХ МОД И ОПТИМАЛЬНЫХ ВОЗМУЩЕНИЙ

Каждой собственной паре $(\lambda ,{\mathbf{\tilde {u}}})$ матрицы $H$ системы (4.2) соответствует (с точностью до погрешности аппроксимации) собственная мода (3.1) исходной системы (2.7). Ее амплитуда в узлах расчетной сетки может быть вычислена по формуле (4.3) с ${\mathbf{u}} = {\mathbf{\tilde {u}}}$.

В силу унитарной инвариантности второй нормы $\left\| {{\mathbf{u}}(t)} \right\|_{2}^{2} = \left\| {{\mathbf{v}}(t)} \right\|_{2}^{2}$. Следовательно, $\left\| {{\mathbf{u}}(t)} \right\|_{2}^{2}$ является дискретным аналогом средней плотности полной энергии соответствующего возмущения вида (2.6). Поскольку произвольное решение системы (4.2) представимо в виде

${\mathbf{u}}(t) = exp\{ tH\} {{{\mathbf{u}}}^{0}},$

с точностью до погрешности аппроксимации

${{\Gamma }^{{\alpha \gamma }}}(t) = \left\| {\exp \{ tH\} } \right\|_{2}^{2},\quad \Gamma _{{max}}^{{\alpha \gamma }} = \mathop {max}\limits_{t \geqslant 0} {{\Gamma }^{{\alpha \gamma }}}(t).$

Таким образом, вычисление максимальной амплификации средней плотности полной энергии возмущений сводится к вычислению для заданной квадратной комплексной матрицы $H$ максимума функции ${{\Gamma }^{{\alpha \gamma }}}(t)$ при $t \geqslant 0$. Для нахождения $t = t_{{{\text{opt}}}}^{{\alpha \gamma }}$, дающего максимум ${{\Gamma }^{{\alpha \gamma }}}(t)$ с заданной относительной точностью, будем использовать эффективный алгоритм, предложенный в [12] и основанный на малоранговой аппроксимации. После того как $t_{{{\text{opt}}}}^{{\alpha \gamma }}$ найдено, вычисляем максимальное сингулярное число ${{\sigma }_{{{\text{opt}}}}}$ и отвечающий ему правый сингулярный вектор ${{{\mathbf{u}}}_{{{\text{opt}}}}}$ матрицы $\exp \{ t_{{{\text{opt}}}}^{{\alpha \gamma }}H\} .$ Максимальная амплификация $\Gamma _{{max}}^{{\alpha \gamma }}$ равна $\sigma _{{{\text{opt}}}}^{2}$, а амплитуда оптимального возмущения в узлах расчетной сетки может быть вычислена по формуле (4.3) с ${\mathbf{u}} = {{{\mathbf{u}}}_{{{\text{opt}}}}}$.

Пусть $\Lambda $ означает некоторое изолированное подмножество спектра матрицы $H$, $\mathcal{U}$ – инвариантное подпространство, отвечающее $\Lambda $, а ${{P}_{\mathcal{U}}}$ – соответствующий спектральный проектор, т.е. проектор на $\mathcal{U}$, коммутирующий с матрицей $H$. Нас будет интересовать квадрат нормы

${{c}_{\mathcal{U}}}(t) = \left\| {{{P}_{\mathcal{U}}}exp\{ tH\} {{{\mathbf{u}}}_{{{\text{opt}}}}}} \right\|_{2}^{2}$

проекции глобального оптимального возмущения на подпространство $\mathcal{U}$ в моменты времени $t = 0$ и $t_{{{\text{opt}}}}^{{\alpha \gamma }}$ и максимальная амплификация

${{\Gamma }_{{\mathcal{U},max}}} = \mathop {max}\limits_{t \geqslant 0} \mathop {max}\limits_{{\mathbf{w}} \in \mathcal{U},{\text{||}}{\mathbf{w}}{\text{|}}{{{\text{|}}}_{2}} = 1} \left\| {exp\{ tH\} {\mathbf{w}}} \right\|_{2}^{2}$

на этом подпространстве. Эти величины мы будем вычислять на основе разложения Шура (см. [20])

(5.1)

$H = [{{Q}_{1}},{{Q}_{2}}]\left[ {\begin{array}{*{20}{c}} {{{S}_{1}}}&{{{S}_{{12}}}} \\ 0&{{{S}_{2}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {Q_{1}^{*}} \\ {Q_{2}^{*}} \end{array}} \right],$

где ${{S}_{1}}$ и ${{S}_{2}}$ – верхние треугольные матрицы, причем $\lambda ({{S}_{1}}) = \Lambda $, а $[{{Q}_{1}},{{Q}_{2}}]$ – унитарная матрица, разбитая на два блока в соответствии с разбиением на блоки формы Шура. Максимальная амплификация на подпространстве $\mathcal{U}$ равна

${{\Gamma }_{{\mathcal{U},max}}} = \mathop {max}\limits_{t \geqslant 0} \left\| {exp\{ t{{S}_{1}}\} } \right\|_{2}^{2},$

${{c}_{\mathcal{U}}}(t) = \left\| {{{Q}_{1}}exp\{ t{{S}_{1}}\} (Q_{1}^{*} - MQ_{2}^{*}){{{\mathbf{u}}}_{{{\text{opt}}}}}} \right\|_{2}^{2},$

где $M$ – решение уравнения Сильвестра

(5.2)

$M{{S}_{2}} - {{S}_{1}}M = {{S}_{{12}}}.$

Отметим, что уравнение (5.2) однозначно разрешимо, так как спектры матриц ${{S}_{1}}$ и ${{S}_{2}}$ не пересекаются.

6. РЕЗУЛЬТАТЫ РАСЧЕТОВ

Профили продольной компоненты скорости $\bar {U}(y)$ и температуры $\bar {T}(y)$ основного течения, а также соответствующие коэффициенты турбулентной вязкости $\bar {\nu }(y)$ и теплопроводности $\bar {\mu }(y)$ возьмем из результатов работы [5], где было выполнено прямое численное моделирование стратифицированного турбулентного течения Куэтта в широких диапазонах чисел Рейнольдса $1 \leqslant {\text{Re}} \times {{10}^{{ - 4}}} \leqslant 6$ и Ричардсона $0 \leqslant {\text{Ri}} \leqslant 0.12$. В данной работе мы ограничимся рассмотрением устойчиво-стратифицированного течения при фиксированном числе Ричардсона ${\text{Ri}} = 0.03$ и значениях числа Рейнольдса ${\text{Re}} = 2 \times {{10}^{4}}$, $4 \times {{10}^{4}}$ и $6 \times {{10}^{4}}$. Профили $\bar {U}(y)$ и $\bar {\nu }(y)$ при указанных значениях $\operatorname{Ri} $ и $\operatorname{Re} $ изображены на фиг. 1. Профили $\bar {T}(y)$ и $\bar {\mu }(y)$ выглядят аналогично.

Фиг. 1.

Профили $\bar {U}(y)$ (а) и $\bar {\nu }(y)$ (б) основного течения при ${\text{Ri}} = 0.03$ и ${\text{Re}} = 2 \times {{10}^{4}}$ (красным), $4 \times {{10}^{4}}$ (синим) и $6 \times {{10}^{4}}$ (зеленым).

6.1. Сравнение ведущих мод и глобальных оптимальных возмущений

Линии уровня $\Gamma _{{max}}^{{\alpha \gamma }}$ при рассматриваемых значениях чисел Ричардсона и Рейнольдса изображены на фиг. 2. Значения глобальной максимальной амплификации ${{\Gamma }_{{max}}}$, оптимальных волновых чисел и оптимального времени $({{\alpha }_{{{\text{opt}}}}},{{\gamma }_{{{\text{opt}}}}},{{t}_{{{\text{opt}}}}})$ приведены в табл. 1. Видно, что с увеличением числа Рейнольдса наблюдается увеличение ${{\Gamma }_{{max}}}$. Глобальная максимальная амплификация достигется при ненулевом поперечном волновом числе $\gamma $ при всех рассмотренных числах Рейнольдса и ненулевом продольном волновом числе $\alpha $ при всех рассмотренных числах Рейнольдса, кроме минимального. Причем,

Фиг. 2.

Линии уровня $\Gamma _{{max}}^{{\alpha \gamma }}$ в плоскости $(\alpha ,\gamma )$ при ${\text{Ri}} = 0.03$ и ${\text{Re}} = 2 \times {{10}^{4}}$ (а), $4 \times {{10}^{4}}$ (б) и $6 \times {{10}^{4}}$ (в).

Таблица 1.

Значения ${{\Gamma }_{{max}}}$, $({{\alpha }_{{{\text{opt}}}}},{{\gamma }_{{{\text{opt}}}}},{{t}_{{{\text{opt}}}}})$ и ${{r}_{{max}}}$ в зависимости от числа Рейнольдса ${\text{Re}}$ при ${\text{Ri}} = 0.03$

Re	${{\Gamma }_{{max}}}$	$({{\alpha }_{{{\text{opt}}}}},{{\gamma }_{{{\text{opt}}}}},{{t}_{{{\text{opt}}}}})$	${{r}_{{max}}}$
2 × 10⁴	24.8383	(0.0000, 1.1165, 57.6)	–0.000847
4 × 10⁴	30.4001	(0.3890, 1.1606, 46.6)	–0.000683
6 × 10⁴	38.0889	(0.4664, 1.1424, 46.8)	–0.000553

с увеличением числа Рейнольдса наблюдается увеличение оптимального продольного волнового числа ${{\alpha }_{{{\text{opt}}}}}$, а оптимальное поперечное волновое число ${{\gamma }_{{{\text{opt}}}}}$ почти не меняется.

Отметим, что согласно [9], [10], оптимальные возмущения при $\alpha = 0$ представляют собой ролики – крупномасштабные вихри приблизительно круглой формы в поперечном сечении канала, а оптимальные возмущения при $\alpha > 0$ – крупномасштабные наклонные слоистые структуры. Причем ролики развиваются во времени заметно медленнее, чем слоистые структуры, что соответствует приведенным в табл. 1 оптимальным временам. Таким образом, с увеличением числа Рейнольдса наблюдается существенное изменение структуры глобального оптимального возмущения.

В табл. 1 также представлены результаты вычисления глобальной максимальной вещественной части ${{r}_{{max}}}$ собственных значений матрицы $H$, т.е. инкременты нарастания глобальных ведущих мод. Видно, что как и глобальная максимальная амплификация, величина ${{r}_{{max}}}$ возрастает с ростом числа Рейнольдса. Однако она остается отрицательной. Отметим, что величина ${{r}_{{max}}}$ достигалась при значениях волновых чисел $\alpha = \gamma = 0$ при всех рассмотренных числах Рейнольдса. Таким образом, глобальные ведущие моды существенно отличаются от глобальных оптимальных возмущений даже по волновым числам, т.е. глобальные ведущие моды не входят в состав глобальных оптимальных возмущений. Для сравнения ведущих мод и оптимальных возмущений при оптимальных значениях волновых чисел обратимся к фиг. 3 и 5 .

Фиг. 3.

Ведущая часть спектра матрицы $H$ при ${\text{Ri}} = 0.03$, ${\text{Re}} = 2 \times {{10}^{4}}$ (а), $4 \times {{10}^{4}}$ (б), $6 \times {{10}^{4}}$ (в) и оптимальных значениях волновых чисел $({{\alpha }_{{{\text{opt}}}}},{{\gamma }_{{{\text{opt}}}}})$, вычисленная на сетках с $n = 100$ (“+”) и $200$ (“o”). Подмножества ${{\Lambda }_{{11}}}$ и ${{\Lambda }_{{12}}}$ выделены красным и синим цветами соответственно.

Фиг. 4.

Результаты при ${\text{Re}} = 4 \times {{10}^{4}}$, ${\text{Ri}} = 0.03$ и оптимальных значениях волновых чисел. Слева: абсолютные величины амплитуд компонент нормированной ведущей вещественной моды (синим) и нормированных ведущих комплексно-сопряженных мод (зеленым и красным – моды, отвечающие собственным значениям с положительными и отрицательными мнимыми частями соответственно); справа: абсолютные величины амплитуд компонент нормированного глобального оптимального возмущения в момент времени $t = 0$ (штрихпунктир) и нормированного глобального оптимального возмущения в моменты времени ${{t}_{{{\text{opt}}}}}$ (сплошная линия).

Фиг. 5.

Абсолютные величины амплитуд компонент нормированного глобального оптимального возмущения (черным) и оптимального на подпространстве ${{\Lambda }_{1}}$ (красным) в моменты времени $t = 0$ (штриховая) и ${{t}_{{{\text{opt}}}}}$ (сплошная линия) при ${\text{Re}} = 2 \times {{10}^{4}}$ (слева) и ${\text{Re}} = 6 \times {{10}^{4}}$ (справа).

На фиг. 3 изображена ведущая часть спектра матрицы $H$ при различных числах Рейнольдса и значениях числа узлов сетки $n = 100$ и $200$. В изображенных ведущих частях спектров кратных собственных значений обнаружено не было. Видно, что имеется сходимость по шагу сетки при $n = 100$. Все дальнейшие результаты расчетов мы будем приводить только для $n = 100$, поскольку расчеты при $n = 200$ давали те же результаты с точностью до приводимых значащих цифр.

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

На фиг. 4 при ${\text{Re}} = 4 \times {{10}^{4}}$ сравниваются абсолютные величины амплитуд компонент нормированной ведущей вещественной моды и нормированных ведущих комплексно-сопряженных мод с абсолютными величинами амплитуд компонент нормированного глобального оптимального возмущения в моменты времени $t = 0$ и ${{t}_{{{\text{opt}}}}}$. Видно, что как ведущая вещественная, так и ведущие комплексные моды существенно отличаются от оптимального возмущения как при $t = 0$, так и при $t = {{t}_{{{\text{opt}}}}}$.

6.2. Спектральный состав глобальных оптимальных возмущений

Для каждого из рассматриваемых значений числа Рейнольдса и соответствующих оптимальных значений волновых чисел будем рассматривать инваринтные подпространства, отвечающие изолированным подмножествам спектра $\lambda (H)$ матрицы $H$ следующего вида:

${{\Lambda }_{1}} = \{ \lambda \in \lambda (H):\operatorname{Real} \lambda > {{r}_{1}}\} ,\quad {{\Lambda }_{2}} = \lambda (H){\backslash }{{\Lambda }_{1}},$

${{\Lambda }_{{11}}} = \{ \lambda \in \lambda (H):\operatorname{Real} \lambda > {{r}_{2}}\} ,\quad {{\Lambda }_{{12}}} = {{\Lambda }_{1}}{\backslash }{{\Lambda }_{{11}}},$

где ${{r}_{1}} < {{r}_{2}} < 0$. Подмножества ${{\Lambda }_{{11}}}$ и ${{\Lambda }_{{12}}}$ выделены на фиг. 3 соответственно красным и синим цветами.

Пусть $\Lambda $ означает одно из описанных выше подмножеств спектра, $\mathcal{U}$ – инвариантное подпространство, отвечающее $\Lambda $, $dim\mathcal{U}$ – его размерность, т.е. суммарная алгебраическая кратность собственных значений, входящих в $\Lambda $. Нас будет интересовать квадрат нормы ${{c}_{\mathcal{U}}}(t)$ проекции глобального оптимального возмущения на $\mathcal{U}$ в моменты времени $t = 0$ и ${{t}_{{{\text{opt}}}}}$ и максимальная амплификация ${{\Gamma }_{{\mathcal{U},max}}}$ на этом подпространстве. Результаты вычисления этих величин для введенных подмножеств спектра приведены в табл. 2.

Таблица 2.

Размерность инвариантного подпространства $\mathcal{U}$, максимальная амплификация ${{\Gamma }_{{\mathcal{U},max}}}$ и квадрат нормы проекции глобального оптимального возмущения на $\mathcal{U}$ в моменты времени $t = 0$ и ${{t}_{{{\text{opt}}}}}$ в зависимости от числа Рейнольдса

Параметры	$\mathcal{U}$	$dim\mathcal{U}$	${{\Gamma }_{{\mathcal{U},max}}}$	${{c}_{\mathcal{U}}}(0)$	${{c}_{\mathcal{U}}}({{t}_{{{\text{opt}}}}})$
$\operatorname{Re} = 2 \times {{10}^{4}}$	${{\mathcal{U}}_{1}}$	21	24.6367	1.0133	24.8382
${{r}_{1}} = - 0.11$	${{\mathcal{U}}_{2}}$	279	1.0000	0.0133	0.0000
${{r}_{2}} = - 0.025$	${{\mathcal{U}}_{{11}}}$	9	24.0670	1.0584	24.8702
	${{\mathcal{U}}_{{12}}}$	12	1.6596	0.0685	0.0004
$\operatorname{Re} = 4 \times {{10}^{4}}$	${{\mathcal{U}}_{1}}$	37	30.3495	1.0216	30.4001
${{r}_{1}} = - 0.24$	${{\mathcal{U}}_{2}}$	263	1.0900	0.0216	0.0000
${{r}_{2}} = - 0.13$	${{\mathcal{U}}_{{11}}}$	26	30.0807	3.2660	30.4001
	${{\mathcal{U}}_{{12}}}$	11	1.0000	2.3251	0.0000
$\operatorname{Re} = 6 \times {{10}^{4}}$	${{\mathcal{U}}_{1}}$	52	38.0686	1.0122	38.0889
${{r}_{1}} = - 0.38$	${{\mathcal{U}}_{2}}$	248	1.1516	0.0122	0.0000
${{r}_{2}} = - 0.15$	${{\mathcal{U}}_{{11}}}$	33	37.8930	10.5533	38.0889
	${{\mathcal{U}}_{{12}}}$	19	1.0000	9.5607	0.0000

Из табл. 2 видно, что глобальное оптимальное возмущение лежит главным образом в подпространстве ${{\mathcal{U}}_{1}}$, поскольку квадрат нормы его проекции на дополнительное инвариантное подпространтво ${{\mathcal{U}}_{2}}$ при $t = 0$ равен величине порядка ${{10}^{{ - 2}}}$, а при ${{t}_{{{\text{opt}}}}}$ становится равным величине порядка ${{10}^{{ - 4}}}$ либо меньше. При этом необходимая размерность подпространства ${{\mathcal{U}}_{1}}$ растет (21, 37 и 52) с ростом числа Рейнольдса.

Учитывая, что ${{\Gamma }_{{U,max}}}$ примерно равно ${{c}_{\mathcal{U}}}({{t}_{{{\text{opt}}}}}){\text{/}}{{c}_{\mathcal{U}}}(0)$ для $\mathcal{U} = {{\mathcal{U}}_{1}}$, можно также заключить, что проекция глобального оптимального возмущения на подпростраство ${{\mathcal{U}}_{1}}$ является в этом подпространстве оптимальным возмущением, т.е. имеет наибольшую амплификацию. В то же время любой вектор из дополнительного подпространства ${{\mathcal{U}}_{2}}$ не может иметь амплификацию, большую максимальной амплификации на этом подпространстве, которая максимальна в случае ${\text{Re}} = 6 \times {{10}^{4}}$ и равна $1.1516$.

Для более детального спектрального анализа глобального оптимального возмущения мы разлагали инвариантное пространство ${{\mathcal{U}}_{1}}$ в прямую сумму двух инвариантных подпространств ${{\mathcal{U}}_{{11}}}$ и ${{\mathcal{U}}_{{12}}}$ и рассматривали проекции глобального оптимального возмущения на каждое из этих подпространств. Результаты некоторых возможных вариантов таких разбиений показаны в остальной части табл. 2. Видно, что подскок глобального оптимального возмущения достигается за счет двух различных факторов. Во-первых, величина квадрата нормы его проекции на подпространство ${{\mathcal{U}}_{{11}}}$ при $t = 0$ больше единицы. Во-вторых, квадрат нормы этой проекции возрастает при $t = {{t}_{{{\text{opt}}}}}$, а проекция глобального оптимального возмущения на подпространство ${{\mathcal{U}}_{{12}}}$ становится малозначимой. При максимальном числе Рейнольдса наиболее значим первый фактор, при минимальном – второй.

Интерес представляет и то, что максимальная амплификация векторов из подпространства ${{\mathcal{U}}_{{11}}}$ немногим меньше амплификации глобального оптимального возмущения. Возникает вопрос, чем отличается глобальное оптимальное возмущение от оптимального возмущения из подпространства ${{\mathcal{U}}_{{11}}}$. Ответ на этот вопрос дает фиг. 5, где сравниваются абсолютные величины компонент этих оптимальных возмущений при $t = 0$ и ${{t}_{{{\text{opt}}}}}$. Видно, что абсолютные величины сравниваемых оптимальных возмущений близки друг к другу в $C$-норме, причем близость увеличивается с ростом числа Рейнольдса, но глобальное оптимальное возмущение – значительно более гладкая функция (кроме его температурной компоненты в пристеночной области) даже при максимальном из рассмотренных значений числа Рейнольдса.

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

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

Авторы благодарны А.В. Глазунову и Е.В. Мортикову за предоставление данных прямого численного моделирования, интерес к данной работе и полезные обсуждения результатов.

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

Drobinski P., Brown R., Flamant P., Pelon J. Evidence of organized large eddies by ground-based doppler lidar, sonic anemometer and sodar // Bound.-Layer Meteor. 1988. V. 88. № 3. P. 343–361.
Глазунов А.В. Численное моделирование устойчиво-стратифицированных турбулентных течений над плоской и городской поверхностями. // Изв. РАН, сер. ФАиО. 2014. Т. 50. № 3. С. 271–281.
Глазунов А.В., Мортиков E.В., Барсков К.В., Каданцев E.В., Зилитинкевич С.С. О слоистой структуре устойчиво-стратифицированных турбулентных течений со сдвигом скорости // Изв. РАН, cер. ФАиО. 2019. Т. 55. № 4. С. 13–26.
Sullivan P.P., Weil J.C., Patton E.G., Jonker H.J., Mironov D.V. Turbulent winds and temperature fronts in large-eddy simulations of the stable atmospheric boundary layer // J. Atmos. Sci. 2016. V. 73. № 4. P. 1815–1840.
Mortikov E.V., Glazunov A.V., Lykosov V.N. Numerical study of plane Couette flow: turbulence statistics and the structure of pressure-strain correlations // Russ. J. Num. Anal. Math. Model. 2019. V. 34. № 2. P.119–132.
Lilly D.K. On the instability of Ekman boundary flow // J. Atmos. Sci. 1966. V. 23. № 5 P. 481–494.
Brown A.R. A secondary flow model for the planetary boundary layer // J. Atmos. Sci. 1970. V. 27. № 5. P. 742–757.
Glazunov A.V., Zasko G.V., Mortikov E.V., Nechepurenko Yu.M. Optimal disturbances of stably stratified turbulent Couette flow // Doklady Physics. 2019. V. 64. № 7. P. 308–312.
Засько Г.В., Глазунов А.В., Мортиков Е.В., Нечепуренко Ю.М. Крупномасштабные структуры стратифицированного турбулентного течения Куэтта и оптимальные возмущения: Препринты ИПМ им. Келдыша. 2019. № 63.
Zasko G.V., Glazunov A.V., Mortikov E.V., Nechepurenko Yu.M. Large-scale structures in stratified turbulent Couette flow and optimal disturbances // Russ. J. Numer. Anal. Math. Model. 2020. V. 35. № 2. P. 37–53.
Бойко А.В., Клюшнев Н.В., Нечепуренко Ю.М. Устойчивость течения жидкости над оребренной поверхностью. М.: ИПМ им. Келдыша, 2016. 123 с.
Nechepurenko Yu.M., Sadkane M. A low-rank approximation for computing the matrix exponential norm // SIAM J. Matrix. Anal. Appl. 2011. V. 32. № 2. P. 349–363.
Schmid P.J., Henningson D.S. Stability and transition in shear flows. Berlin: Springer–Verlag, 2000.
Romanov V.A. Stability of plane-parallel Couette flow // Func. Anal. Appl. 1973. V. 7. P. 137–146.
Butler K.M., Farrell B.F. Optimal perturbations and streak spacing in wall-bounded turbulent shear flow // Phys. Fluids A: Fluid Dynamics. 1993. V. 5. № 3. P. 774–777.
Canuto C., Hussaini M.Y., Quarteroni A., Zang T.A. Spectral methods. Fundamentals in single domains. Berlin: Springer, 2006.
Canuto C., Hussaini M.Y., Quarteroni A., Zang T.A. Spectral methods. Evolution to complex geometries and applications to fluid dynamics. Berlin: Springer, 2007.
Weideman J.A.C., Reddy S.C. A MATLAB Differentiation matrix suite // ACM Transact. Math. Soft. 2000. V. 26. № 4. P. 465–519.
Нечепуренко Ю.М. О редукции линейных дифференциально-алгебраических систем управления // Докл. АН. 2012. Т. 445. № 1. С. 17–19.
Golub G.H., van Loan C.F. Matrix computations. London: The John Hopkins Univer. Press, 1991.

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

Инструменты

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

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

Архивы выпусков Информация о журнале Отправить рукопись в журнал