«МАТЕМАТИЧЕСКАЯ МОДЕЛЬ СВОБОДНОЙ ТУРБУЛЕНТНОСТИ НА ОСНОВЕ ПРИНЦИПА МАКСИМУМА ...»

где CS – эмпирическая константа (константа Смагоринского). В отличие от турбулентной вязкости, подсеточная вязкость зависит не только от параметров отфильтрованного течения (компонент тензора скоростей деформаций), но и от размеров фильтра (от сетки). При измельчении сетки дополнительные по уменьшаются, и решение LES асимптотически стремится к решению DNS. В этом состоит принципиальное отличие метода LES от метода RANS, в котором измельчение сетки приводит лишь к получению «точных» (независящих от сетки) решений уравнений Рейнольдса.

Как правило, значение параметра Смагоринского выбирается из диапазона Cs=0.06-0.28, в частности, Cs=0.1, 0.15 и 0.18 для течения в канале, течения в свободном слое смешения и однородной изотропной турбулентности. Однако, важной особенностью подсеточных моделей является то обстоятельство, что используемого для решения задачи численного алгоритма. Это объясняется тем, что точность разрешения крупномасштабных вихревых структур в LES зависит не только от сетки, но и от свойств метода, в частности, от присущей ему численной диссипации. Иными словами, численная диссипация метода сама по себе играет роль своеобразной подсеточной модели. Если эта диссипация велика, то уменьшена, а если мала, то, наоборот, увеличена. В связи с этим, строго говоря, для каждого численного метода должна проводиться индивидуальная калибровка константы Смагоринского.

Модель Смагоринского обладает значительной диффузией и диссипацией, что позволяет стабилизировать численные расчеты. Один из способов сделать её менее диссипативной состоит в том, чтобы добавлять диссипацию только в тех областях течения, которые характеризуются существенными изменениями величины завихренности и её направления [15] (Selective Smagorinsky, SS). В другой модификации модели Смагоринского – модели структурной функции (Structure Function, SF) подсеточная вязкость вычисляется по формуле [16]:

Она играет существенную роль в теории турбулентности и связывается со спектральной плотностью кинетической энергии турбулентности. Угловые скобки означают осреднение по точкам, для которых r. С вычислительной точки зрения, преимущество модели структурной функции является использование приращений скорости вместо её производных. Модель даёт более точные результаты для переходных течений, чем модель Смагоринского.

Модель Смагоринского дает ненулевое значение подсеточной вязкости на стенке, где существует ненулевой градиент скорости. Вместе с тем, подавление флуктуаций вблизи стенки приводит к тому, что sgs 0 при y 0. Модели WALE [17] и RNG [11] позволяют учесть влияние стенки без введения пристеночных демпфирующих функций.

В динамических моделях (Dynamic Model, DM) для оценки параметра Смагоринского используется информация, содержащаяся в разрешимых масштабах скорости. Помимо фильтрации с сеточным фильтром, используется вторичная фильтрация поля скорости с фильтром, полоса пропускания которого превышает ширину сеточного фильтра. К примеру, в динамической модели Германо [18] параметр Смагоринского находится на каждом шаге из соотношения Здесь Lij Tij ij – компоненты тензора леонардовых напряжений, Tij и ij – компоненты тензоров подсеточных напряжений, соответствующих сеточному Сагоринского, рассчитанное с помощью динамической процедуры, сильно колеблется в пространстве и времени. Специфическая трудность, являющаяся результатом таких колебаний, состоит в том, что вихревая вязкость может стать отрицательной. Это означает перенос энергии от подсеточных масштабов к разрешенным. В принципе, в нестационарном потоке такой процесс может иметь место, приводя обычно к вычислительной неустойчивости. Для устранения этого недостатка на диапазон изменения модельного параметра накладываются дополнительные ограничения, самое простое из которых СD 0.

динамическая модель [19] (Lagrangian Dynamic Model, LDM), где осреднение модельного параметра производится вдоль траектории движения жидкой частицы, локальная динамическая модель [20] (Dynamical Localization Model, DLM), векторная динамическая модель [21] (Vector Dynamic Model, VDM) и прочие. Все эти модели довольно дороги с вычислительной точки зрения, но улучшают исходную модель Смагоринского, делают её более универсальной.

Существуют так же, как и в RANS подходе, дифференциальные подсеточные модели, которые классифицируются по числу дополнительных уравнений, решаемых в дополнение к фильтрованной по пространству системе уравнений Навье-Стокса. Наибольшее распространение получили модели с одним которая связывается с подсеточной вихревой вязкостью при помощи соотношения sgs C ksgs, где C 0.067 – эмпирический коэффициент. Уравнение для подсеточной кинетической энергии получается при помощи вычитания уравнения кинетической энергии. Формально оно совпадает с уравнением для кинетической энергии турбулентности в k модели турбулентности (отличие заключается в форме записи диссипативного члена) Дифференциальные модели позволяют учесть предысторию потока и нелокальные эффекты, а также обратный переход энергии от подсеточных к разрешенным масштабам. Иногда для преодоления вычислительной неустойчивости уравнение для кинетической энергии используется для отключения вихревой вязкости, когда кинетическая энергия падает до нуля. Во многих случаях привлечение дополнительного уравнения не приводит к более точным результатам по сравнению с алгебраическими моделями подсеточной вихревой вязкости, в связи с чем широкого распространения на практике такие модели не получили.

Существуют ещё модели подобных масштабов (Scale Similarity Model, SSM) [22], которые не используют концепцию вихревой вязкости и не требуют предположения о локальном равновесии порождения и диссипации кинетической энергии, которое используется в модели Смагоринского, а также учитывают возможность обратной передачи энергии от подсеточных к разрешенным масштабам. Недостаток моделей такого типа состоит в том, что они не обладают свойством сеточной диссипации. Недостаточная диссипативность приводит к проблемам с устойчивостью вычислительной процедуры.

В работе [23] проводится детальное сравнение характеристик различных подсеточных моделей. Сравнение результатов позволяет сделать вывод о том, что все подсеточные модели имеют приблизительно одинаковую точность, включая модель с постоянной подсеточной вязкостью, а модификации модели Смагоринского не приводят к значительному улучшению численных оценок.

Квазипрямое численное моделирование Существуют также подходы, которые не используют подсеточные модели в явном виде. Это т.н. неявное моделирование крупных вихрей (Implicit LES, ILES).

Ещё выделяют монотонное интегрирование крупных вихрей – MILES.

Применяемые в ILES/MILES разностные схемы уже имеют достаточную диффузию и диссипацию, которая имитирует стабилизирующее действие подсеточных моделей. В работе [24] показано, что численная диссипация при явном LES моделировании может быть одного порядка с диссипацией подсеточной модели, и вообще может маскировать её. Это значит, что усилия по совершенствованию вычислительных моделей турбулентности можно (может быть, даже с большим успехом) сосредоточить на разработке разностных схем для уравнений Навье-Стокса с улучшенными диссипативными и дисперсионными свойствами.

Впервые техника MILES была предложена в 1992 году группой авторов во главе с Boris J.P. [25]. Идея заключалась в следующем: так как турбулентность характеризуется высоким уровнем флуктуаций завихренности, и возникающие большие градиенты завихренности чем-то схожи с сильными разрывами в сжимаемой жидкости, то для её моделирования можно попытаться использовать т.н. схемы высокой разрешающей способности (shock-capturing schemes), использующиеся для моделирования невязких течений с ударными волнами.

Одним из способов построения таких схем является использование нелинейных разностных схем, например, TVD-схем с «лимитером». Фильтрация «подсеточных» пульсаций осуществляется в них неявно при нелинейной коррекции, без введения полуэмпирических моделей. Дальнейшее развитие TVDподхода привело к созданию популярных на данный момент схем ENO, WENO, которые способны обеспечить достаточную монотонность решения при произвольном порядке точности в областях гладкости.

Следует отметить, что несмотря на вычислительную эффективность и имеющиеся примеры использования [26], подходы ILES/MILES не имеют должного физического обоснования, а полученные результаты нуждаются в тщательной проверке. В существенной степени успех применения подхода определяется удачным выбором разностной схемы для дискретизации конвективных потоков. Например в статье [27] авторы (Eric Garnier et. al.) проводят численное моделирование свободно затухающей однородной изотропной турбулентности с использованием популярных на данный момент эйлеровских схем высокой разрешающей способности (Jameson, MUSCL-TVD, ENO, WENO, MENO). Цель расчетов заключается в оценке возможности использования данного класса схем для моделирования крупных вихрей. По удовлетворять одному из требований: или численная диссипация схемы должна быть на много меньше физической «подсеточной» диссипации (C1), или численная диссипация схемы способна «иммитировать» «подсеточную» (sgs) модель (C2).

Для оценки диссипации схем высокой разрешающей способности свободная турбулентность моделировалась без молекулярной вязкости и явных sgs-моделей, на грубых сетках 643 и 1283 (MILES). Численная диссипация интерпретировалась в терминах sgs-диссипации. Был введен эффективный эквивалент константы в модели Смагоринского – «обобщенная константа Смагоринского» C gs :

где sgs ui j ij u – диссипация в модели Смагоринского, ij u 2 sgs Sij u – тензор подсеточных напряжений, sgs 2 S u – коэффициент вязкости.

S u 2 Sij u Sij u. Как оказалось, для всех вышеперечисленных схем вне зависимости от времени «обобщенная константа Смагоринского» значительно превосходила классическую константу Смагоринского Cs 0.18 0.2 (рис. 0.1, взято из статьи).

Рис. 0.1. Обобщенная константа Смагоринского в зависимости от времени в MILES моделировании однородной изотропной турбулентности на грубых сетках 643, При дополнительном включении sgs-модели, на фоне собственной численной диссипации схем, «подсеточная» диссипация была незначительна и не оказывала существенного влияния на получаемые решения. Наличие столь сильной диссипации оказывало существенное негативное влияние на малые масштабы турбулентности. Это отражалось на энергетических спектрах (сильный загиб спектральной кривой вниз в высокочастотной области) (рис. 0.2).

Рис. 0.2. Энергетические спектры в MILES моделировании однородной изотропной Проанализировав результаты, авторы пришли к выводу: для эйлеровых схем высокой разрешающей способности не выполняется ни первое требование C ( num sgs ), ни второе C2 ( num не имитирует sgs ). Подход MILES, с использованием вышеперечисленных схем высокой разрешающей способности, «червеобразную» вихревую структуру и колмогоровский наклон k 5/3 на некотором участке спектра при достаточном разрешении сетки. Однако значительная часть течения испытывает сильное численное демпфирование.

Особенно это сказывается на малых структурах, определенных менее чем на узлах сетки, для всех схем без исключения.

Концепция идеального LES алгоритма Таким образом, диссипативно-устойчивые разностные схемы могут достаточно хорошо описывать крупные вихри и ослаблять или полностью гасить ненужные мелкомасштабные возмущения. Вместе с тем, возможна ситуация, когда для данного масштаба l диссипация так велика, что вихрь сам гасится схемой. Поэтому схемы с большой диссипацией требуют аккуратного выбора сеток. На рисунке 0.3 представлены возможные виды спектров турбулентности в методах LES/ILES в сравнении c DNS, которые можно получить при использовании схем с большой диссипацией (в), с недостаточной диссипацией (а), с оптимальной диссипацией (г). Также можно вообще не получить правильный наклон инерционного интервала (б).

Рис. 0.3. Возможные виды энергетических спектров турбулентности В случае г) диссипация оптимальна в том смысле, что каскад энергии, направленный от структур с большими размерами к структурам с меньшими размерами, уводится на подсеточный уровень наиболее правильно, без чрезмерной диссипации и без накопления энергии в коротковолновой части. Если при измельчении расчетных сеток, такой алгоритм, поэтапно включая новые гармоники, будет сохранять это свойство (рис. 0.4), то можно говорить о неком идеальном LES (Perfect LES, PLES) методе.

Рис. 0.4. Спектры турбулентности на последовательности сеток при моделировании с помощью Степень близости к идеальному LES была бы неплохим критерием при построении модели турбулентности. В данной работе этот критерий взят за основу.

Цели и структура диссертационной работы Основной целью данной диссертационной работы является разработка математической модели, относящейся к классу Implicit LES, для расчета течений со свободной турбулентностью. Исходной точкой для построения математической модели стал выбор базового вычислительного алгоритма (разностной схемы).

Основным требованием к исходному алгоритму стала его минимальная внутренняя (схемная) диссипация и максимально компактный вычислительный шаблон.

Диссипативность в схеме может достигаться, например, при использовании ориентированных разностей в соответствии с направлением потока, как в схеме симметричных схемах диссипация отсутствует, а вместе с ней отсутствует и механизм ослабления вредных для описания крупных вихрей коротковолновых ошибок. В данной работе исследуется оба способа аппроксимации конвективных членов. Направленная аппроксимация используется в схеме «КАБАРЕ» [28-30] («Upwind Leapfrog»), а центральная в схеме «двухслойный крест» [31, 32] («Central Leapfrog»). Схема «двухслойный крест», является баланснохарактеристическим представлением классической схемы «Крест» с разнесенными консервативными и потоковыми переменными, как в схеме «КАБАРЕ». Преимущество такой модификации схемы «Крест» заключается в компактности её шаблона и в возможности введения коррекции потоков. Обе пространственно-временных сетках и являются явными, максимально компактными, консервативными, бездиссипативными и обратимыми по времени. Помимо схем «КАБАРЕ» и «двухслойного креста» в работе исследуются схемы, полученные линейной гибридизацией схем «КАБАРЕ» и «Крест», с параметром гибридизации 0 1 ( 0 – «КАБАРЕ», 1 – «Крест»), изменяя который можно регулировать диссипацию и дисперсию в схеме.

Далее, нужно было сконструировать сеточный диссипативный механизм, который сохранял бы статистические характеристики турбулентных течений в инерционном интервале для вихрей всех размеров, представимых на заданной расчетной сетке. Исходным было предположение, что таким механизмом может оказаться нелинейная коррекция потоков на основе принципа максимума [29], предложенная ранее для схемы КАБАРЕ для обеспечения монотонности газодинамических течений. И оно в полной мере подтвердилось.

Диссертация состоит из введения, четырёх глав и заключения.

Первая глава содержит результаты спектрально анализа исследуемых схем в линейном и нелинейном случаях на примере простейшего одномерного уравнения переноса. В линейном случае (без коррекции) анализ диссипативных и дисперсионных ошибок проводился аналитически, а при наличии коррекции в качестве оценки использовался численный метод. Далее рассматривается нелинейное уравнение Бюргерса и задача об одномерной турбулентности (бюргюленции). Результаты моделирования на грубых сетках по схемам «КАБАРЕ», «двухслойный крест» и их гибридам сравниваются с DNS расчетом и с результатами по другими схемами («WENO-5» (MILES) и «S2-SMA» (LES)).

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

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

Четвертая глава посвящена трехмерным вихревым течениям.

Рассматриваются устойчивый колоннообразный вихрь Рэнкина и неустойчивые вихри Хилла и Тейлора-Грина. Результаты моделирования вихря Тейлора-Грина по всем рассматриваемым схемам сравниваются с DNS расчетом Брэчета (1991).

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

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

Глава1. Спектральный анализ Схемы «КАБАРЕ» и «двухслойный крест» для одномерного уравнения переноса Пусть задана одномерная область G x : x 0; L, покрытая равномерной сеткой h xi : i 1,..., N, где xi – координаты узлов сетки. Одномерное уравнение переноса t c x 0, c 0, можно аппроксимировать с помощью центральных разностей по времени и пространству Эта схема, известная как схема «Крест» (в иностранной литературе схема Leapfrog [32]), имеет второй порядок аппроксимации по пространству и времени в силу симметричной структуры её шаблона (рис. 1.1. (a)). При числах Куранта CFL 1 схема устойчива и бездиссипативна.

Рис. 1.1. Шаблон схемы «крест» (a) – классический, (b) – с введением дополнительных Одним из недостатков схемы «крест» является отсутствие компактности её шаблона (схема трехслойная по времени и определена на двух пространственных ячейках). Для устранения этого недостатка в работе [31] была предложена двухслойная консервативная модификация схемы «Крест».

Введем новые переменные в центрах пространственно-временных ячеек (рис. 1.1. (b)) и аппроксимируем уравнение переноса так:

Если сложить (1.2) и (1.3) переменные и между собой то получим классический трехслойный «Крест» (1.1).

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

На неравномерных пространственных сетках процедура экстраполяции (1.5) дает только первый порядок аппроксимации. Для сохранения второго порядка необходимо подправить её следующим образом:

Чтобы использовать такую схему на неравномерных временных сетках, введем ещё одни дополнительные переменные в центрах ячеек на целых временных слоях in1/2 и перепишем «двухслойный крест» в трехэтапной форме:

Если сложить первое и третье уравнение системы (1.7), то мы получим консервативную аппроксимацию уравнения переноса с симметричным шаблоном в пределах одной пространственно-временной ячейки (рис 1.2):

Отличие схемы «двухслойный крест» от схемы «КАБАРЕ» [28, 29, 33] возникает только во второй фазе алгоритма. Если в схеме «двухслойный крест»

экстраполяция во второй фазе центральная (Central Leapfrog) (1.6), то в схеме «КАБАРЕ» направленная (Upwind Leapfrog):

Сравнение диссипативных и дисперсионных свойств схем «КАБАРЕ» и «двухслойный крест»

Проанализируем диссипативные и дисперсионные свойства схемы (1.7). Для этого представим решение в виде бегущих волн (для равномерной сетки):

где A, B и C – произвольные амплитуды.

Подставляя частное решение (1.10) в (1.7), приходим к системе уравнений:

Для того чтобы эта система имела решение при любых амплитудах A, B, D нужно, чтобы ее определитель был равен нулю:

Решая полученное уравнение, получаем:

Если подкоренные выражения неотрицательны, то легко видеть, что Получим условие не отрицательности подкоренного выражения:

Т.к. max sin kh 1, то подкоренное выражение будет не отрицательным при r 1. Т.о. при числах Куранта CFL 1 схема «двухслойный крест», как и схема классический «Крест», устойчива и бездиссипативна, что, очевидно, учитывая, что новая схема это всего лишь компактная форма записи схемы «Крест».

На рисунке 1.3 приведены дисперсионная / kc и диссипативная q поверхности схемы «двухслойный крест» для модуля перехода q1. Второй корень q2 имеет «паразитный» характер и проявляет себя только при несогласованном задании начальных данных 0 и 0.

Рис. 1.3. Дисперсионная (a) и диссипативная (b) поверхности схемы «двухслойный Аналогичные результаты для схемы «КАБАРЕ» [33] приведены на рисунке 1.4.

Рис. 1.4. Дисперсионная (a) и диссипативная (b) поверхности схемы «КАБАРЕ»

В схеме «КАБАРЕ» дисперсионная поверхность имеет особенности при числе Куранта CFL 0. Дисперсионная поверхность схемы «двухслойный крест»

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

на месте. Попытки оптимизировать дисперсионные свойства расчетной схемы естественным образом привели к построению гибридных схем как линейной комбинации схем, использующих центральную («Крест») и направленную («КАБАРЕ») процедуру экстраполяции при вычислении потоковых переменных.

Гибридные схемы Введем дополнительный параметр 0 1 – параметр гибридизации.

Вторую фазу гибридной схемы будем вычислять как линейную комбинацию центральной и направленной экстраполяций:

Построим дисперсионные и диссипативные поверхности схемы при разных значениях параметра гибридизации (рис. 1.5).

Рис. 1.5. Дисперсионные (слева) и диссипационные (справа) поверхности гибридной Гибридные схемы при 0 1 устойчивы при числах Куранта CFL 1, однако они уже не являются бездиссипативными, как схемы «КАБАРЕ» и «Крест». Максимум схемной диссипации достигается при 0.5.

При 0 особенность в дисперсионной поверхности, присущая схеме «КАБАРЕ», пропадает. На рисунке 1.6 приведено сравнение профилей дисперсионных и диссипативных поверхностей по схемам «КАБАРЕ», «двухслойный крест», «hybrid_01», «hybrid_03», «hybrid_05», «hybrid_07» и «hybrid_09» для трех чисел Куранта CFL 0.1, CFL 0.3 и CFL 0.8.

Рис. 1.6. Дисперсионные (слева) и диссипационные (справа) кривые для схем «КАБАРЕ», Нелинейная коррекция потоковых переменных Согласно теореме Годунова [34] среди линейных разностных схем для уравнения переноса с порядком аппроксимации выше первого нет монотонных.

Для улучшения монотонности решений можно ввести дополнительную процедуру коррекции потоковых переменных как в схемах высокой разрешающей способности. Это могут быть различные TVD-лимитеры, фильтры высоких частот, или, как введенная для балансно-характеристических схем, нелинейная коррекция потоковых переменных на основе принципа максимума [28, 35].

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

Анализ диссипативных и дисперсионных свойств нелинейных схем В статьях [36, 37] был развит метод для изучения дисперсионных и диссипативных свойств нелинейных схем высокой разрешающей способности (TVD, WENO и т.п.). Воспользуемся данным методом для оценки реальной диссипации в схемах при наличии коррекции по принципу максимума.

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

где u – переменная, c u – поток, c 0. Решением задачи (1.20) будет являться u x, t u0 x c t. Тогда разложение в ряд Фурье дает способ нахождения модифицированного волнового числа k :

Где u – численное решение, – фаза решения. Для спектральных схем k k.

диссипативных и дисперсионных ошибок, связанных с конечно-разностной дискретизацией схемы.

На серии рисунков 1.7 показаны зависимости действительной части модифицированного волнового числа (дисперсионные характеристики) для чисел Куранта 0.1, 0.3 и 0.8. А также проведено исследование на сходимость при увеличении количества шагов по времени nt 5, 10, 20...

Re(k`)/kmax характерно наличие аномальной дисперсии (фазовые скорости гармоник превышают скорость распространения возмущений c). В нелинейном случае аномальная дисперсия также возникает при малых числах Куранта CFL 0.3 (рис.

1.7). При числах Куранта CFL ~ 0.3 все расчетные точки с высокой точностью ложатся на прямую k k. При больших значениях чисел Куранта схема «КАБАРЕ» дает нормальную дисперсию.

модифицированного волнового числа (диссипативные характеристики) для чисел Куранта 0.1, 0.3 и 0.8. А также проведено исследование на сходимость при увеличении количества шагов по времени.

Im(k`)/kmax Наибольшая диссипация, обусловленная действием принципа максимума, соответствует высокочастотным гармоникам и малым числам Куранта.

Аналогичные результаты для схемы «двухслойный крест» с коррекцией приведены на рисунках 1.9-1.10.

Re(k`)/kmax Рис. 1.9. Дисперсионные кривые для схемы «двухслойный крест» с коррекцией по является нормальной (фазовая скорость меньше скорости распространения возмущений).

Im(k`)/kmax Рис. 1.10. Диссипативные кривые для схемы «двухслойный крест» с коррекцией по то для гибридной схемы будем строить только диссипативные характеристики. На рисунке 1.11 приведены диссипативные характеристики для схемы «hybrid_05».

Im(k`)/kmax Рис. 1.11. Диссипативные кривые для схемы «hybrid_05» с коррекцией по принципу Im(k`)/kmax Рис. 1.12. Диссипативные кривые для схемы «hybrid_09» с коррекцией по принципу «WENO5-LF-RK3» [38, 39], являющийся очень популярной в настоящее время схемой высокой разрешающей способности, и со схемой «центральная разность»

второго порядка «CD2_020» с вязкостью фон Неймана-Рихтмайера [40].

Коэффициент вязкости взят равным 0.2, как в модели Смагоринского для расчета однородной-изотропной турбулентности [1, 14]. На рисунке 1.13 приводятся результаты расчета дисперсии с числами Куранта CFL=0,1 и CFL=0,3.

CABARET CABARET

Im(k`)/kmax Рис. 1.14. Сравнение диссипативных кривых по схемам «КАБАРЕ», «двухслойный крест», «hybrid_05», «hybrid_09», «WENO5_LF_RK3», «CD2_020». Число Куранта Таким образом, можно сделать вывод:

для длинных и средних волн k 0.6 kmax (3-4 узла сетки на длину волны) наиболее диссипативна схема «WENO5-LF-RK3», а для высоких частот гибридная схема «hybrid_05». Наименее диссипативна схема «двухслойный крест». Если сравнить её диссипацию с вязкой диссипацией в LES-модели Смагоринского («CD2_020»), то она примерно вдвое меньше. Подбирая параметр гибридизации, можно добиться того, чтобы диссипация в гибридной схеме была бы примерно равна диссипации LES-модели Смагоринского. Это выполняется в схеме «hybrid_09» (рис. 1.14);

дисперсия меньше всего в схеме «КАБАРЕ», а самая большая, наоборот, в схеме «двухслойный крест». Дисперсия в схеме «hybrid_09» примерно совпадает с дисперсией «CD2_020».

Перенос профиля на неравномерной сетке Для исследования транспортных свойств рассматриваемых схем на неравномерных пространственных сетках зададимся степенью неравномерности, равной 1 / 4, и сгенерируем случайную, распределенную по нормальному закону, последовательность отрезков hi размер области, N – число узлов. Координаты узлов будем определять из соотношений x1 0, xi 1 xi hi, i 1, N 1. Чтобы выполнялось соотношение определенной таким образом расчетной сетке (рис. 1.15) будем решать задачу переносимой функции определим следующим образом [33]:

На серии рисунков 1.16 приведены результаты расчетов переноса начального профиля для числа Куранта CFL=0,1 через 10 и 50 пролетов.

Рис. 1.16. Перенос начального профиля на неравномерной сетке для числа Куранта CFL=0.1 без коррекции (слева) и с коррекцией (справа) Анализируя полученные результаты можно отметить:

наименьшие амплитудные ошибки получаются по схеме «hybrid_05», причем коррекция заметного влияния не оказывает, а наименьшие фазовые ошибки в схеме «КАБАРЕ»;

коррекция по принципу максимума не делает решения монотонными для гибридных схем и схемы «двухслойный крест»;

Уравнение Бюргерса Широко распространено мнение, что одномерным аналогом уравнений Навье-Стокса является уравнение Бюргерса:

Действительно, их роднит квадратичная нелинейность в адвективном слагаемом и линейная вязкость в правой части. Что касается свойств решений, то они совершенно разные. У уравнения Бюргерса, при коэффициенте вязкости, стремящемся к нулю, формируются как сильные (ударные волны), так и слабые разрывы, в то время как решения уравнений Навье-Стокса такими особенностями не обладают. Уравнение (1.24) с помощью преобразования Хопфа сводится к линейному уравнению параболического типа, что позволяет проанализировать все особенности решения при произвольных начальных данных. Отсюда, в частности, следует, что все стохастические свойства ансамбля решений полностью определяются стохастическими свойствами начальных данных.

Тем не менее, уравнение Бюргерса является объектом пристального внимания исследователей, изучающих свойства турбулентности [41-43]. Задача, обычно ставится следующим образом: на отрезке x 0,2 задаются случайные начальные данные с заданными статистическими свойствами, и решается нестационарное уравнение (1.24) с периодическими граничными условиями и заданной случайной функцией f.

Начальные данные задаются рядом Фурье с заданными коэффициентами ck и случайными фазами k, равномерно распределенными на отрезке 0,2 :

Коэффициенты ck вычисляются как:

где A – нормировочный множитель,, k0 – заданные константы.

Целью расчетов является изучение поведения во времени некоторых средних по ансамблю величин, например, полной кинетической энергии, спектральных функций этих величин и т.н. структурных функций случайного поля скоростей. Если внешняя случайная сила исследуется затухание однородной «бюргуленции», в противном случае – влияние «форсинга» (от англ. forcing) на статистически установившееся состояние. Известно, что коэффициент наклона спектральной кривой на логарифмической плоскости в одномерном случае равен не «-5/3», как в трехмерном случае, а «-2». Это обстоятельство обычно используется для тестирования трехмерных LES алгоритмов.

Прежде чем переходить к моделированию «бюргюленции», рассмотрим задачу Коши для закона сохранения t u x u 2 / 2 0 с периодическими граничными условиями и начальным профилем в виде прямоугольника («ступенька»). Точное решение состоит из ударной волны распространяющейся направо и линейной волны разрежения с левой стороны «ступеньки». Схема «КАБАРЕ» и гибридные схемы легко обобщаются на случай переменной скорости переноса [28]. На рисунке 1.17 приведены результаты расчета по всем схемам для числа Куранта CFL 0.3.

Рис. 1.17. Расчет «ступеньки» для числа Куранта CFL=0.3 без коррекции (слева) и с Анализируя полученные результаты, можно сделать следующие выводы:

схема «КАБАРЕ» с коррекцией монотонна и размазывает фронт ударной волны на 1 ячейку;

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

Для сравнения приведем результаты по схемам «WENO5-LF-RK3» (рис. 1. (a)) и полностью консервативной кососимметричной центральной схеме второго порядка «S2» [44]. Схема «S2» сохраняет не только полный импульс, но и полную энергию системы, и поэтому является устойчивой без введения искусственной вязкости. Разностная аппроксимация уравнения Бюргерса в схеме «S2» выглядит следующим образом:

Пространственная аппроксимация схемы «S2» в линейном случае, совпадает диссипативных и дисперсионных свойств схемы «S2» будет совпадать со схемой «CD2». Т.к. схема второго порядка, то она не монотонна (рис. 1.18 (b)), поэтому введём на нижнем временном слое вязкость фон Неймана-Рихтмайера с коэффициетом 0.2 (как в модели Смагоринского) «S2-020» (рис. 1.18 (c)).

Рис. 1.18. Расчет «ступеньки» для числа Куранта CFL=0.3 по схеме «WENO5-LF-RK3»

DNS моделирование затухающей «бюргюленции»

Перейдём теперь к рассмотрению задачи о затухающей «бюргюленции»

(1.24)-(1.25). Прямое численное моделирование подразумевает, что сеткой разрешаются все масштабы, включая колмогоровский масштаб l. Оценим одномерном случае роль диссипации энергии играет скорость диссипации максимальной вариации скорости:

имеющая размерность м/с2. Будем считать скорость диссипации константой, не зависящей ни от времени, ни координаты. Кинематическая вязкость имеет размерность м2/с, следовательно колмогоровский масштаб из размерностных соображений равен l 1/3 2/3. Пусть начальная вариация скорости равна u0.

Характерный масштаб времени T L / u0. Тогда оценка для скорости диссипации u0 / L. Таким образом, коэффициент кинематической вязкости:

где константа порядка единицы.

Возьмем расчетную сетку из 4096 ячеек. Пусть колмогоровский масштаб Нормировочный множитель определится из условия нормировки на энергию E u0 4096, где u0 – среднеквадратичная величина пульсаций скорости начального распределения. Пусть u0 0.1. Заданному начальному спектру распределений пульсаций скорости. Пример одного случайного распределения пульсаций с заданным спектром и его эволюция с течением времени представлена на рисунке 1.19.

Рис. 1.19. Эволюция поля скорости в задача о моделировании «бюргюленции». Схема Энергетические спектры турбулентных пульсаций для DNS-расчета по схемам «КАБАРЕ», «двухслойный крест», «hybrid_05», «S2», «WENO5-LF-RK3», осредненный по 100 реализациям, приводится на рисунке 1.20.

Рис. 1.20. Энергетический спектр «бюргюленции» по всем рассматриваемым схемам Спектральная кривая сохраняет наклон «-2» в инерционном интервале частот k 400 (примерно 10 точек на волну). При частотах k 400 спектр попадает в диссипативный интервал и загибается вниз. Различие между схемами проявляется только на высоких гармониках при k 500 (рис. 1.20 последний график). Наиболее диссипативна схема «WENO5-LF-RK3», наименее – схема «двухслойный крест». Далее, в качестве «эталонного» DNS будем использовать расчет по схеме «S2».

Расчет «бюргюленции» на грубых расчетных сетках На рисунке 1.21 приведены расчеты «бюргюленции» на грубых сетках, не разрешающих колмогоровский масштаб ( N 32, 64, 128, 256, 512 ).

Рис. 1.21. Энергетический спектр «бюргюленции» на последовательности грубых Результаты по схеме «WENO5-LF-RK3» приведены на рисунке 1.22.

Рис. 1.22. Энергетический спектр «бюргюленции» на последовательности грубых Результаты по схеме «S2-020» приведены на рис. 1.23.

Рис. 1.23. Энергетический спектр «бюргюленции» на последовательности грубых Результаты по схеме «двухслойный крест» приведены на рисунке 1.24.

Рис. 1.24. Энергетический спектр «бюргюленции» на последовательности грубых Из рисунков видно, что вязкости в схеме «двухслойный крест» не хватает для отвода на подсеточный уровень турбулентных пульсаций. Гибридные схемы более диссипативны. Приведём пример для схемы «hybrid_05» (рис. 1.25).

Рис. 1.25. Энергетический спектр «бюргюленции» на последовательности грубых Для гибридной схемы «hybrid_03» результаты максимально приближены к DNS (рис. 1.26).

Рис. 1.26. Энергетический спектр «бюргюленции» на последовательности грубых Таким образом, на грубых сетках, диссипация в схеме «КАБАРЕ», подсеточный уровень ровно столько энергии, сколько её приходит с более низких гармоник, как это и должно происходить в инерционном интервале спектра. Как в PLES алгоритме, нет ни накопления энергии на высоких частотах в окрестности частоты Найквиста, ни чрезмерного подавления высоких гармоник.

Схема «КАБАРЕ» и некоторые её гибридные модификации не уникальные в этом плане. Например, если в схеме «WENO» проводить интерполяцию из центров ячеек в узлы не функции потока ( f 0.5u 2 ), а сами скорости (u ), а потом восстанавливать по скоростям поток, то можно получить спектры, практически такие же, как и по схеме «КАБАРЕ» (рис. 1.27). Назовём эту схему «WENO5-GRK3», где G означает, что при решении частичной задачи Римана использовался поток Годунова [45], а не поток Лакса-Фридрихса, как в оригинальных статьях разработчиков схемы «WENO». Т.е. можно утверждать, что множество PLES не пусто.

Полученные результаты говорят в пользу дальнейшего построения математической модели свободной турбулентности на базе разностной схемы «КАБАРЕ» в двумерном и трехмерном случаях.

Рис. 1.30. Энергетический спектр «бюргюленции» на последовательности грубых Глава 2. Турбулентность (краткий обзор) Модель Навье-Стокса несжимаемой жидкости Опыт показывает, что в довольно широком классе течений многих жидкостей даже большие изменения давления не приводит к существенному изменению плотности. Поэтому в таком классе плотность можно считать константой. Давление перестает быть термодинамическим параметром состояния, поскольку перестает участвовать в основном термодинамическом тождестве:

const V 1 / const dV 0, и слагаемое pdV исчезает. Поэтому, основное термодинамическое тождество принимает вид свидетельствующий о том, что приток тепла в среду идет только на увеличение ее внутренней энергии. Следовательно, тепловые потоки в среде не влияют на движение среды и могут быть найдены уже после нахождения u и p. Это позволяет выделить уравнение притока тепла из модели и решать его независимо.

Если в задаче не требуется находить тепловые характеристики, то его можно отбросить. Далее, так как const, уравнение неразрывности принимает вид div u 0, поэтому объемная вязкость перестает играть какую-либо роль в модели. Таким образом, остается только сдвиговая вязкость. В общем случае коэффициент кинематической вязкости может довольно сильно зависеть от температуры T (например, в магматических расплавах вязкость в зависимости от температуры может отличаться на несколько порядков). Однако в простейшей модели, рассматриваемой нами, мы будем считать, что const.

Такое ограничение оставляет класс описываемых жидкостей достаточно широким.

Уравнение неразрывности и уравнение сохранения импульса составляют математическую модель вязкой несжимаемой жидкости:

Эта система уравнений называется системой уравнений Навье–Стокса и представляет собой систему из четырех скалярных уравнений для четырех скалярных неизвестных u и p. Модель Навье-Стокса – одна из наиболее широко применяющихся моделей жидкости. Нахождение общего аналитического решения системы уравнений Навье-Стокса осложняется тем, что оно нелинейное и сильно зависит от граничных и начальных данных. В настоящее время решения найдены лишь в некоторых частных случаях с простой геометрией (течение Пуазейля, используется численное моделирование.

Остановимся теперь на интегралах движения [46], т.е. величинах, сохраняемых уравнениями при невязкой эволюции. Уравнение движения (2.1) перепишем в переменных Лагранжа:

Умножим уравнение (2.2) скалярно на u и проинтегрируем по объему.

С учетом уравнения неразрывности div u 0, первый интеграл в правой части равен нулю:

Второй интеграл преобразуем так:

энстрофию (интеграл от квадрата завихренности по всему объёму). Тогда закон сохранения кинетической энергии примет вид:

Т.о., скорость диссипации кинетической энергии определяется произведением вязкости на энстрофию.

Выпишем теперь уравнение, описывающее эволюцию энстрофии. Если учесть векторное тождество u, u rot на уравнение движения, то получим:

a, b, c b a, c c a, b, приходим к т.н. уравнению Гельмгольца:

При отсутствии вязкости в уравнении (2.8) завихренность как бы вморожена в жидкость, перемещаясь и деформируясь вместе с ней. Это следует из аналогии уравнения (2.8) и уравнения описывающее расстояние между двумя соседними бесконечно близкими точками в жидкости dt l l, u. Можно переписать правую часть (2.8) в виде:

антисимметричный тензор вращения. Подставляя полученное выражение в (2.8):

Умножим уравнение (2.10) на завихренность и проинтегрируем по объёму:

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

В трехмерных уравнениях Навье-Стокса имеется еще один инвариант – спиральность, введенный Моффатом в 1969 году [47]:

Спиральность характеризует степень связности вихревых линий в потоке [48].

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

Умножая (2.8) на скорость u, не трудно получить уравнение баланса для спиральности:

Каскад энергии Для описания баланса энергии в одном отдельно взятом масштабе, нужно записать уравнение Навье-Стокса в пространстве Фурье. Пусть течение занимает ограниченное пространство, затухая на бесконечности, и все входящие в уравнение Навье-Стокса величины допускают представление в виде интеграла Фурье [46]:

где Выразим величины в уравнении (2.1) через Фурье образы:

Уравнение (2.16) умножается на eikr и интегрируется по dr. Учитывая, что Уравнение неразрывности в пространстве Фурье имеет простой вид и может быть использовано для исключения из уравнения (2.17) члена с давлением. После простых преобразований приходим к уравнению Далее, умножим уравнение (2.19) на u k :

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

Т.к. целью проводимых преобразований является уравнение для энергии, заключенной в данных масштабах, которая получается путем интегрирования квадрата модуля Фурье-компонент поля скорости по всем волновым векторам с заданным значением модуля k k :

то проинтегрировав уравнение (2.20) по поверхности сферы, радиуса k, получим следующую структуру, описывающую баланс энергии в системе:

Здесь член T k получается из нелинейного слагаемого в уравнении (2.20) и описывает перенос энергии в заданный масштаб k в результате тройных взаимодействий пульсаций скорости, D k 2 k 2 E k описывает диссипацию за счет действия молекулярной вязкости, и F k – форсинг, характеризует приток энергии за счет сил, поддерживающих турбулентное движение (работа внешних сил).

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

Рассмотрим величину E k E k dk, равную энергии, заключенной во всех масштабах, больших данного. Соответствующее уравнение получается путем интегрирования уравнения (2.22) от 0 до k :

Если рассмотреть масштаб, принадлежащий инерционному интервалу, и считать течение стационарным, то k – скорость спектрального переноса, по сути, есть поток энергии через текущий масштаб k. Этот поток равен суммарной энергии вносимой в течение за единицу времени, а также он равен суммарной скорости диссипации энергии, т.е.

энергии передающейся в тепло за единицу времени.

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

1-я гипотеза состоит в том, что статистические свойства в инерционном и диссипативном интервале не зависят от способа возбуждения турбулентности и диссипации энергии, кинематической вязкостью и самим масштабом k.

2-я гипотеза состоит в том, что статистические свойства турбулентности в инерционном интервале универсальны и зависят только от скорости диссипации энергии и масштаба k.

Далее делается самое сильное предположение, являющееся, по сути, главной гипотезой теории К41. Оно заключается в том, что скорость диссипации энергии является универсальной константой для заданного однородного и изотропного турбулентного течения, т.е. в любой момент времени и в любой точке пространства диссипация энергии в единицу времени равна.

Соображения размерности позволяют получить форму спектра пульсаций скорости. Если спектр энергии может зависеть только от и k, то единственно возможная комбинация есть Формулу (2.26) называют законом Колмогорова, а входящую в него безразмерную константу C – константой Колмогорова.

Очевидно, закон подобия k 5/3 обеспечивает постоянство потока энергии в инерционном интервале k const (действительно: если k const и из соображений размерности E 3/2 k 5/2, то E k 5/3 ).

Теория KLB (Kraichnan-Leith-Batchelor) двумерной турбулентности Турбулентность – это явление существенно трехмерное, и переход к плоской геометрии приводит к качественным изменениям свойств течений [46].

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

Несмотря на это, двумерная турбулентность привлекла к себе внимание значительного числа исследователей, которое не ослабевает и по сегодняшний день. Объясняется это несколькими причинами. Во-первых, качественное своеобразие двумерной турбулентности дает прекрасные возможности для опробования различных моделей турбулентности (модель, претендующая на адекватное описание турбулентности должна быть чувствительной к изменению размерности пространства и правильно отражать её свойства в случае трех и двух измерений). Во-вторых, двумерная турбулентность стала доступной для прямых численных экспериментов уже в 70-х годах (в 80-х с появлением ЭВМ типа «Cray» удалось выйти на сетки 1024х1024, достаточные для приличного воспроизведения инерционных интервалов), а такое разрешение для трехмерных потоков стало возможным только в последние годы. Третья причина состоит в том, что, хотя строго двумерных турбулентных течений в природе не существует, многие черты двумерной турбулентности проявляют крупномасштабные геофизические и астрофизические течения.

Крэйчнан [49], Лейтс [50] и Бэтчелор [51] были пионерами в создании теории двумерной турбулентности. Теория KLB двумерной турбулентности, ставшая классической, основана на статистическом подходе, развитым Колмогоровым при изучении свойств мелкомасштабной трехмерной турбулентности.

Главное отличие в свободной эволюции двумерной турбулентности от эволюции трехмерной турбулентности следует из анализа уравнения (2.8). В двумерном случае механизм растяжения вихрей, который обеспечивал прямой энергетический каскад и рост энстрофии в трехмерном течении, блокирован.

Энстрофия в невязком случае является вторым инвариантом движения:

Это приводит к тому, что, если в трехмерном случае уменьшение вязкости, при сохранении условий генерации турбулентности, сопровождается ростом среднеквадратичной завихренности, то в двумерном случае этого не происходит, т.к. dt 0. Т.о., при больших числах Рейнольдса кинетическая энергия в двумерном турбулентном течении остается практически постоянной. Каскадный процесс передачи энергии к высоким гармоникам, который затрагивает значительную долю от общей кинетической энергии потока, в двумерном случае исключен. С другой стороны, каскадный процесс переноса завихренности к большим волновым числам возможен. Чтобы показать это запишем уравнение относительно величины B rot :

Уравнение (2.28) фактически совпадает с трехмерным уравнением Гельмгольца (2.8). Если умножить (2.28) скалярно на B и проинтегрировать по объёму, то мы получим уравнение описывающее эволюцию палинстрофии:

Первый член в правой части уравнения (2.29) описывает скорость дизавихренности. Увеличение или уменьшение градиентов завихренности зависит от знака среднего значения квадратичной формы B, SB. По структуре уравнение (2.29), описывающее эволюцию палинстрофии в двумерном случае, ни чем не отличается от уравнения (2.11), описывающее эволюцию энстрофии в трехмерном случае. Следовательно и процессы растяжения вихрей в трехмерном случае и линий дизавихренности в двумерном случае имеют одну и ту же природу. Именно поэтому можно считать, что двумерная турбулентность имеет право на существование.

В теории KLB вводится понятия скорости диссипации среднеквадратичной завихренности B. Используя предположение, что const можно применить подход Колмогорова к изучению двумерной турбулентности. Из размерностных соображений будем иметь:

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

где C – универсальная константа.

Т.к. в двумерной турбулентности поток энергии в среднем направлен в сторону больших масштабов, а поток энстрофии в сторону малых масштабов, то и процессы переноса будут определяться двумя величинами – скоростью диссипации энергии и – скоростью диссипации энстрофии. Следовательно, если энергия и энстрофия вносятся в поток на некотором промежуточном масштабе k L k I k, то в двумерной турбулентности следует ожидать появления двух инерционных интервалов. В области k I k k формируется прямой энстрофийный каскад с законом подобия E k C 2/3k 3, а в области соответствии с теорией Колмогорова, должен обладать законом подобия Действительно, энстрофия и при kmax, тогда как закон сохранения (2.27) запрещает её рост. Чтобы избежать этой расходимости, к спектрам (2.32) надо ввести логарифмические поправки. Для этого положим [49] условие k const может обеспечить только «закон пяти третей», и при этом основной вклад в интеграл будут вносить волновые числа демонстрирует локальность спектрального переноса энергии. Для обеспечения условия k const необходимо положить При этом основной вклад в (2.34) будут вносить волновые числа k из интервала k k k0, так что спектральный перенос энстрофии оказывается нелокальным процессом.

Обзор результатов моделирования двумерной турбулентности Фрик [46] отмечает, что одним из наиболее удачных экспериментов по наблюдению двумерной турбулентности можно считать работу Соммериа года [52]. В этой работе исследовался обратный каскад энергии в плоском течении в тонком слое ртути, возбуждаемом электромагнитными силами на малых масштабах. Эксперимент проводился в вертикальном магнитном поле, достигавшем 1Тл. Такое сильное магнитное поле практически подавляет вертикальные движения и приводит к горизонтальному течению, с вертикальным профилем, описываемым известным профилем Гартмана f z. Полагая, что поле проинтегрировав уравнение Навье-Стокса по z, получим:

где коэффициенты и зависят от конкретного профиля течения в слое.

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

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

В эксперименте удалось получить турбулентные режимы, в которых наблюдалось формирование обратного каскада энергии со спектром «-5/3».

На рисунке 2.1 показан экспериментальный спектр пульсаций скорости.

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

Рис. 2.1. Энергитический спектр в двумерной турбулентности. Эксперимент Соммериа 1986 г.

Все численные эксперименты можно разделить на две группы. Первая группа – это эксперименты по свободному вырождению турбулентности, вторая – по стационарно возбуждаемой турбулентности. Свободное вырождение подразумевает отсутствие внешних сил. В этом случае решение зависит только от начальных условий. С точки зрения динамики инерционных интервалов более интересны эксперименты по моделированию стационарной турбулентности. Для обеспечения стационарности необходимо обеспечить подвод энергии. В двумерной турбулентности интересны динамические процессы по обе стороны от масштабов возбуждения, поэтому сила f записывается в пространстве Фурье таким образом, что она поддерживает на заданном уровне энергию гармоник с заданным модулем волнового числа k k I.

Что касается диссипативного слагаемого, то, во-первых, диссипация должна обеспечить отвод энергии на больших масштабах. Во-вторых, для получения модифицируют и характер трения и на малых масштабах. При написании обычного диссипативного слагаемого в Фурье-представлении получаем член вида D k ~ k 2. В численных экспериментах искусственным образом повышают степень волнового числа и записывают диссипацию в виде с типичным значением показателем степени n m 8. Диссипативный член вида (2.36) приводит к тому, что действие диссипации концентрируется в узких интервалах вблизи граничных условий.

В своем обзоре Боффета (2011) [53] отмечает, что первые численные расчеты обратного энергетического каскада были получены при DNSмоделировании Лилли (1969) [54], затем Фришем и Сулемом (1984) [55], Херингом и Мак Вильямсом (1985) [56]. Точность численных расчетов увеличивалась с ростом вычислительной мощности, подтверждая закон «-5/3» и уточняя значение безразмерной константы. Используя разрешение 20482, Смит и Яхот (1993) [57] получили значение C 7 и показали, что функция распределения флуктуаций скорости u r близка к распределению Гаусса, что говорит об отсутствии влияния перемежаемости на обратный каскад.

Рис. 2.2. a) Спектр энергии E ( k ) при DNS-моделировании обратного энергетического каскада с разрешением 20482. Форсинг действует на гармонике k 600, линейное трение демпфирует энергию около k 6. b) Функция распределения вероятности (PDF) пульсаций продольной Рис. 2.3. a) Поле завихренности после образования конденсата в численном эксперименте, b) эволюция энергетического спектра при отсутствии вязкости на больших масштабах, c) траектории маркеров в эксперименте (EML) d),e) энергетические спектры в EML Численные результаты были подтверждены экспериментально в серии электромагнитном слое (EML) Парет и Табелинг (1997) [58] получили зависимость энергетического спектра k 5/3 и нашли константу C 6.5 1.0. В 1998 году [59] эти же авторы провели эксперимент (EML) по изучению аналога «бозе-эйнштейновской конденсации» в двумерной турбулентности, т.е. эффекта, связанного с накоплением энергии на самой низкой моде в ограниченной области за счет действия обратного энергетического каскада при отсутствии диссипации энергии на больших масштабах.

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

Большой неожиданностью стал очень крутой спектр в инерционном интервале переноса энстрофии. Вместо закона «-3» численные эксперименты давали результаты от «-3.5» до «-5» и для свободно затухающей турбулентности (Herring 1974 [60], McWilliams 1984 [61]), и при наличии форсинга (Basdevant 1981 [62] и Legras 1988 [63]). Напомним, что в трехмерной турбулентности перемежаемость дает поправки к закону «-5/3» порядка нескольких сотых, а в двумерном случае расхождение составило единицы. В качестве примера рассмотрим результаты численных расчетов, взятых из работы Babiano 1995 [64] (рис. 2.4). На левом рисунке моделируется прямой энстрофийный каскад. Используется сетка 10242, случайная сила действует на волновое число k I 10, диссипативный член используется в форме (2.36). На правом рисунке делается попытка одновременно получить оба инерционных интервала. Используется сетка 17282, случайная сила действует на волновое число k I 40. Наклоны спектров в обоих случаях примерно равны «-3.5».

Рис. 2.4. Результаты численного эксперимента Babiano. a) k I 10, сетка 10242 b) k I 40, Расхождение результатов связано с существованием в турбулентном потоке долгоживущих когерентных структур – изолированных интенсивных вихрей, слабовзаимодействующих с окружающим их турбулентным потоком (рис. 2.4). В работе Babiano 1987 [65] был проведен интересный эксперимент. Изолируемые вихри разрушались искусственно таким образом, что при этом не изменялось распределение энергии по спектру (это делалось путем введения случайных фаз в Фурье-компоненты). В результате спектральное распределение возвращалось к виду (2.32).

Рис. 2.5. Когерентные структуры в двумерном турбулентном потоке. (слева) Поле концентрации пассивной примеси. (справа) Babiano 1987.

Уравнение переноса завихренности t u по форме совпадает с уравнением переноса пассивной примеси (в качестве примеси может выступать, например, температура) Однако аналогия между этими уравнениями не работает. На рисунке 2. показано поле концентрации пассивной примеси, полученное в том же эксперименте, что и поле завихренности на предыдущем рисунке. Существенное отличие заключается в том, что в поле пассивной примеси нет столь выраженных изолированных структур. След от каждого вихря можно ясно увидеть и в поле пассивной примеси, и это кажется естественным, но при этом не наблюдается интенсивный рост концентрации к центру вихря, как это имеет место в поле завихренности. На рисунке 2.6 приведены спектры пульсаций завихренности и пассивной примеси. Видно, что спектр пассивной примеси соответствует закону k 1, в то время как спектр завихренности (энстрофии) после сравнительно короткого участка с наклоном «-1», дает крутой спад с законом близким к «-3» («для энергетического спектра). Различие в результатах обусловлено тем, что если температура (пассивная примесь) действительно не зависит от скорости, то завихренность жестко связана со скоростью соотношением rot u.

Рис. 2.6. Спектры пульсаций завихренности и концентрации пассивной примеси. Babiano 1987.

Довольно интересна работа Ohkitani 1991 [66]. Он использовал т.н.

декомпозицию Вейса 1981 [67] для выделения гиперболических областей (h), в которых доминирующую роль играет тензор деформаций, и эллиптических областей (e), в которых главную роль имеет тензор вращений. Эллиптические области будут соответствовать когерентным структурам, а гиперболические промежуточной турбулентной пелене. Соответственно поле завихренности можно представить в виде отфильтровывав вклад когерентных структур.

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

Рис. 2.7. a) Спектры пульсаций завихренности от гиперболических областей. Сетки 128, 256, 512. Сплошная линия – наклон «-1», пунктирная – «-2» («-3» и «-4» соответственно для энергетического спектра); b) спектры от эллиптических областей, пунктирная – «-2», Ohkitani В спектре от гиперболических областей есть область с наклоном «-1», что соответствует классическому спектру энергии KLB «-3». Эллиптическая же часть приводит к укручению спектра (наклон «-2»).

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

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

продемонстрировано, что при разрушении всех когерентных вихрей с помощью сильной инфракрасной диссипации устанавливается спектр k 3. В более поздних работах [69], [70], [71] также были получены спектры типа KLB, и даже с учетом логарифмической поправки из-за нелокальности взаимодействия в Фурьепространстве. Достигалось это путем использования особого вида диссипации, уменьшающей влияние когерентных структур на спектры турбулентности. К примеру, в работе Ishihira и Kaneda 2001 [69] использовалась диссипация вида:

где k – завихренность (в фурье-пространстве), const 0 при k K, и 0 при k K, rms – rms завихренности, вычисляется на каждом шаге по времени, – настроечный параметр, K max зависит от мелкости сетки, n обычно равен 4. Форсинг организуется с помощью марковского процесса:

где A – параметр, отличный от нуля в интервале частот k K f min, K f max, – случайная фаза, R – функция от шага t и корреляционного времени форсинга.

На рисунке 2.8 приведены результаты расчетов энергетического спектра и потока энстрофии для двух разных наборов значений параметров входящих в (2.41)-(2.42) на сетках 10242 и 20482. Наклон в спектрах изменяется от «-3.14» до «-3.17».

Рис. 2.8. a) Энергетический спектр, b) поток энстрофии. Ishihira and Kaneda (2001) Также можно отметить работу Brachet et al. 1988 [72], где было продемонстрировано, что наклоны спектральных характеристик в свободно затухающей турбулентности могут зависеть от времени наблюдения. Численный эксперимент проводился на сетке с разрешением 20482, в начальный момент были турбулентности формировался спектр ~ k 4, и постепенно наклон уменьшался до значения ~ k 3, когда турбулентность становилась более «старой».

В своем обзоре двумерной турбулентности Табелинг 2002 [73] отмечает, что результаты моделирования прямого каскада энстрофии сильно зависят от начальных условий, размеров и распределения когерентных вихрей, времени наблюдения, структуры диссипации и форсинга в системе.

Что касается физических экспериментов, то первые попытки касались наблюдения энстрофийного каскада в тонком слое ртути в сильном магнитном поле (EML) (рис. 2.10). Спектр k 3 был получен, но при интерпретации диссипации, а также стенок. Позже эксперименты в мыльных пленках (рис. 4.12) показали, что энстрофийный каскад может развиваться в соответствии с классической теорией в реальных физических системах [74], [75], [76], [77], [78], [79]. Преимущество экспериментов с мыльными пленками по сравнению с EML заключается в том, что толщина мыльной пленки может быть в 100 раз меньше эффективного слоя электролита в магнитном поле. Это позволяет получать и исследовать турбулентность на гораздо более малых масштабах.

Рис. 2.10. Схема эксперимента со слоем электролита в магнитном поле. PIV метод. Tabeling P., Рис. 2.11. Схема эксперимента с мыльными пленками. LDV/PIV методы. Couder Y., Goldburg На рисунке 2.12 приводятся результаты экспериментов в мыльных пленках (взято из обзора Boffeta 2011 [53]) Рис. 2.12. a) E(k) в эксперименте в мыльных пленках (LDV-метод), Belmonte et al. 1999; b) поле завихренности в аналогичном эксперименте (PIV-метод); c) E(k) в эксперименте в мыльных В работе Belmonte et al. в 1999 году был проведен эксперимент в горизонтальной мыльной пленке с использованием специального «гребешка» для генерации турбулентности. Наклон спектра и инерционном интервале был близок к спектру KLB («-3.3»). В аналогичном эксперименте Rivera et al. 2003 (рис. 2. с) были получены спектры в зависимости от расстояния по потоку. Наклоны спектров близки к классическому («-3»).

Основной мотивацией для развития теории двумерной турбулентности, следуя Крейчнану, была надежда в то, что данная модель будет полезной для описания атмосферной турбулентности. Первые наблюдения за атмосферными спектрами показали их степенной характер. Анализ ветровых и температурных измерений, выполненных в программе GASP (Global Atmospheric Sampling Programm) Nastrom et al. [80, 81], подтвердил существование закона k 3 в области 1000-3000 км, а в области от 600 км и ниже k 5/3 (рис. 2.13).

Рис. 2.13.Спектральная плотность атмосферной турбулентности. Nastrom and Gage 1984.

В свете теории KLB двумерной турбулентности, наклон «-3» был интерпретирован как прямой энстрофийный каскад, обусловленный вводом энстрофии при бароклинной неустойчивости. Появление же спектра k 5/ длительное время оставался не ясным, т.к. было широко распространено мнение, что поток энергии к малым масштабам в двумерной турбулентности запрещен.

Можно объяснить это либо отражением обратного каскада при введении энергии на малых масштабах [82], либо как возникновение положительного потока энергии от больших структур к малым в результате распада длинных гравитационных волн. Однако последние наблюдения за атмосферой показали, что в ней могут существовать и прямой поток энстрофии, и прямой поток энергии на мезомасштабах от нескольких десятков километров, до нескольких тысяч километров [83]. Как отмечали Tung и Orlando [84], при некоторых условиях нисходящий поток энергии, хотя и небольшой по сравнению с восходящим, может обнаруживать себя в энергетическом спектре. Это и приводит к возникновению спектра k 5/3, который наблюдался в атмосфере в области высоких волновых чисел.

Спектры турбулентности порождаемые сингулярностями В связи с трудностями, возникшими при моделировании прямого каскада, появлялись также и альтернативные теоретические работы двумерной турбулентности. В 1971 году Саффман [85] предложил динамический сценарий формирования «филаментов» завихренности (линий с большими градиентами завихренности). Из его теории следовало, что наклоны энергетических спектров должны быть равны «-4», что хорошо подтверждалось в ряде численных работ, отмеченных в предыдущем разделе. Позже Моффат [47], подчеркивая важность спиралеобразных структур, предложил спектр «-11/3». Также следует отметить теорию конформного поля Полякова 1992 [86].

В своих работах [87-89] Кузнецов Е.А. отмечает: «Хорошо известно, что особенности, возникающие в результате нелинейной эволюции в непрерывных средах, продуцируют в коротковолновой области степенные хвосты в спектре турбулентности. Эта идея была использована Филлипсом [90] для нахождения спектра турбулентности волн на воде за счет образования барашков, т.е.

особенностей поверхности z x, y, t в виде заострений. В этом случае вторая производная вертикального отклонения содержит -функционные особенности, благодаря чему Фурье-образ от при больших волновых числах падает как k 2.

Отсюда, по Филлипсу, спектр E k оказывается пропорциональным k 3 …».

Применяя эти же рассуждения к двумерной турбулентности, где ожидаемые особенности могут быть связаны с резкими градиентами завихренности, можно прийти к распределению Т.е. с точностью до логарифмической поправки данный результат в точности совпадает со спектром в теории KLB. Однако при выводе этого соотношения неявно предполагалось, что особенности точечные, хотя в двумерной турбулентности они скорее оказываются распределенными вдоль линий. В частности, если предположить, что скачки завихренности распределены изотропно, то спектр турбулентности будет сильно отличаться от распределения (2.43). Как было показано Саффманом в 1971 году [85], спектр в этом случае имеет вид Появление резких градиентов завихренности, т.е. когда величина завихренности r испытывает скачки, толщина которых мала по сравнению с характерными масштабами турбулентности, наблюдалось в большом количестве численных экспериментов [61, 63, 91], моделирующих турбулентность при больших числах Рейнольдса – в режиме, когда в нулевом приближении можно вместо уравнений Навье-Стокса рассматривать уравнение Эйлера. Поэтому главным является вопрос: возможен ли такой процесс в рамках уравнения Эйлера? Кузнецов Е.А. с соавторами приводит следующие соображения в пользу следует, что B изменяется только благодаря нормальной компоненте скорости un. Тангенциальная компонента скорости u при этом играет пассивную роль, обеспечивая выполнение условия несжимаемости div u 0. Введем новые траектории, определяемые нормальной компонентой скорости:

Решение этой системы уравнений задает отображение r r a, t. В силу того, что траектории задаются не самой скоростью, а её нормальной компонентой, то Якобиан J этого отображения не фиксирован, он может принимать произвольные значения. Этот факт непосредственно следует из уравнения которое получается из (2.45) дифференцированием по переменным a. С помощью интегрирование В виду произвольных значений J, из уравнения (2.47) следует, что в поле предполагать, следуя Саффману, что в двумерной турбулентности при больших числах Рейнольдса происходит формирование скачков завихренности. Ширина скачков определяется из баланса инерционных и вязких сил. Нас будет интересовать область коротких длин, когда k лежит в интервале Для того чтобы найти спектр от разрывов, вначале запишем выражение для градиента от одного разрыва, предполагая его ориентированным вдоль y :

где G x непрерывная функция, отличная от нуля на интервале [ x1, x2 ] при G x1,2 0, а затем просуммируем по всем разрывам.

От одного разрыва Фурье-амплитуда k примет вид:

от ансамбля разрывов:

Чтобы найти спектр энстрофии необходимо произвести осреднение величины k по всем случайным переменным. Считая, что положения разрывов равномерно распределены в пространстве, получим где N – число разрывов в области S, угловые скобки означают осреднение и по x1, x2 и по угловому распределению.

Т.к. нас интересует коротковолновая часть спектра kL 1, то под интегралом в уравнении (2.52) стоит быстро осциллирующая функция от x.

Следовательно, интеграл можно оценить по методу стационарной фазы. Этот метод можно применить для всех углов k 0 ( k – угол между векторами k и n ), где 0. В области k 0 k 0 и плотность энергетического спектра имеет вид где n – плотность разрывов, Gl G x dx.

Для остальной области по методу стационарной фазы с учетом условия G x1,2 0 получим Эта формула имеет особенности при углах k близких к 0 и /2. При малых углах это выражение нужно заменить на (2.53). А при углах близких к / интеграл в (2.54) следует ограничить из за пересечения с линией разрыва.