«А.В. Гарбарук, М.Х. Стрелец, М.Л. Шур МОДЕЛИРОВАНИЕ ТУРБУЛЕНТНОСТИ В РАСЧЕТАХ СЛОЖНЫХ ТЕЧЕНИЙ Учебное пособие (рекомендовано Учебно-методическим объединением вузов РФ по образованию в области прикладных математики и ...»
Министерство образования и науки Российской Федерации
Санкт-Петербургский государственный политехнический университет
А.В. Гарбарук, М.Х. Стрелец, М.Л. Шур
МОДЕЛИРОВАНИЕ ТУРБУЛЕНТНОСТИ
В РАСЧЕТАХ СЛОЖНЫХ ТЕЧЕНИЙ
Учебное пособие
(рекомендовано Учебно-методическим объединением вузов РФ по
образованию в области прикладных математики и физики в качестве учебного пособия для студентов высших учебных заведений, обучающихся по направлению подготовки «Прикладные математика и физика») Санкт-Петербург Издательство Политехнического университета 2012 УДК 532+533 ББК Рецензенты:
Доктор физико-математических наук, профессор СПбГУ Н. В. Егоров Доктор технических наук, профессор СПбГПУ А. Ю. Снегирев Гарбарук А.В. Моделирование турбулентности в расчетах сложных течений:
учебное пособие / А.В. Гарбарук, М.Х. Стрелец, М.Л. Шур – СПб: Изд-во Политехн. унта, 2012. – 88 с.
Пособие содержит краткий очерк существующих подходов к моделированию турбулентности для аэродинамических расчетов, анализ их физического содержания и областей применимости. Приводится формулировка наиболее распространенных моделей:
осредненные по Рейнольдсу уравнения Навье-Стокса (RANS), метод моделирования крупных вихрей (LES), гибридные RANS-LES модели. Изложение сопровождается примерами применения этих моделей для расчета течений различных типов. Особое внимание уделяется гибридному подходу, получившему название метода моделирования отсоединенных вихрей (DES), который в настоящее время наиболее широко используется при решении прикладных задач аэрогидродинамики.
Предназначено для студентов высших учебных заведений, обучающихся по образовательным программам бакалаврской и магистерской подготовки направления «Прикладные математика и физика». Пособие может быть полезно и для студентов других направлений подготовки: «Механика и математическое моделирование», «Прикладная механика», «Авиастроение», «Энергетическое машиностроение», «Ядерная энергетика и теплофизика».
Печатается по решению редакционно-издательского совета Санкт-Петербургского государственного политехнического университета.
© Гарбарук А.В., Стрелец М.Х., Шур М.Л., © Санкт-Петербургский государственный ISBN политехнический университет,
ОГЛАВЛЕНИЕ
Стр.ВВЕДЕНИЕ
I. СОВРЕМЕННОЕ СОСТОЯНИЕ И ОБЛАСТИ ПРИМЕНИМОСТИ
РАЗЛИЧНЫХ ПОДХОДОВ К МОДЕЛИРОВАНИЮ
ТУРБУЛЕНТНЫХ ТЕЧЕНИЙ1.1. ОБЩАЯ ХАРАКТЕРИСТИКА СУЩЕСТВУЮЩИХ ПОДХОДОВ
1.2. ПОЛУЭМПИРИЧЕСКИЕ МОДЕЛИ ТУРБУЛЕНТНОСТИ
1.3. МЕТОД МОДЕЛИРОВАНИЯ ОТСОЕДИНЕННЫХ ВИХРЕЙ (DES) И ЕГО
МОДИФИКАЦИИ (DDES И IDDES)1.3.1. Общая характеристика DES и его возможностей
1.3.2. Недостатки оригинальной версии DES и модификации, предложенные для их устранения
1.3.2.1. Лавинообразное падение турбулентной вязкости при измельчении сетки или уменьшении числа Рейнольдса ниже порогового значения
1.3.2.2. Переключение DES в LES-режим внутри пограничного слоя, обусловленное измельчением шагов сетки в направлениях, касательных к обтекаемой поверхности.
1.3.2.3. Формирование двух логарифмических участков на профиле скорости при расчете присоединенных пограничных слоев
II. МАТЕМАТИЧЕСКАЯ ФОРМУЛИРОВКА НАИБОЛЕЕ
РАСПРОСТРАНЕННЫХ ПОДХОДОВ К МОДЕЛИРОВАНИЮ
ТУРБУЛЕНТНОСТИ2.1. УРАВНЕНИЯ РЕЙНОЛЬДСА И ПОЛУЭМПИРИЧЕСКИЕ МОДЕЛИ СПАЛАРТААЛЛМАРАСА И МЕНТЕРА
2.1.1. Уравнения Рейнольдса для сжимаемого газа
2.1.2 Модель Спаларта-Аллмараса (SA модель)
2.1.2.1. Оригинальная версия модели
2.1.2.2. Поправка к SA модели на кривизну линий тока и вращение (SARC модель)
2.1.2.3. Поправка к SA модели на сжимаемость потока (SACC)
2.1.3 Модель Ментера (SST модель)
2.2. МЕТОД LES
2.2.1. Основные уравнения
2.2.2. Подсеточная модель Смагоринского
2.3. МЕТОД DES
2.3.1. Общий принцип построения DES на основе различных RANS моделей
2.3.2. DES на основе SA RANS модели
2.3.2.1. Подсеточная версия SA модели
2.3.2.2. Объединение SA RANS и SA LES моделей.
2.3.3. Поправка к SA DES на низкие числа Рейнольдса
2.3.4. DES на основе SST модели
2.4. МЕТОД DDES
2.5. МЕТОД IDDES
2.5.1. Формулировка IDDES на основе SA и SST RANS моделей
2.5.2. Примеры применения IDDES
2.5.2.1. Установившееся течение в плоском канале
2.5.2.2. Течение в плоском канале с обратным уступом на нижней стенке
СПИСОК ЛИТЕРАТУРЫ
ВВЕДЕНИЕ
Поиск приемлемых для практики форм математического описания турбулентных течений или, как принято говорить, моделей турбулентности вот уже на протяжение более 100 лет (начиная с классических работ О. Рейнольдса) занимает умы многих выдающихся математиков и механиков. Это объясняется как исключительной сложностью турбулентности как физического явления, так и тем обстоятельством, что именно турбулентная форма движения газов и жидкостей наиболее часто реализуется в природе и в различных технических приложениях. Специфика ситуации состоит в том, что в отличие от многих других физических явлений, для описания которых строгие математические модели отсутствуют, турбулентные течения, согласно современным представлениям, подчиняются классическим уравнениям Навье-Стокса, и в этом смысле проблема может считаться давно решенной. Однако, несмотря на фантастический прогресс вычислительной техники, наблюдаемый в последние десятилетия, ее возможности все еще недостаточны для решения этих уравнений при представляющих практический интерес высоких числах Рейнольдса и, даже по самым оптимистичным прогнозам, будут оставаться таковыми по крайней мере вплоть до второй половины, а то и до конца XXI века. В связи с этим, как и ранее, ключевым вопросом в рассматриваемой области является поиск приемлемого компромисса между физической адекватностью модели и приемлемым для практического применения уровнем ее сложности. Естественно, что определение “приемлемости” существенным образом зависит от возможностей вычислительной техники, с одной стороны, и от потребностей промышленности, являющейся основным потребителем результатов расчетов турбулентных течений, с другой стороны. Поскольку и те, и другие достаточно быстро растут, представление о приемлемости компромисса быстро эволюционирует в сторону все более сложных моделей, что, в свою очередь, требует проведения новых исследований, направленных на определение границ их применимости. Таким образом, несмотря на свою более чем вековую историю, математическое моделирование турбулентных течений по-прежнему остается исключительно актуальной и быстро развивающейся областью классической и вычислительной гидро- и аэродинамики.Данное пособие ни в коей мере не претендует на замену многочисленных учебников, посвященных последовательному изложению этой области науки. Его цель гораздо скромнее и состоит лишь в том, чтобы дать студентам соответствующих специальностей общее представление о физической сущности, возможностях и ограничениях современных подходов к моделированию турбулентности, что, по мнению авторов, является не только необходимой предпосылкой, но и стимулом к осознанному более глубокому изучению различных аспектов проблемы моделирования турбулентности.
Пособие состоит из двух глав.
В первой главе дана общая характеристика существующих подходов к моделированию турбулентности. При этом акцент сделан на лежащих в их основе предпосылках и вытекающих из них физических и вычислительных ограничениях. Во второй главе приведены математические формулировки наиболее популярных сегодня моделей, некоторые из которых ранее не были представлены в учебной литературе на русском языке. В обеих главах изложение сопровождается многочисленными примерами расчетов с использованием тех или иных подходов и сопоставлением результатов этих расчетов с экспериментальными данными.
Завершая данное краткое введение, авторы считают приятным долгом выразить искреннюю признательность своим учителям – профессорам кафедры гидроаэродинамики СПбГПУ Льву Герасимовичу Лойцянскому, Юрию Викторовичу Лапину и Николаю Ивановичу Акатнову, заразившим их своим энтузиазмом и интересом к изучению турбулентности и тем самым в значительной степени определившим их дальнейшую научную деятельность. Авторы также благодарны профессору А.Ю. Снегиреву, который внимательно ознакомился с рукописью пособия и сделал ряд полезных замечаний.
I. Современное состояние и области применимости различных подходов к моделированию турбулентных течений 1.1. Общая характеристика существующих подходов Несмотря на бурный (экспоненциальный) рост производительности компьютеров и значительные успехи, достигнутые в последние годы в области построения эффективных численных алгоритмов для решения задач аэродинамики и теплообмена, расчет турбулентных течений, как и на протяжении многих предшествующих десятилетий, является одной из наиболее сложных проблем вычислительной аэродинамики. Более того, надежное предсказание характеристик турбулентных потоков все еще остается скорее исключением, чем правилом, что объясняется исключительной физической сложностью турбулентности, в частности ее стохастической природой, принципиально трехмерным нестационарным характером и широким спектром пространственно-временных масштабов.
Вместе с тем, общий прогресс вычислительной аэродинамики, разумеется, не мог не сказаться и на состоянии проблемы моделирования турбулентности. В частности, в последние годы все большее применение находят подходы к моделированию турбулентности, базирующиеся на первых принципах аэродинамики (метод прямого численного моделирования - в англоязычной литературе Direct Numerical Simulation или DNS и метод моделирования крупных вихрей - Large Eddy Simulation или LES). Однако, из-за крайней вычислительной трудоемкости этих подходов их широкое практическое использование при решения сложных задач аэродинамики может начаться лишь в конце заимствованная из работы [1], опубликованной в 2000 году. В таблице представлены оценки вычислительных ресурсов, необходимых для расчета обтекания типичного гражданского самолета или автомобиля с использованием всех известных методов расчета турбулентных течений, начиная от полуэмпирических методов, базирующихся на осредненных по Рейнольдсу уравнениях Навье-Стокса (Reynolds Averaged Navier-Stokes – Это не означает, конечно, что уже сегодня DNS и, в особенности, LES не могут применяться при расчете более простых течений.
RANS) и кончая полностью свободным от эмпиризма методом DNS. Данные этой таблицы для RANS основаны на реальном опыте использования соответствующих методов, имевшемся в 2000 г., а прогноз готовности методов DNS, LES и DES (Метод моделирования отсоединенных вихрей - Detached-Eddy Simulation) сделан на основе весьма оптимистичной оценки темпов роста производительности компьютеров (в два раза каждые пять лет). При этом оценки вычислительных ресурсов, необходимых для DNS, LES и DES основываются на общих представлениях о характеристиках турбулентности и свойствах указанных методов.
Вычислительные ресурсы и перспективы практического применения различных подходов к моделированию турбулентных течений [1] Под готовностью подразумевается возможность расчета одного варианта в течение суток на самых мощных из доступных компьютеров **) Имеется в виду LES c пристеночным RANS моделированием; в случае LES вплоть до твердых стенок, затраты оказываются сопоставимыми с затратами DNS ***) На компьютере с производительностью 1 терафлоп, время расчета составляет 5000 лет!
Так, единственное (общепринятое в настоящее время) допущение, на котором базируется DNS, состоит в том, что уравнения Навье-Стокса адекватно описывают не только ламинарные, но и турбулентные течения. Соответственно, в рамках этого подхода расчет турбулентных течений производится путем непосредственного (без какого-либо предварительного осреднения) численного решения уравнений Навье-Стокса. При этом независимо от характера осредненного течения (то есть, независимо от того является ли оно двумерным или трехмерным, стационарным или нестационарным) должны использоваться трехмерные нестационарные уравнения Навье-Стокса, поскольку турбулентность является принципиально трехмерным и нестационарным явлением. Кроме пространственно-временных масштабов турбулентности. Наглядное представление об этих масштабах дает рисунок 1.1, на котором изображен типичный энергетический спектр турбулентности (зависимость кинетической энергии турбулентности от волнового числа) при достаточно высоких числах Рейнольдса. Этот спектр имеет три области.
Область I соответствует крупномасштабным “энергонесущим” турбулентным вихрям с размерами порядка интегрального линейного масштаба рассматриваемого течения L (ему отвечает волновое число kI=2/L), черпающим энергию из осредненного течения.
Колмогоровского масштаба турбулентности) и волновыми числами k kd 2 /, вязкая диссипация которых переводит кинетическую энергию турбулентности в тепло.
Наконец, область II или инерционная область спектра, лежащая между областями I и Ш, соответствует вихревым структурам с размерами в диапазоне L < l<. Влияние вязкости в этой области никак не проявляется, и энергия турбулентности не генерируется и не диссипирует, а лишь передается от более крупных вихрей к менее крупным (так называемый энергетический каскад). В результате, энергетический спектр в области II оказывается универсальными и описывается законом Колмогорова E~k-5/3.
Рисунок 1.1. Различные области энергетического спектра турбулентности при высоких значениях числа Рейнольдса турбулентности L/ пропорционально числу Рейнольдса в степени 3/4, L/=Re3/4, в результате чего размер пространственной сетки, необходимой для проведения расчетов с помощью DNS, растет с увеличением числа Рейнольдса как Re9/4. Наряду с этим, с ростом числа Re увеличивается также и отношение интегрального I и минимального (соответствующего Колмогоровским вихрям) временных масштабов ( / )1 / 2, определяющее число шагов по времени, необходимое для проведения расчета:
I / ~ Re1 / 2. В итоге, суммарные затраты на проведение DNS растут с ростом числа Рейнольдса как Re11/4.
Именно эти оценки и определяют приведенные в Таблице 1.1 данные о числе узлов сетки и числе временных шагов, необходимых для DNS реальных течений.
В силу указанных обстоятельств, в настоящее время DNS может использоваться лишь для расчета течений с относительно низкими числами Рейнольдса и применяется, главным образом, в фундаментальных исследованиях, целью которых является получение детальной информации о структуре и основных закономерностях турбулентности. Это нисколько не умаляет важности данного подхода, поскольку полученные с помощью DNS результаты, наряду с экспериментальными данными, составляют основу для калибровки и тестирования полуэмпирических моделей турбулентности. Более того, следует иметь в виду, что в будущем DNS станет, по-видимому, доминирующим подходом не только в аэродинамике, но и в смежных областях техники, например, в химической технологии, атомной энергетике и т. д.
Вторым по трудоемкости из существующих подходов к моделированию турбулентности является метод моделирования крупных вихрей (LES). Этот подход сформировался в начале 80-х годов прошлого века (см., например, обзорную работу [2]).
Его идея состоит в замене “глобального” осреднения характеристик реального турбулентного течения по времени, на котором базируется вывод уравнений Рейнольдса (см. раздел 2.1), “фильтрацией” этих характеристик от коротковолновых неоднородностей или, иными словами, их пространственным осреднением по областям с размерами порядка размера фильтра (вывод уравнений LES кратко обсуждается в разделе 2.2).
Данную процедуру иллюстрируют рисунок 1.2. В его левой части показан результат фильтрации одномерного сигнала с помощью фильтра единичной ширины, а в правой некоторое реальное течение, сетка, используемая для его локального пространственного осреднения, и теряемая при этом информация о вихрях с размерами меньше размеров ячеек данной сетки.
Рисунок 1.2. Сравнение исходного и отфильтрованного сигналов Для вывода уравнений LES актуальные переменные f в уравнениях Навье-Стокса заменяются на сумму соответствующих “отфильтрованных” и “подсеточных” переменных f f f, а затем к полученным уравнениям применяется операция фильтрации. При этом величина f определяется выражением где G ( x x, ) - функция фильтра, x – координата рассматриваемой точки потока, а – ширина фильтра (два примера часто используемых фильтров представлены на рисунке 1.3).
Система уравнений LES по форме аналогична системе уравнений RANS (см. разделы 2.1.1 и 2.2.1). Однако физическое содержание этих двух систем различно. Так, дополнительные (содержащие напряжения Рейнольдса) члены RANS описывают влияние всех турбулентных неоднородностей на осредненное решение, в то время как аналогичные члены уравнений LES (“подсеточные” напряжения) описывают влияние только относительно мелких (с размерами меньшими размера фильтра ) вихрей на зависящее от времени решение отфильтрованных уравнений. Иными словами, в рамках LES вихревые структуры с размерами, превышающими размеры фильтра, разрешаются “точно”2, а моделируются лишь вихревые структуры меньших размеров. Для того, чтобы подчеркнуть это, модели турбулентности для LES называют “подсеточными”.
Из приведенного описания становится ясным, что если размеру фильтра соответствует волновое число k, лежащее в универсальной (“инерционной”) области энергетического спектра турбулентности, то есть, если k I k k d (см. рисунок 1.4), то моделированию подлежат только относительно универсальные (не зависящие от конкретной геометрии и граничных условий) вихри. В результате, роль подсеточной модели в LES сводится к обеспечению правильной скорости каскадной передачи энергии турбулентности от крупных к мелким вихрям в пределах инерционного интервала волновых чисел или, иными словами, - правильной скорости диссипации наименьших из “разрешенных” вихрей.
Рисунок 1.4. Масштабы турбулентных структур, разрешаемые и моделируемые в рамках LES Именно в этом состоит принципиальное преимущество LES перед RANS подходом, в энергосодержащих вихрей, не подчиняющихся каким-либо универсальным законам. С практической точки зрения это преимущество означает существенное снижение Под этим подразумевается только отсутствие ошибок, связанных с физическим моделированием (некоторая погрешность численного решения, разумеется, всегда имеет место).
требований к подсеточным моделям для LES по сравнению с моделями турбулентности для замыкания уравнений RANS. В частности, опыт применения LES убедительно свидетельствует о том, что при выполнении условия k I k k d, он обеспечивает высокую точность расчета не только осредненных и основных статистических, но и пульсационных характеристик турбулентности даже при использовании простейших подсеточных моделей, например, классической алгебраической модели Смагоринского (аналог модели пути смешения Прандтля для RANS – см. раздел 2.2.2). Это преимущество LES играет чрезвычайно важную роль при решении задач аэроакустики, аэроупругости и при расчете турбулентных химически реагирующих течений. Отметим также, что при расчете свободных турбулентных течений с использованием монотонных численных методов оказывается возможным вообще отказаться от моделирования подсеточной турбулентности, так как в этом случае функции подсеточной модели с успехом выполняет присущая таким методам численная диссипация (см., например, [3], [4]). Это направление в LES получило название Monotonically Integrated LES (MILES) или Implicit LES (ILES)3.
Естественной платой за описанные преимущества LES является то, что он требует несопоставимо больших вычислительных ресурсов, чем RANS. Это связано с необходимостью, как и в случае DNS, проведения трехмерных нестационарных расчетов на достаточно мелких сетках даже в тех случаях, когда целью расчета является определение только параметров осредненного течения. С другой стороны, по понятным причинам (мелкомасштабная часть спектра моделируется, а не рассчитывается “точно”) ресурсы, необходимые для реализации LES, оказываются намного меньшими, чем для DNS. Так, для расчета турбулентности вдали от твердых стенок число ячеек сетки, необходимой для проведения LES, увеличивается с ростом числа Рейнольдса намного медленнее, чем в случае DNS: пропорционально Re0.4, а не Re2.25. Однако вблизи стенок, где даже энергонесущие вихри имеют очень маленькие размеры, требования к сеткам для LES существенно ужесточаются и приближаются к аналогичным требованиям для DNS:
число ячеек, необходимых для LES таких течений пропорционально Re1.8. В качестве иллюстрации на рис. 1.5 показаны зависимости от числа Рейнольдса числа ячеек сетки, которые должны находиться во внутренней и внешней областях пограничного слоя на пластине при его расчете в рамках LES [5]. Из рисунка видно, в частности, что при Re= Данный термин часто используется также в другом смысле, а именно, применительно к наиболее распространенной на практике версии LES, в которой явная фильтрация исходных уравнений не производится, и роль фильтра играет вычислительная сетка.
полное число ячеек составляет ~108.5, причем 99% из них находятся во внутренней части пограничного слоя.
Рисунок 1.5. Зависимость числа ячеек сетки во внутренней и внешней областях пограничного слоя на плоской пластине и общего числа ячеек, необходимых для LES данного течения, от числа Рейнольдса Именно это обстоятельство, которое делает LES сложных пристеночных течений при представляющих практический интерес высоких числах Рейнольдса невозможным не только в настоящее время, но и в обозримом будущем (см. Таблицу 1.1), послужило стимулом для создания гибридных RANS-LES подходов, первым и наиболее развитым из которых является метод DES. Следует подчеркнуть, что первый расчет обтекания самолета с помощью этого метода, в полном соответствии с прогнозом, приведенным в Таблице 1.1, был выполнен вскоре после его опубликования в 2000 г. [6] (см. раздел 1.3.1, рис. 1.23), что свидетельствует также в пользу достаточной надежности оценок сроков “готовности” методов, приведенных в этой таблице.
Таким образом, в течение ближайших тридцати лет основным рабочим инструментом для решения прикладных задач аэродинамики будут оставаться полуэмпирические методы, базирующиеся на использовании RANS в сочетании с различными полуэмпирическими моделями турбулентности, и метод DES, который также в значительной степени опирается на эти модели. В связи с этим становимся на них более подробно.
Характеризуя современное состояние полуэмпирических моделей турбулентности, предназначенных для замыкания уравнений Рейнольдса, приходится констатировать, что надежды на создание универсальной модели такого типа (то есть, модели, пригодной для расчета всех или, по крайней мере, большинства аэродинамических течений), казавшиеся вполне реальными в 70-80-х годах прошлого века, в настоящее время представляются неосуществимыми. Эта пессимистичная оценка базируется на многолетнем опыте применения полуэмпирических моделей турбулентности и на отсутствии реального прогресса в области построения таких моделей в течение последних 20 лет [7].
Рисунок 1.6. Примеры когерентных вихревых структур: а – сверхкритическое обтекание цилиндра; b – турбулентный пограничный слой при наличии положительного продольного градиента давления; с – “след” острова в океане; d – извержение вулкана, e – галактические облака; f – Солнце Причина этого состоит в том, что наряду с “универсальной” мелкомасштабной турбулентностью, существенное влияние на параметры турбулентных потоков оказывают относительно устойчивые крупномасштабные (с размерами порядка макро-масштаба течения), принципиально трехмерные нестационарные турбулентные структуры (“вихри”), некоторые примеры которых показаны на рис. 1.6. Характеристики этих Для того, чтобы подчеркнуть устойчивость этих структур в литературе их называют когерентными структурами.
структур зависят от конкретной геометрии рассматриваемого течения и граничных условий и, следовательно, не могут быть описаны в рамках полуэмпирических моделей, не учитывающих этих важных обстоятельств5.
Следует отметить, что признаки определенного “кризиса” полуэмпирической теории турбулентности стали проявляться уже в 70-х - 80-х годах прошлого века, когда начался лавинообразный рост числа публикаций, посвященных модификации известных и разработке новых полуэмпирических моделей, что само по себе свидетельствует о неудовлетворенности многочисленными моделями такого типа, уже существовавшими в то время. В результате было предложено огромное число моделей, и весьма остро встал вопрос об обоснованном выборе той или иной из них при расчете конкретных течений. В связи с этим в США, а затем и в Западной Европе начались систематические исследования, направленные на тестирование и определение границ применимости полуэмпирических моделей турбулентности. Эта работа проводилась и продолжает проводиться как отдельными исследователями и исследовательскими группами, так и в рамках масштабных международных программ, координируемых NASA и Стэнфордским университетом в США, Комиссией ЕС по развитию научных исследований (6-ая и 7-ая рамочные программы) и Европейским Сообществом по Течениям, Турбулентности и Горению (ERCOFTAC). В частности, значительный вклад в решение данной проблемы внесли Стэндфордские международные конференции 1968, 1980 и 1990 гг., получившие неофициальное название “олимпиад моделей турбулентности”. Аналогичные цели преследуются регулярно проводимыми в Западной Европе и в США специальными рабочими семинарами. На этих семинарах результаты расчетов различных типов турбулентных течений, полученные с применением различных моделей турбулентности несколькими исследовательскими группами, сравниваются между собой и со специально отобранными (достаточно надежными и информативными) экспериментальными данными, а также (при наличии такой возможности) – с результатами расчетов в рамках LES и DNS. Наконец, большой объем исследований рассматриваемой проблемы выполнен в рамках Европейских проектов ECARP [8] и FLOMANIA [9].
Это в равной степени относится как к простейшим алгебраическим и дифференциальным моделям, базирующимся на понятии турбулентной вязкости (гипотеза Буссинеска о линейной связи между тензорами Рейнольдсовых напряжений и скоростей деформаций), так и к наиболее “строгим” моделям переноса Рейнольдсовых напряжений.
В результате этих целенаправленных усилий удалось накопить обширную и, что особенно важно, объективную (практически свободную от вычислительных неточностей) информацию о возможностях различных полуэмпирических моделей турбулентности при расчете тех или иных типов турбулентных течений. Анализ этой информации свидетельствует о том, что применительно к задачам внешней аэродинамики, связанным с расчетом безотрывных течений и течений с ограниченными отрывными зонами, наиболее высокий “рейтинг” имеют две модели турбулентности: модель Спаларта и Аллмараса (SA модель) [10], и модель Ментера ( k - Shear Stress Transport или SST модель) [11], [12].
Первая из них содержит только одно дифференциальное уравнение переноса (для модифицированной турбулентной вязкости - см. раздел 2.1.2) а вторая – два таких уравнения (для кинетической энергии турбулентности k и удельной скорости ее диссипации - см. раздел 2.1.3). К сожалению, отдать однозначное предпочтение одной из этих моделей не представляется возможным, вследствие чего выбор той или иной из них при расчете того или иного течения является в значительной степени “делом вкуса” пользователя. С определенной уверенностью можно констатировать лишь то, что SA модель, как правило, несколько “затягивает” отрыв пограничного слоя, индуцированный неблагоприятным продольным градиентом давления. Кроме того, обе модели занижают темп релаксации пограничного слоя, формирующегося вниз по потоку от точки присоединения, к своему “равновесному” состоянию6 и значительно завышают размеры так называемого “углового отрыва”, то есть, отрыва от поверхности двухгранного угла при наличии неблагоприятного градиента давления (например, угла между крылом и фюзеляжем самолета). Следует также отметить, что даже для некоторых присоединенных течений ни та, ни другая модель не позволяют получить результаты, удовлетворяющие современным требованиям к точности расчета аэродинамических характеристик. К таким течениям, в частности, относятся течения с сильной кривизной линий тока и вращением, для которых предпочтительными считаются модели Рейнольдсовых напряжений (см, например, [13]).
Приведем некоторые конкретные примеры, иллюстрирующие кратко описанные выше свойства SA и SST моделей.
Отметим, что этот недостаток в той или иной степени присущ всем известным полуэмпирическим моделям турбулентности.
Для расчета таких течений с использованием SA модели разработана специальная поправка (см.
раздел 2.1.2.2).
Рисунок 1.7. Сравнение результатов расчета коэффициента трения и профиля скорости в турбулентном пограничном слое на плоской пластине, полученных с использовании SA и SST моделей, с теорией и с экспериментальными корреляциями Так, рисунок 1.7, заимствованный из базы данных NASA по моделям турбулентности [14], свидетельствует о том, что при расчете турбулентного пограничного слоя на пластине обе модели обеспечивают практически одинаковую и достаточно высокую точность предсказания основных характеристик этого течения.
Рисунок 1.8. Сравнение результатов расчета подъемной силы аэродинамического профиля NACA0012 при угле атаки =15о, рассчитанной с помощью различных CFD кодов с использованием SA и SST моделей, с экспериментальными данными (см. [14]) Также близкие и достаточно точные результаты дают обе модели при расчете обтекания профиля NACA0012 даже при достаточно больших (предотрывных) углах атаки (см. рисунок 1.8) и при моделировании многих других аналогичных течений.
Однако при расчете отрывных течений ситуация заметно изменяется.
Например, как видно из рисунка 1.9, SA и SST модели предсказывают заметно отличающиеся друг от друга координаты точки отрыва на наклонной стенке ассиметричного диффузора (так называемый диффузор Obi), причем, согласно SA модели, точка отрыва находится ниже по потоку и ближе к экспериментальному значению [15], чем согласно SST модели. С другой стороны, SST модель существенно лучше, чем SA модель, описывает изменение трения вдоль верхней (прямой) стенки диффузора. Кроме того, данный рисунок, а также рисунок 1.10 служат иллюстрацией упомянутого выше факта завышения обеими моделями длины участка релаксации течения к равновесному состоянию вниз по потоку от точки присоединения8.
Рисунок 1.9. Геометрия диффузора [15] (a) и сравнение экспериментальных распределений трения вдоль наклонной и прямой стенок асимметричного плоского диффузора [15] с расчетными распределениями, полученными с использованием SA и SST моделей (b) Еще одним важным классом течений, при расчете которых SA и SST модели приводят к существенно отличающимся друг от друга результатам, является трансзвуковое обтекание аэродинамических профилей. В качестве примера на рисунке 1.11 сравниваются полученные с помощью обеих моделей продольные распределения коэффициентов давления и трения по поверхности профиля RAE Расчеты, результаты которых представлены на рис. 1.9 и ниже (если это не оговорено особо), выполнены авторами настоящего пособия и их коллегами с помощью кода NTS (Numerical Turbulence Simulation) [16, 17].
(Режим 10 – число Маха M=0.75, угол атаки =3.19о [18]). Характерной особенностью данного течения является отрыв потока за скачком уплотнения, замыкающим сверхзвуковую зону, формирующуюся на верхней поверхности профиля. Из рисунка видно, в частности, что согласно SA модели, скачок находится несколько ниже по потоку и является более интенсивным, чем согласно STT модели. В результате, при использовании SST модели оторвавшийся пограничный слой вновь присоединяется к поверхности профиля (трение на поверхности профиля вновь становится положительным – см. рис.1.11), а при использовании SA модели отрывной пузырь простирается вплоть до задней кромки профиля. При этом результаты расчетов по SST модели несколько лучше согласуются с экспериментом.
Рисунок 1.10. Сравнение экспериментальных профилей продольной скорости (верхний ряд) и касательных напряжений (нижний ряд) в различных сечениях диффузора с результатами расчетов с использованием SA и SST моделей Рисунок 1.11. Сравнение расчетных распределений трения и давления по поверхности профиля RAEA2822 при трансзвуковом режиме течения (Режим 10), полученных с использованием SA и SST моделей, с результатами измерений [18] На рисунках 1.12, 1.13 показаны результаты расчетов обтекания обратного уступа, исследовавшегося в экспериментах [19]. Видно, что хотя как SA, так и SST модель качественно правильно описывают данное течение, ни одна из них не обеспечивает достаточно высокой точности предсказания его основных характеристик. В частности, обе модели завышают интенсивность возвратного течения в отрывной зоне и, как отмечалось выше, предсказывают слишком медленное восстановление параметров пограничного слоя вниз по потоку от точки присоединения. Тем не менее, в данном случае, как и при расчете трансзвукового обтекания профиля, SST модель, в целом, заметно превосходит по точности SA модель.
Рисунок 1.12. Схема канала с обратным уступом, расчетная сетка и сравнение расчетных распределений трения вдоль стенки канала за уступом, полученных с использованием SA и SST моделей, с экспериментальными данными [19] Рисунок 1.13. Сравнение расчетных профилей скорости, полученных с использованием SA и SST моделей, с экспериментальными профилями [19] в различных сечениях канала с обратным уступом на нижней стенке На рисунках 1.14-1.18 приведено еще несколько примеров турбулентных течений, для которых ни одна из двух рассматриваемых моделей не обеспечивает удовлетворительной точности расчета.
Рисунок 1.14. Геометрия крыла конечного размаха NACA0012 со скругленной кромкой (a), поле модуля завихренности (b) и сравнение расчетных и экспериментальных профилей поперечной и продольной составляющих скорости в различных сечениях следа за крылом (c) Первое из них представляет собой вихревой след за крылом конечного размаха с профилем NACA 0012 и скругленной боковой кромкой (рис.1.14а). Характерной особенностью данного и многих других аналогичных течений является сильная кривизна линий тока. Как отмечалось выше, расчет такого рода течений с помощью обычных полуэмпирических моделей турбулентности, базирующихся на гипотезе Буссинеска (см.
раздел 2.1.1), приводит к большим погрешностям. Это хорошо видно из рисунка 1.14:
согласование с экспериментом поля скорости в следе, рассчитанного с помощью SA и SST моделей, резко ухудшается вниз по потоку, причем обе модели сильно завышают скорость диссипации кромочного вихря. Для устранения этого дефекта необходимо введение в модели специальных поправок9.
Второй из упомянутых примеров – обтекание модели автомобиля (Ahmed car body [20] – см. рисунок 1.15а) – является известным международным тестом для оценки возможностей моделей турбулентности при расчете отрывных течений (см., например, [9]).
Рисунок 1.15. Сравнение результатов расчетов обтекания упрощенной модели автомобиля, полученных с использованием SA и SST моделей, с экспериментом [20] На рисунке 1.15 представлены некоторые результаты расчетов этого течения, полученные с использованием SA и SST моделей для случая, когда угол наклона заднего Такая поправка к SA модели представлена в разделе 2.1.2.2.
скоса “крыши” составляет 25о. Из этого рисунка видно, что SA модель предсказывает безотрывное течение в плоскости симметрии скоса, в то время как согласно SST модели в этой плоскости имеет место полностью отрывное течение. При этом ни та, ни другая модель не воспроизводит измеренное поле скорости, в котором в указанной плоскости имеет место отрыв потока, сопровождающийся его присоединением к поверхности “скоса” и последующим отрывом от его задней кромки.
Последний пример “проблемных” течений, расчет которых с использованием SA и SST моделей приводит к большим погрешностям, представлен на рисунке 1.16. Это течение внутри двугранного угла, образованного поверхностью аэродинамического профиля ONERA (A-airfoil) и плоскими боковыми стенками секции аэродинамической трубы, в которой проводились эксперименты [21]. На рисунке представлено сравнение расчетных (полученных с использованием обеих моделей) и экспериментальной картин обтекания профиля при предотрывном угле атаки 13о. Из него видно, что обе модели предсказывают формирование в двугранных углах обширных отрывных зон, в то время как в эксперименте угловой отрыв не наблюдается.
Рисунок 1.16. Сравнение экспериментальных линий тока на поверхности аэродинамического профиля ONERA (A-airfoil), полученных методом масляной пленки при угле атаки атаки 13о (левая рамка), с расчетными линиями тока (направление потока – сверху вниз) Приведенные примеры достаточно убедительно свидетельствуют о том, что при расчете течений с ограниченными отрывными зонами SA и SST модели, как правило, приводят к существенно отличающимся друг от друга результатам, а в некоторых случаях обе модели оказываются неспособными с приемлемой точностью предсказать характеристики таких течений.
Что касается течений с протяженными отрывными зонами, то в настоящее время преобладающей является точка зрения, согласно которой RANS модели в принципе неспособны обеспечить приемлемую для практики точность расчета таких течений, и в этом смысле SA и SST модели не являются исключениями10.
Типичный пример, иллюстрирующий неудовлетворительную работу этих моделей даже при расчете относительно “простого” течения такого типа (обтекание двумерного аэродинамического профиля NACA0012 при больших углах атаки), представлен на рисунке 1.1811.
Рисунок 1.17. Сравнение расчетных зависимостей коэффициента сопротивления профиля NACA0012 от угла атаки при Re=105, полученных с использованием SA и SST RANS моделей, с экспериментальными данными [22] (наличие двух экспериментальных кривых при углах атаки около 15о связано с явлением гистерезиса) Все это привело в конце прошлого века к смещению акцентов в исследованиях, посвященных моделированию отрывных турбулентных течений, в сторону поиска отличных от RANS подходов к решению этой сложной задачи. В результате, в 1997г. был предложен метод моделирования отсоединенных вихрей (DES) [23], положивший начало “эры” гибридных RANS-LES методов моделирования пристеночной турбулентности, рассмотрению которых посвящен следующий раздел данного пособия.
Это связано с важной ролью крупномасштабных когерентных вихревых структур в таких течениях (см. рис. 1.6).
Отметим, что в данном случае стационарное решение уравнений Рейнольдса получить не удается, и приведенные на рисунке результаты получены путем осреднения по времени нестационарного периодического решения, полученного при численном интегрировании нестационарных уравнений Рейнольдса (Unsteady RANS – URANS).
1.3. Метод моделирования отсоединенных вихрей (DES) В данном разделе проанализированы предпосылки, послужившие стимулом для создания DES, объяснена основная идея этого метода и приведены примеры, иллюстрирующие как открываемые им новые возможности, так и недостатки, выявленные в процессе широкого практического использования DES на протяжение уже более чем десяти лет. Наряду с этим, в общих чертах описаны более поздние версии этого метода (DDES и IDDES), позволяющие в той или иной степени устранить эти недостатки.
Детальные математические формулировки DES и его модификаций представлены в Главе II пособия.
1.3.1. Общая характеристика DES и его возможностей Как уже отмечалось, метод DES был предложен в качестве альтернативы RANS и LES методам при расчете пристеночных течений с обширными отрывными зонами, для которых RANS модели не способны обеспечить приемлемую точность, а LES требует чрезмерно больших вычислительных ресурсов, что отодвигает возможность его широкого практического использования до середины нынешнего века. При этом учитывалось то обстоятельство, что львиная доля вычислительных затрат LES связана с расчетом пристеночной части присоединенных пограничных слоев, населенных энергонесущими вихрями малых размеров (см. раздел 1.1 и рисунок 1.5), и что расчет именно таких течений с помощью RANS является достаточно надежным и экономичным (см. предыдущий раздел). Это естественным образом привело авторов DES к мысли о создании комбинированной RANS-LES модели, которая сочетала бы в себе лучшие качества обоих методов, а именно, надежность и вычислительную эффективность RANS в присоединенных пограничных слоях с понятной физикой, высокой точностью и приемлемыми вычислительными затратами LES вдали от стенок.
Именно такая комбинация RANS и LES, предложенная в 1997г., получила название метода моделирования отсоединенных вихрей (Detached-Eddy Simulation). Это название подчеркивает принципиальное отличие DES от LES. Оно состоит в том, что в рамках DES “точно” рассчитываются не все энергонесущие вихри, а лишь “отсоединенные” вихри которые населяют отрывные зоны, а вихри, населяющие области присоединенных пограничных слоев, описываются обычными полуэмпирическими RANS моделями. Еще одна важная особенность DES состоит в том, что в рамках этого подхода в RANS и LES областях используется одна и та же “базовая” модель турбулентности, которая функционирует как RANS модель внутри пристеночного пограничного слоя и как ее подсеточный (Sub - Grid Scale или SGS) аналог вдали от твердых стенок. В первой версии DES [23] в качестве такой базовой модели использовалась SA модель, однако в дальнейшем были предложены модификации DES, базирующиеся на других RANS моделях, в частности, на SST модели [24].
Рисунок 1.18. Изоповерхность параметра 2 (модуль второго собственного числа тензора градиента скорости) из расчетов докритического обтекания кругового цилиндра, выполненных в рамках различных подходов к моделированию турбулентности: a) – двумерные стационарные уравнения Рейнольдса, SA модель; b) - двумерные нестационарные уравнения Рейнольдса, SA модель; с) - трехмерные нестационарные уравнения Рейнольдса, SA модель; d) - SA DES на грубой сетке; e) - SA DES на мелкой сетке Рисунок 1.18 дает общее представление о возможностях, открываемых DES по сравнению с другими подходами к моделированию отрывных течений на примере обтекания кругового цилиндра при числе Рейнольдса Re 1.4 10 5, соответствующем докритическому режиму обтекания (отрыв ламинарного пограничного слоя и переход к турбулентности в оторвавшемся слое смешения). Из него видно, в частности, что при расчете данного течения в рамках стационарных двумерных уравнений Рейнольдса (2D RANS), поле скорости оказывается симметричным относительно геометрической плоскости симметрии цилиндра.12 В случае решения той же задачи в рамках нестационарных двумерных уравнений Рейнольдса (2D URANS) в следе за цилиндром формируется двумерная вихревая дорожка Кармана, а при использовании трехмерных уравнений (3D RANS) на нее накладываются продольные вихревые структуры. Однако, при всех отмеченных принципиальных отличиях, ни один из этих трех подходов не позволяют с приемлемой степенью точности предсказать такие важные параметры осредненного течения как сопротивление цилиндра и длина зоны рециркуляции в его следе. В противоположность этому, при использовании DES все эти параметры удается рассчитать с достаточно высокой точностью даже на относительно грубых сетках (в рассматриваемом случае - порядка полумиллиона ячеек), причем при измельчении сетки погрешность расчета уменьшается и в пределе стремится к нулю (DES переходит в DNS).
Об этом убедительно свидетельствуют рисунки 1.19 и 1.20, на которых представлено сравнение результатов расчетов данного течения с использованием DES на разных сетках с соответствующими данными 2D URANS, которые являются сошедшимися по сетке (практически не изменяется при ее дальнейшем измельчении), и с результатами измерений [25], [26].
Рисунок 1.19. Сравнение распределений коэффициента давления по поверхности цилиндра, полученных путем осреднения по времени SA DES и SA URANS решений, с экспериментальными данными [25] (в “легенде” указаны также соответствующие значения коэффициента сопротивления) Это решение, вообще говоря, неустойчиво, но может быть получено при проведении расчета в половине расчетной области с использованием условий симметрии.
Рисунок 1.20. Сравнение результатов DES c измерениями c помощью PIV [26]: а – линии тока осредненного течения; b – осредненное по фазе мгновенное поле поперечной составляющей завихренности Благодаря исключительной простоте реализации DES и открываемой им принципиальной возможности достаточно точного расчета широкого круга сложных отрывных турбулентных течений, этот метод быстро получил широкое распространение и, в частности, был внедрен во всех основных коммерческих CFD кодах. Это, а также большое число исследований, выполненных с помощью DES как отдельными научными группами во всем мире, так и рамках трех специальных Европейских проектов (FLOMANIA [9], DESider [27] и ATAAC [28]), позволило накопить большой объем информации обо всех аспектах данного метода, включая его достоинства, недостатки и вычислительные особенности.
Остановимся более подробно на возможностях и недостатках исторически первых версий DES [23, 24], которые базируются на SA и SST RANS моделях турбулентности.
Как уже отмечалось, первоначально DES предназначался для расчета аэродинамических течений с обширным отрывом при высоких числах Рейнольдса, типичными примерами которых являются течения у плохообтекаемых тел (цилиндр, сфера, крыло под большим углом атаки и т. п.). В настоящее время в литературе имеется множество примеров, демонстрирующих высокую точность DES при расчете таких течений.
Так, на рисунке 1.21 показаны результаты одного из первых исследований, проведенных с помощью DES, а именно, результаты расчета обтекания профиля NACA0012. Для наглядности на этом рисунке представлены также уже обсуждавшиеся выше результаты расчетов этого течения с помощью SA и SST RANS моделей.
Рисунок 1.21. Мгновенное поле модуля завихренности из расчета обтекания профиля NACA под углом атаки 45о методом DES (а) и сравнение с экспериментом [22] расчетных зависимостей коэффициента сопротивления этого профиля от угла атаки при Re=105, полученных с использованием SA и SST RANS моделей и с помощью основанных на этих моделях версий DES (b) Рисунок убедительно демонстрирует существенное превосходство DES над RANS. В частности, из него видно, что при относительно низких значениях угла атаки ( 30 o ), RANS и DES дают близкие результаты, достаточно хорошо согласующиеся с экспериментом. Однако при больших углах атаки, при которых за профилем формируется обширная отрывная зона, результаты расчетов по RANS моделям начинают сильно отличаться между собой и далеки от эксперимента, в то время как обе версии DES предсказывают очень близкие и прекрасно согласующиеся с экспериментом значения сопротивления. Следует подчеркнуть, что хотя расчетная сетка в трехмерных DES расчетах была примерно в 30 раз больше двумерной сетки, использовавшейся в URANS, ее общее число узлов составляло около 400,000, что по современным меркам является весьма скромным (проведение описанных расчетов даже на персональных компьютерах занимает всего несколько десятков часов).
Один из ярких примеров успешного применения DES для расчета сложных внутренних течений и теплообмена представлен на рисунке 1.22. Речь идет о развитом турбулентном течении во вращающемся канале квадратного сечения с оребренными нагретыми стенками [29] (см. левую верхнюю рамку на рисунке 1.22). Число Рейнольдса, построенное по среднерасходной скорости и стороне канала, составляло 20000, а число Россби Ro L / U b, ( U b – среднерасходная скорость потока, L – сторона канала, а – угловая частота его вращения) изменялось в диапазоне от 0.18 до 0.67. Как видно из рисунка, DES позволяет рассчитать данное течение и теплообмен с точностью, не уступающей точности LES, на сетке, содержащей примерно в 3.4 раза (1.53) меньше узлов.
При этом, благодаря возможности увеличения шага интегрирования по времени при проведении расчетов на более грубой сетке, общие затраты машинного времени оказываются примерно в 5 раз меньше, чем в LES.
Рисунок 1.22. Расчетная область, поверхностная сетка и результаты DES развитого течения во вращающемся канале с оребренными стенками при различных скоростях вращения (верхний ряд) и результаты расчета теплообмена при отсутствии вращения [29] (Nu0 – число Нуссельта при развитом течении в канале без оребрения [30]) Приведем, наконец, еще несколько примеров применения DES, которые весьма убедительно свидетельствуют о том, что уже в настоящее время этот подход может с успехом применяется для решения исключительно сложных задач аэродинамики, представляющих непосредственный практический интерес, а также для решения задач аэроакустики.
Рисунок 1.23. Мгновенная изоповерхность Q -критерия ( Q 2 0.5 2 S2 ) из SA DES обтекания истребителя F-15 под углом атаки 65° (цветом показано давление) и сопоставление результатов расчетов коэффициентов интегральных сил и момента, рассчитанных с помощью SA RANS и SA DES [6] На рисунке 1.23 представлены результаты расчета обтекания истребителя F-15 под углом атаки 65о, полученные с помощью SA DES в работе [6] с использованием коммерческого CFD кода COBALT на трех последовательно измельченных сетках. Видно, что на самой мелкой из них (она содержит около 6 миллионов ячеек) погрешность расчета интегральных характеристик течения (подъемной силы, силы сопротивления и крутящего момента) оказывается меньше 6%, в то время как погрешность расчета этих характеристик в рамках SA RANS на той же сетке превышает 10%.
рисунках 1.24, 1.25.
Первый из них относится к расчету обтекания автомобиля Форд [31] (расчет выполнен с помощью коммерческого CFD кода STAR-CD) и показывает, что основными источниками генерируемого при этом шума являются боковое зеркало заднего вида, углубления для “дворников”, передняя стойка кузова, арки колес, покрышки и след за автомобилем.
Рисунок 1.24. Мгновенная изоповерхность интенсивности источников шума на поверхности автомобиля Форд из расчета с помощью SA DES [31] Второй пример (расчет обтекания конфигурации цилиндр – аэродинамический профиль [32]), дает количественное представление о возможностях DES при расчете шума, генерируемого в дальнем поле сложным турбулентным потоком, формирующимся при взаимодействии следа за тонким цилиндром с аэродинамическим профилем.
Рисунок 1.25. Изоповерхность завихренности из расчета обтекания профиля в следе за цилиндром методом DES и сравнение расчетного и экспериментального энергетических спектров шума, генерируемого системой в дальнем поле [32] 1.3.2. Недостатки оригинальной версии DES и модификации, предложенные для их устранения Как уже отмечалось, наряду с широким возможностями и важными достоинствами DES по сравнению с RANS при расчете течений с обширным отрывом, в процессе интенсивной эксплуатации DES были выявлены также некоторые недостатки данного метода. Остановимся на них несколько подробнее.
1.3.2.1. Лавинообразное падение турбулентной вязкости при измельчении сетки или уменьшении числа Рейнольдса ниже порогового значения Данный недостаток проявляется только при использовании тех версий DES, которые базируются на RANS моделях, содержащих специальные низкорейнольдсовые члены, например, при использовании SA-DES (см. раздел 2.3.2)13, и состоит в следующем.
При уменьшении числа Рейнольдса и/или при сильном измельчении расчетной сетки уровень вихревой (“подсеточной”) вязкости в LES-области DES уменьшается, что полностью соответствует физике LES. Однако, при достижении ею некоторой пороговой величины, примерно в 10 раз превышающей значений молекулярной вязкости, это Отметим, что SST RANS модель таких членов не содержит падение воспринимается SA моделью как “сигнал” о близости твердой стенки, что автоматически влечет за собой активацию низкорейнольдсовых членов модели. Это, в свою очередь, приводит к дальнейшему быстрому падению вихревой вязкости практически до нуля, и расчет превращается в “DNS на грубой сетке”, что естественно, ухудшает получаемое решение.
Рассмотрим два примера, иллюстрирующих данный недостаток SA-DES.
Рисунок 1.26. Мгновенные поля подсеточной вязкости, полученные при расчете DIHT на сетке 643 с помощью подсеточной версии SA модели и модели Смагоринского (см. разделы 2.2.2 и 2.3.2) и сравнение соответствующих расчетных спектров кинетической энергии турбулентности с экспериментом [33] и с теоретическим законом “-5/3” Первый из них – задача о расчете процесса затухания однородной изотропной турбулентности (Decay of Homogeneous Isotropic Turbulence или DIHT) при сравнительно низком числе Рейнольдса, которая является классическим тестом для оценки и калибровки подсеточных моделей турбулентности. Как видно из рисунка 1.26, на котором представлены результаты решения данной задачи с помощью подробно рассмотренных ниже (см. разделы 2.2.2 и 2.3.2.1) подсеточной версии SA модели и модели Смагоринского на сетке, содержащей 643 узлов. В первом случае подсеточная вязкость действительно оказывается очень низкой ( t / 1 ), в то время как во втором случае она достигает значения 5. В результате подсеточная версия SA модели предсказывает совершенно неверный спектр турбулентности (синяя кривая на рис. 1.26), а спектр, рассчитанный по модели Смагоринского (красная кривая), хорошо согласуется с теорией и экспериментальными данными. Отметим, что на сетке 323, вследствие значительного (примерно в 4 раза) увеличения уровня подсеточной вязкости, описанный недостаток подсеточной версии SA модели практически не проявляется: в пределах диапазона волновых чисел, разрешаемых на этой сетке, полученный с ее помощью спектр кинетической энергии турбулентности хорошо согласуется с результатами расчетов по модели Смагоринского и с экспериментом.
Второй пример, иллюстрирующий проявление данного недостатка DES, заимствован из работы [34], в которой с помощью версии DES, базирующейся на низкорейнольдсовой k- модели Уилкокса [35] выполнен расчет обтекания тандема кубов, расположенных на плоской стенке (см. рисунок 1.27).
Рисунок 1.27. Визуализация результатов расчета обтекания тандема кубов на плоской стенке с помощью k- DES [34] Рисунок 1.28. Поля вихревой вязкости и скорости в плоскости симметрии тандема кубических тел на плоской поверхности, рассчитанные с помощью SA DES на грубой (а) и мелкой (b) сетках [34], и экспериментальное поле скорости [36] (c) На рисунке 1.28 показаны осредненные по времени поля вихревой вязкости и скорости в плоскости симметрии потока, полученные в этом расчете при использовании грубой (cg) и мелкой (fg) сеток, и результаты измерений поля скорости в этом течении из работы [36].
Из рисунка видно, что при измельчении сетки максимальный уровень подсеточной вязкости падает в 20 раз (от 20 до 1), и что при этом происходит значительное ухудшение согласования результатов расчетов с экспериментом. В частности, расчет на мелкой сетке приводит к занижению интенсивности возвратного течения в зазоре между кубами и не воспроизводит наблюдаемого в эксперименте быстрого (уже при x=8) восстановления профилей продольной скорости вниз по потоку от второго куба.
Рисунок 1.29. То же, что и на рисунке 1.28, для горизонтального сечения потока при y=0.375h Кроме того, при измельчении сетки наблюдается также заметное ухудшение согласования расчетного и экспериментального полей скорости в горизонтальных плоскостях течения. В частности, как видно из рисунка 1.29, в решении на мелкой сетке появляются ложные вихри в окрестности передних угловых точек второго куба.
Очевидно, что ухудшение качества решения при измельчении сетки, наглядно продемонстрированное двумя рассмотренными примерами, является крайне нежелательным явлением в рамках любого подхода к моделированию турбулентности. В связи с этим в работе [37] был выполнен детальный анализ причин данного дефекта DES и предложен простой путь его устранения (см. раздел 2.3.3 данного пособия).
1.3.2.2. Переключение DES в LES-режим внутри пограничного слоя, обусловленное измельчением шагов сетки в направлениях, касательных к обтекаемой поверхности.
Как и дефект, рассмотренный выше, данный недостаток DES проявляется при измельчении расчетной сетки, однако он имеет совсем другую природу.
Как уже отмечалось, для правильной работы DES необходимо, чтобы весь пограничный слой находился в RANS области DES. Однако, оригинальная формулировка DES (см. раздел 2.3) гарантирует выполнение этого требования лишь при условии, что шаги сетки || в направлениях параллельных стенке удовлетворяют неравенству || ( -толщина пограничного слоя – см. верхнюю рамку на рисунке 1.30).
Рисунок 1.30. Три возможных типа продольных сеток в пограничном слое: a – типичная сетка для DES, b – промежуточная сетка, c – типичная сетка для LES Такие сетки являются типичными для RANS расчетов обтекания тел и, как правило, позволяют получить сошедшиеся (практически не изменяющиеся при дальнейшем измельчении сетки) решения. Однако в некоторых случаях, например, вследствие локальных особенностей формы обтекаемой поверхности для получения точного решения RANS требуется уменьшение ||, которое приводит к нарушению неравенства || и возникновению “промежуточной сетки”, показанной в нижней рамке рисунка 1.3014.
представленной в следующей главе данного пособия (см. раздел 2.3), использование сеток промежуточного типа приводит к тому, что переключение DES из RANS в LES режим происходит внутри пограничного слоя. Однако при этом величина || недостаточно мала для LES пограничного слоя (для этого требуется “LES сетка” с существенно меньшими значениями ||, показанная в нижней правой рамке на рисунке 1.30). В результате, модельная часть турбулентных напряжений занижается, а их разрешенная часть оказывается недостаточной для компенсации этого эффекта, что приводит к существенному снижению точности расчета, а при наличии положительного градиента давления может привести даже к ложному (“индуцированному сеткой”) отрыву пограничного слоя.
Это явление иллюстрирует рисунок 1.31, на котором сравниваются решения, полученные с помощью SA RANS и SA DES при расчете обтекания профиля и течения в сверхзвуковом сопле, полученные при использовании сеток промежуточного типа. Видно, что в рамках DES отрыв потока в обоих случаях происходит существенно раньше, чем при использовании RANS15, что противоречит исходной концепции DES, согласно которой весь присоединенный пограничный слой, включая точку отрыва, должен описываться RANS ветвью DES.
Для устранения индуцированного сеткой отрыва (Grid Induced Separation) в работе [37] была предложена соответствующая модификация исходной версии DES, получившая название Delayed DES (Задержанный DES). Этот метод и примеры, демонстрирующие его эффективность, представлены в разделе 2.4.1 второй главы данного пособия.
Такая же ситуация может возникнуть, например, из-за наличия сильных градиентов параметров потока в направлениях параллельных стенке, как это имеет место при взаимодействии с пограничным слоем скачков уплотнения, или при фиксированном шаге по продольной координате - из-за увеличения толщины пограничного слоя под воздействием положительного градиента давления или при уменьшении числа Рейнольдса.
В англоязычной литературе это явление получило название Grid-Induced Separation или GIS.
Рисунок 1.31. Ложный отрыв пограничного слоя от поверхности профиля (а) и от стенок сопла (b) при проведении SA-DES на сетках промежуточного типа [38] 1.3.2.3. Формирование двух логарифмических участков на профиле скорости при расчете присоединенных пограничных слоев Как уже неоднократно подчеркивалось, представленные выше оригинальный метод DES и его модификации первоначально предназначались только для расчета отрывных течений, и применение DES для LES с пристеночным моделированием (Wall Modeling LES или WMLES) присоединенных течений авторами DES не предполагалось. Однако попытка такого рода использования SA DES, предпринятая в 2000г. в работе [39], оказалась достаточно успешной: полученные в ней результаты расчета установившегося течения в плоском канале были сопоставимы по точности с аналогичными результатами “специализированных” WMLES подходов. В частности, было показано, что при увеличении числа Рейнольдса DES не требует измельчения сетки в направлениях параллельных стенке или, иными словами, допускает использование неограниченных значений соответствующих шагов сетки, выраженных в единицах “закона стенки” дефект, а именно вместо одного логарифмического участка, формируются два таких участка, имеющих одинаковый (правильный) наклон, но несколько сдвинутых друг относительно друга.
Рисунок 1.32. Профиль скорости при установившемся течении в плоском канале, полученный с помощью SA DES в работе [39] Данный дефект, хорошо видный на рис. 1.32, не только приводит к большой (до 15%) ошибке в расчете трения на стенках канала, но и к сильному (до 65%) завышению локального наклона профиля скорости du / dy в окрестности точки естественно, недопустимо для такого простого течения.
Возникновение двух сдвинутых относительно друг друга логарифмических участков на профиле скорости (этот эффект получил в литературе название Log-Layer Mismatch LLM) характерно не только для DES, но и для других гибридных WMLES моделей, и объясняется следующим образом. Внутренний логарифмический участок формируется в RANS области, благодаря тому, что все “хорошие” RANS модели настроены таким образом, чтобы гарантировать наличие такого участка. Внешняя же часть профиля скорости формируется в области действия LES, который начинает работать “правильно” (то есть воспроизводить логарифмический закон) лишь тогда, когда расстояние от стенки до рассматриваемой точки становится много большим, чем все три пространственных шага сетки. Между двумя этими логарифмическими участками возникает переходный участок, и они оказываются сдвинутыми друг относительно друга.
Устранение LLM является исключительно важной задачей, поскольку ее решение открывает возможность применения WMLES при расчете сложных турбулентных течений на доступных уже в настоящее время компьютерах. Еще более привлекательной представляется идея создания комбинированного подхода, который функционировал бы как DES (DDES) при расчете течений с отрывом и автоматически переходил в режим WMLES, не приводя при этом к LLM, при расчете присоединенных течений. Эта идея реализована в методе DDES с улучшенным пристеночным LES моделированием (DDES with Improved wall-modeling capabilities или IDDES [40]). Данный метод и примеры его применения представлены в разделе 2.5.
II. Математическая формулировка наиболее распространенных подходов к моделированию турбулентности В данной главе приводится математическая формулировка кратко представленных и проанализированных выше подходов к моделированию турбулентных течений, а именно RANS, LES и DES подходов.
2.1. Уравнения Рейнольдса и полуэмпирические модели Спаларта-Аллмараса и 2.1.1. Уравнения Рейнольдса для сжимаемого газа В случае несжимаемой жидкости уравнения Рейнольдса (Reynolds Averaged NavierStokes или RANS) могут быть легко получены из уравнений Навье-Стокса с использованием следующей процедуры осреднения, предложенной Рейнольдсом в его классической работе [41]:
предполагается достаточно большим по сравнению с временными масштабами всех турбулентных неоднородностей, присутствующих в рассматриваемом течении, и достаточно малым по сравнению с характерным временным масштабом осредненного течения.
Как указывается в [41], процедура осреднения предполагает выполнение следующих естественных условий:
где f и g - произвольные функции, которые могут быть представлены как суммы средних и пульсационных переменных константа, а s - пространственная координата или время.
В случае сжимаемого газа предпочтительным является другой способ осреднения (осреднение по Фавру), при котором плотность и давление p осредняются по Рейнольдсу, а для остальных переменных вводятся так называемые средневзвешенные значения В результате осредненные уравнения Навье-Стокса для сжимаемого совершенного газа могут быть представлены в виде (знаки осреднения для простоты опущены):
Здесь u - вектор скорости осредненного течения с компонентами u, v и w, m и t E CvT 0.5 u 2 v 2 w 2 - полная энергия газа, H E p / C pT 0.5 u 2 v 2 w 2 его полная энтальпия, qm и qt - молекулярная и турбулентная составляющие вектора плотности теплового потока, T – температура, Cv C p R m – удельная теплоемкость газа при постоянном объеме, C p - удельная теплоемкость газа при постоянном давлении, R =8.31434 Дж/(моль·K) – универсальная газовая постоянная, а m – молярная масса газа.
Величины молекулярных составляющих тензора напряжений и вектора плотности теплового потока в (2.4) определяются соответственно с помощью реологического закона Ньютона и закона Фурье где S и (T ) - коэффициенты молекулярной динамической вязкости и теплопроводности.
турбулентными составляющими тензора напряжений t и вектора плотности теплового потока qt с параметрами осредненного течения, являющимися основными переменными этой системы, неизвестна и должна быть определена с помощью дополнительных соотношений, которые собственно и составляют модель турбулентности.
В случае использования так называемых линейных моделей турбулентности, а именно к этому классу относятся рассматриваемые в дальнейшем модели СпалартаАллмараса (SA модель) и Ментера (SST модель), предполагается, что справедливы так называемые обобщенная гипотеза Буссинеска и закон Фурье турбулентности16, а t - турбулентная теплопроводность.
Таким образом, в данном случае роль модели турбулентности сводится к определению связи между t t и k с параметрами осредненного течения. При этом при определении турбулентной теплопроводности предполагается, что она может быть выражена через турбулентную вязкость с помощью соотношения где Prt – турбулентный аналог числа Прандтля, обычно полагаемый равным константе ( Prt =0.9).
Наконец, для окончательной формулировки задач о расчете турбулентных течений в замыкающими соотношениями, к этой системе необходимо поставить соответствующие граничные (а при решении нестационарных задач и начальные) условия.
На твердых непроницаемых стенках в качестве граничных условий задаются условия прилипания и непроницаемости для скорости и условия первого или второго рода по температуре где n – направление нормали к стенке.
При использовании SA модели кроме того предполагается, что слагаемое с k в соотношении (2.6) можно опустить.
Кроме того, в тех случаях, когда этого требует используемый численный метод, необходимо задать граничное условие на стенке для давления. В качестве такового в задачах о расчете турбулентных течений обычно используется условие равенства нулю производной от давления по нормали к поверхности [42] На “проницаемых” границах расчетной области (“входе” и “выходе”) в различных задачах используются самые разнообразные граничные условия, однако наиболее распространенными при решении задач внешней аэродинамики являются так называемые характеристические граничные условия [43], формулируемые относительно инвариантов Римана. Связь этих инвариантов с основными переменными, u, v, w и T определяется следующими соотношениями:
Здесь a RT / m - скорость звука, C p Cv – показатель адиабаты, а Vn, V 1 и V 2 - локальные значения нормальной и касательных к границе составляющих вектора скорости.
Конкретная форма граничных условий на проницаемых границах зависит от знака и величины проекции вектора скорости на направление внешней нормали к границе Vn, то есть от того, является ли рассматриваемая точка границы “входом” или “выходом” расчетной области, и от локальной величины числа Маха M n Vn / a в этой точке.
Так, если граница является “дозвуковой” ( Vn a ), то на ее входных участках ( Vn 0 ) задаются значения инвариантов I1, I3, I4 и I5, определенные по известным параметрам набегающего потока, а инвариант I2 линейно экстраполируется на границу из внутренних точек расчетной области. На выходных участках дозвуковых границ ( Vn 0 ) задается значение только одного инварианта, I1, а инварианты I2, I3, I4 и I5 экстраполируются из внутренних точек расчетной области.
В случае сверхзвуковых границ ( Vn a ) на входных участках задаются значения всех набегающем потоке), а на выходных – все инварианты экстраполируются на границу из внутренних точек области.
2.1.2 Модель Спаларта-Аллмараса (SA модель) Как уже отмечалось в главе I, модель Спаларта-Аллмареса была опубликована в 1992 г. [10] и, благодаря своей простоте и достаточно высокой точности при расчете многих аэродинамических течений (см. раздел 1.2), быстро стала одной из наиболее популярных моделей турбулентности.
2.1.2.1. Оригинальная версия модели Модель содержит одно дифференциальное уравнение переноса Модифицированная кинематическая турбулентная вязкость, относительно которой записано это уравнение, и “истинная” турбулентная вязкость t t /, входящая в выражение для Рейнольдсовых напряжений (2.6), связаны соотношением:
где а - молекулярная вязкость газа.
Генерационный и деструктивный (диссипативный) члены в правой части (2.12), P и D, определяются следующими выражениями:
в которых d w - ближайшее расстояние до твердой стенки, - постоянная Кармана, величина представляет собой модуль тензора завихренности а функция f t 2, обеспечивающая подавление так называемого “спонтанного” или “численного” перехода от ламинарного режима течения в пограничном слое к турбулентному, определяется выражением Наконец, последний член в правой части (2.12), предназначенный для инициирования ламинарно-турбулентного перехода в заданной точке (в случае двумерных течений) или на заданной линии (в случае трехмерных течений) на обтекаемой поверхности (“tripterm”)17, рассчитывается с помощью формул:
g t min 0.1, Здесь нижний индекс “trip” относится к величинам, определяемым в той точке на линии перехода, которая находится на минимальном расстоянии от рассматриваемой точки потока, lt ( x 2 z 2 )1 / 2 – диагональ ячейки сетки на обтекаемой поверхности в Отметим, что наличие этого члена заметно облегчает имитацию перехода к турбулентности в рамках полуэмпирических моделей турбулентности, и что SA модель является единственной моделью, в которую с самого начала “официально” был введен такой член.
этой точке, а величина U равна модулю разности векторов скорости в рассматриваемой точке и в точке перехода (в случае неподвижной стенки U u ).
Эмпирические константы модели равны:
2 3, 0.41, cb1 0.1355, cb 2 0.622, Граничные условия к уравнению (2.12) задаются следующим образом.
На твердых стенках модифицированная турбулентная вязкость, как и истинная, должна обращаться в нуль 0.
На входных участках внешней границы расчетной области для должно быть задано условие первого рода, а на выходных - величина экстраполируется на границу из внутренних точек области.
При расчете аэродинамических течений конкретный выбор величины, задаваемой на внешней входной границе, осуществляется следующим образом.
Если точка (линия) перехода от ламинарного к турбулентному режиму течения известна (из эксперимента или предварительного расчета, например, с помощью метода Дрела “ e n ” [44]), то на внешней входной границе расчетной области полагается равной некоторой малой величине (обычно порядка 10-3 молекулярной вязкости), конкретный выбор которой не играет роли. При этом, благодаря наличию в уравнении (2.12) “trip term”, переход от ламинарного к турбулентному режиму течения происходит в малой окрестности заданной точки (линии) перехода.
В тех случаях, когда положение точки (линии) перехода заранее неизвестно, но число Рейнольдса достаточно велико, так что ламинарным участком пограничного слоя на обтекаемом теле можно пренебречь, можно использовать так называемый “полностью турбулентный” (Fully Turbulent или FT) подход к описанию перехода. В рамках этого подхода на входном участке внешней границы расчетной области задается достаточно высокий уровень ( (1 5) ). При этом уравнение (2.12) обеспечивает быстрый (на протяжении нескольких шагов сетки) переход к развитому турбулентному течению в пограничном слое18.
При использовании FT подхода “trip term” в уравнении (2.12) полагается равным нулю.
Наконец, в тех случаях, когда переход к турбулентности происходит не в присоединенном пограничном слое, а в оторвавшемся от обтекаемой поверхности слое смешения19, возможно использование следующей двухстадийной процедуры расчета, получившей название “Trip-Less” (TL) подхода [45].
Расчет начинается в рамках FT-подхода, а затем (после формирования развитой отрывной зоны) производится переход к первому из рассмотренных выше сценариев (то есть относительно большая величина на входных границах заменяется малой величиной) с той лишь разницей, что “trip term” не активируется. В результате, при продолжении расчета турбулентная вязкость “вымывается” из области, расположенной вверх по потоку от отрывной зоны, конвекцией, и течение в этой области становится ламинарным, а в отрывной зоне (и вниз по потоку от нее), благодаря “рециркуляции” турбулентности, течение остается турбулентным.
Отметим, наконец, что для “диагностики” состояния пограничного слоя при проведении расчетов полезно использовать так называемый индекс турбулентности, величина которого в развитом турбулентном пограничном слое близка к 1. Индекс турбулентности вычисляется по формуле [10] где “скорость трения” на стенке u w / в общем случае трехмерного течения рассчитывается как u ( - модуль завихренности в рассматриваемой точке стенки).
Завершая описание оригинальной формулировки SA модели, следует отметить, что в течение почти двадцатилетней истории ее интенсивной эксплуатации, эта модель оставалась практически неизменной: лишь сравнительно недавно к ней было предложено несколько поправок, которые расширили возможность применения модели для течений, расчет которых с ее помощью первоначально не предполагался. Наиболее существенными из этих поправок являются рассмотренные ниже поправка на кривизну линий тока и вращение [46], [47] и поправка на сжимаемость [48].
Такая ситуация может иметь место, например, при докритических режимах обтекания затупленных тел 2.1.2.2. Поправка к SA модели на кривизну линий тока и вращение (SARC модель) Хорошо известно [49], что кривизна линий тока и вращение потока оказывают существенной влияние на характеристики турбулентности и могут приводить как к ее существенной интенсификации, так и к подавлению. Эти эффекты, по крайней мере, в принципе, автоматически описывается в рамках моделей переноса Рейнольдсовых напряжений. В рамках же линейных моделей для их учета необходимо введение специальных поправок. Наиболее удачной из таких поправок является, по-видимому, поправка Спаларта-Шура [46], [47], учитывающая единую природу эффектов кривизны линий тока и вращения и применимая, вообще говоря, к любым линейным моделям.
Однако конкретная калибровка поправки, названной “Rotation-Curvature correction” (RC), была выполнена для SA модели. Соответствующая (при наличии поправки) версия последней получила название SARC модели.
В рамках SARC модели генерационный член P оригинальной SA модели (2.15) модифицируется путем домножения на эмпирическую функцию которая зависит от двух параметров, r * и ~, являющихся мерой кривизны линий тока и вращения. Выражения для этих параметров имеют следующий вид:
Здесь через ( DS / Dt ) ij обозначена соответствующая компонента полной (частной плюс конвективной) производной по времени от тензора скоростей деформаций, а величина D определяется соотношением Три новые константы SARC модели, входящие в (2.24), определяются как Убедительным свидетельством возможностей SARC модели являются результаты ее применения для расчета развитого турбулентного течения во вращающемся плоском канале, для которого оригинальная версия модели, по существу, непригодна (дает те же результаты, что и при отсутствии вращения). В качестве примера на рис. 2.1 сравниваются профили скорости в таком течении, рассчитанные с использованием SA и SARC моделей и в рамках DNS при числе Ренольдса Re=5800 и числе Россби Ro (2H ) / U b 2 ( U b среднерасходная скорость, H - ширина канала, - угловая частота вращения).
Рисунок 2.1. Схема течения в плоском вращающемся канале (слева) и сравнение профилей скорости в таком канале, рассчитанных с использованием SA и SARC моделей, с результатами прямого численного моделирования (DNS) Еще один пример применения SARC модели, позволяющий судить о ее эффективности при расчете течений с сильной кривизной линий тока, представлен на рис. 2.2, где результаты расчета распределения коэффициента трения вдоль внешней стенки U-образного канала, рассчитанные с использованием SA и SARC моделей, сравниваются с экспериментальными данными.
Рисунок 2.2. Схема U-образного канала с расчетной сеткой и сравнение расчетных распределений трения на внешней стенке, рассчитанных с использованием SA и SARC моделей, с экспериментальными данными [50] рассматривавшегося ранее (см. рис. 1.14) вихревого следа за крылом конечного размаха с профилем NACA 0012 и скругленной боковой кромкой. Из рисунка хорошо видно, что в отличие от оригинальной SA модели, которая вследствие завышения турбулентной вязкости предсказывает чрезмерно быструю диссипацию вихрей в следе, SARC модель обеспечивает очень хорошее согласование результатов расчета этого весьма сложного течения с экспериментом.
SA SARC
Рисунок 2.3. Сравнение полей завихренности (верхний ряд) и профилей угловой (средний ряд) и продольной (нижний ряд) составляющих скорости в различных сечениях концевого вихря в следе за крылом конечного размаха с профилем NACA0012 (см. геометрию крыла на рис.1.14а) 2.1.2.3. Поправка к SA модели на сжимаемость потока (SACC) Согласно гипотезе Морковина (см. [51]), при расчете пристеночных пограничных слоев непосредственное влияние сжимаемости на турбулентность начинает сказываться лишь тогда, когда среднеквадратичные пульсации плотности перестают быть малыми по сравнению со средней плотностью потока. Это происходит при числах Маха больше 5, а при меньших числах Маха достаточным оказывается учет эффектов сжимаемости только через учет переменности плотности в уравнениях Рейнольдса. По этой причине оригинальная SA модель, предназначенная, главным образом, для расчета сравнительно низкоскоростных течений, не содержит каких-либо специальных членов, учитывающих влияние сжимаемости потока на характеристики турбулентности. Однако, опыт показывает, что в некоторых течениях, даже при не слишком высоких числах Маха, эффекты сжимаемости все же могут оказывать существенной влияние на турбулентность (приводить к ее подавлению) и, тем самым, на важные характеристики осредненного течения. В связи с этим, в 2000 г. в SA модель была введена поправка на сжимаемость [48], аналогичная такой же поправке в модели А.Н. Секундова t-92 [52]. Она состоит во введении в правую часть уравнения (2.12) следующего дополнительного диссипативного слагаемого где a – местная скорость звука, а эмпирическая константа C5=3.5.Рисунок 2.4. Сравнение распределений давления по радиусу донного среза при сверхзвуковом (M=2.54) обтекании круглого цилиндра, рассчитанных с помощью моделей SA и SACC [54], с экспериментальными данными [53] Характер влияния указанной поправки иллюстрирует рисунок 2.4, на котором представлены некоторые результаты расчетов сверхзвукового (M=2.54) продольного обтекания круглого цилиндра с плоским донным срезом. В частности, из рисунка видно, что введение поправки на сжимаемость позволяет заметно повысить точность расчета уровня донного давления, который является наиболее важной с практической точки зрения характеристикой рассматриваемого течения20.
2.1.3 Модель Ментера (SST модель) Эта модель предложена в 1993 году почти одновременно с SA моделью и, таким образом, опыт ее эксплуатации также насчитывает уже почти 20 лет. Этот опыт свидетельствует о том, что данная модель, получившая название Shear Stress Transport (SST) модели, по совокупности своих качеств является одной из лучших, если не лучшей, среди существующих RANS моделей турбулентности (см. об этом подробнее в разделе 1.1).
SST модель представляет собой комбинацию k- и k- моделей, обеспечивающую сочетание лучших качества этих давно известных моделей. Так, k- модель хорошо зарекомендовала себя при расчете свободных и струйных сдвиговых течений, для анализа которых собственно и была предназначена ее первая версия, предложенная Харлоу [55], а k- модель обеспечивает существенно более точное описание пристеночных пограничных слоев [35]. С учетом этих обстоятельств, Ментером [11] было предложено объединить эти модели с использованием специально сконструированной для этого эмпирической функции F1 (см. далее), которая обеспечивает близость суммарной модели к модели kвдали от твердых стенок и к модели k- в пристеночной части потока.
В первой редакции [11] формулировка этой “гибридной” модели, записанной в терминах k (кинетическая энергия турбулентности) и (удельная скорость ее диссипации), выглядит следующим образом:
Расчет генерационного члена в уравнениях переноса (2.29), (2.30) производится по формуле Предсказываемый SACC моделью волнистый характер изменения давления по радиусу донного среза исчезает при использовании DES (см. [54]).
а последний член в правой части уравнении переноса (так называемый член с перекрестной диффузией – cross-diffusion term) определяется соотношением Для определения турбулентной вязкости по известным значениям k и в SST модели используется не стандартное соотношение t k /, а выражение, базирующееся на известной гипотезе Брэдшоу [56] о пропорциональности напряжение сдвига в пристеночной части пограничного слоя энергии турбулентных пульсаций которое позволяет избежать характерного для k- моделей затягивания отрыва21.
Эмпирическая функция F2, входящая в (2.34), рассчитывается по формуле рассматриваемой точки до ближайшей точки твердой поверхности.
Наконец, эмпирические константы модели определяются через соответствующие константы k- и k- моделей с помощью упоминавшейся выше эмпирической “весовой” функции F1 :
где а CDk max{Dk, 10 20 }.
Индексы “1” и “2” в (2.35) относятся соответственно к константам k - и k- моделей:
Именно введение этого ограничителя (Shear Stress Transport или SST-ограничителя) определило название модели Ментера.
k1 0.85, 1 0.5, 1 = 0.075, а остальные константы равны Граничные условия к уравнениям SST модели (2.29), (2.30) задаются следующим образом.
На твердой стенке кинетическая энергия турбулентности полагается равной нулю, а ее удельная диссипация определяется по формуле где – молекулярная кинематическая вязкость, 1 = 0.075, а y1 – величина первого пристеночного шага сетки.
На входных участках внешней границы расчетной области задается значение удельной диссипации где U и L – характерные для данного течения скоростной и линейный масштабы, а рекомендованные в [11] значения константы C лежат в диапазоне 1 10.
Что касается величины кинетической энергии турбулентности на входных границах, то ее значение k либо задается непосредственно (если оно известно из эксперимента или определено, исходя из каких-то иных физических соображений), либо рассчитывается по величине кинематической турбулентной вязкости на входной границе ( t ), которая предполагается заданной: k ( t ). Как и SA модель, SST модель на протяжение многих лет не подвергалась каким-либо изменениям, и лишь относительно недавно в статье [12], подводящей итоги десятилетнего Значение ( t ), в принципе, может задаваться так же, как и в рамках SA модели (см. раздел 2.1.2.1). Однако, в тех случаях, когда оно мало, модель SST рекомендуется использовать в сочетании с моделью ламинарно-турбулентного перехода [57]. В противном случае, модель SST может приводить к так называемому “вычислительному переходу”, положение которого никак не связано с физическим переходом и зависит от параметров используемой сетки и численного метода.
опыта эксплуатации данной модели, в нее были внесены некоторые незначительные изменения.
Так, при вычислении генерационного члена и члена с перекрестной диффузией в уравнении переноса кинетической энергии турбулентности (2.32) предлагается заменить ограничители в выражении для функции F1 (2.37):
а также несколько видоизменить формулу для расчета турбулентной вязкости (2.33), заменив в ней величину завихренности на скорость деформации S :
Следует, однако, отметить, что при расчете большинства течений эти изменения не оказывают сколько-нибудь значительного влияния на получаемые результаты.
2.2.1. Основные уравнения Уравнения, лежащие в основе LES, могут быть выведены из уравнений Навье-Стокса путем представления всех переменных в виде суммы крупно- и мелкомасштабной составляющих и применения к полученным в результате уравнениям процедуры фильтрации (1.1). Опуская детали вывода, которые можно найти во многих монографиях и учебниках (см., например, [58], [59]), приведем лишь окончательную форму основных уравнений LES для сжимаемого газа:
Единственное формальное различие между уравнениями (2.44) и уравнениями Рейнольдса (2.4) состоит в различии обозначений дополнительных слагаемых, появляющихся в правых частях уравнений движения и энергии в результате пространственной фильтрации конвективных членов уравнений Навье-Стокса (LES) и их временного осреднения (RANS): в уравнениях LES эти члены имеют индекс “SGS” (subgrid), а в уравнениях RANS – индекс “t” (turbulent). Однако, за этими разными обозначениями скрывается принципиально разное физическое содержание LES и RANS подходов к описанию турбулентных течений (см. раздел 1.1), которое непосредственно проявляется на стадии замыкания системы (2.44), то есть при построении моделей подсеточными).
Как и модели для RANS, подсеточные модели, как правило, базируются на использовании обобщенных гипотезы Буссинеска и закона Фурье однако на этом сходство между двумя подходами заканчивается, поскольку соотношения, используемые для определения подсеточной вязкости, принципиально отличаются от моделей турбулентной вязкости для уравнений Рейнольдса.
В качестве примера подсеточной модели в следующем разделе рассматривается первая и до сих пор наиболее популярная модель подсеточной вязкости Смагоринского [60].
2.2.2. Подсеточная модель Смагоринского Данная модель является “подсеточным аналогом” алгебраической модели пути смешения Прандтля для уравнений Рейнольдса. В ее основе лежит предположение (во многом сходное с соображениями Колмогорова) о том, что подсеточная вязкость v SGS определяется средним значением скорости диссипации энергии турбулентности, приходящейся на единицу объема. В этом случае из соображений размерности следует, что где – характерный линейный масштаб фильтра.
Величина скорости диссипации непосредственно не известна, но в случае, когда в энергетическом спектре турбулентности (см. рис. 1.1) имеется отчетливый инерционный интервал, она также может быть выражена с использованием соображений размерности через линейный масштаб фильтра и среднюю скорость деформации S 2 Sij Sij :
Формулировка модели Смагоринского получается подстановкой (2.47) в (2.46):
где CS – эмпирическая константа (константа Смагоринского).
Из выражения (2.48) видно, что в отличие от турбулентной вязкости, подсеточная вязкость зависит не только от параметров отфильтрованного течения (компонент тензора скоростей деформаций), но и от размера фильтра. Как отмечалось в разделе 1.1, в большинстве практических приложений LES явная процедура фильтрации не применяется, а роль фильтра играет используемая расчетная сетка. При этом в качестве обычно используется величина корня кубического из объема ячейки сетки или другие комбинации локальных шагов сетки в трех пространственных направлениях. При использовании сеток близких к изотропным (именно такие сетки являются в большинстве случаев оптимальными для LES) конкретный способ определения не играет роли, так как все они оказываются практически эквивалентными. Однако на практике, в силу особенностей геометрии рассматриваемых течений, часто приходится использовать сильно анизотропные сетки (в первую очередь, это относится к структурированным сеткам), и выбор того или иного определения становится весьма важным. Тем не менее, при измельчении расчетной сетки дополнительные по сравнению с уравнениями НавьеСтокса слагаемые в уравнениях LES (2.44) в любом случае уменьшаются, и решение LES асимптотически стремится к решению DNS. В этом состоит принципиальное отличие метода LES от метода RANS, в котором измельчение сетки приводит лишь к получению “точных” (независящих от сетки) решений уравнений Рейнольдса.
Еще одной важной особенностью подсеточных моделей для LES является то обстоятельство, что входящие в них эмпирические константы (например, константа Смагоринского в модели (2.48)) зависят, вообще говоря, от используемого для решения задачи численного метода. Это объясняется тем, что точность разрешения крупномасштабных вихревых структур в LES зависит не только от сетки, но и от свойств метода, в частности, от присущей ему численной диссипации. Иными словами, численная диссипация метода сама по себе играет роль своеобразной подсеточной модели. Таким образом, если эта диссипация велика, то константа подсеточной модели должна быть соответственным образом уменьшена, а если мала, то, наоборот, увеличена.
В связи с этим, строго говоря, для каждого численного метода должна проводиться индивидуальная калибровка константы Смагоринского или аналогичных ей констант других подсеточных моделей. Такая калибровка осуществляется обычно путем решения задачи о вырождении однородной изотропной турбулентности с использованием различных значений константы и подбора такого ее значения, при котором расчетный спектр разрешенной кинетической энергии турбулентности наилучшим образом согласуется с экспериментальными данными или результатами DNS и подчиняется закону “–5/3” в инерционном диапазоне волновых чисел. В качестве примера, на рис. 2. показаны энергетические спектры, полученные при решении данной задачи с использованием схем с различным порядком аппроксимации невязких составляющих газодинамических потоков и фиксированном значении константы C S.
Рисунок 2.5. Спектры разрешенной энергии турбулентности в задаче о вырождении однородной изотропной турбулентности, полученные из LES с моделью Смагоринского с использованием различных численных методов при CS 0. Из него видно, что численный метод оказывает сильное влияние на результаты LES, использовании только одной из трех рассматриваемых схем аппроксимации невязких потоков, а именно симметричной схемы четвертого порядка точности. В случае использовании противопоточной схемы коротковолновые моды спектра практически не разрешаются из-за слишком высокой численной диссипации даже при достаточно высоком (пятом) порядке точности схемы. В противоположность этому, при использовании симметричной схемы второго порядка точности (для невязких течений эта схема является неустойчивой) с тем же значением константы диссипация модели и метода оказывается недостаточной для обеспечения правильного описания каскада турбулентности (энергия турбулентности “генерируется” и накапливается в коротковолновых модах спектра).
Рисунок 2.6. Спектры разрешенной энергии турбулентности в задаче о вырождении однородной изотропной турбулентности, полученные из LES с моделью Смагоринского: a) – влияние выбора константы на результаты расчета на сетке 323; b) – влияние сетки на результаты расчетов при CS 0. На рис.2.6 представлены энергетические спектры, полученные при различных значениях C S с помощью NTS кода, в котором для LES используется симметричная аппроксимация невязких потоков с четвертым, а вязких – со вторым порядком точности.
Из него следует, что оптимальным значением константы Смагоринского для этого кода является C S 0.2.
Отметим в заключении, что модель Смагоринского (2.48) не обеспечивает равенства нулю подсеточной вязкости на твердой поверхности и поэтому не может непосредственно применяться для расчета пристеночных течений. Для устранения этого недостатка в нее вводится демпфирующий множитель [61], являющийся аналогом множителя Ван Дриста в модели Прандтля для RANS:
Кроме того, опыт показывает, что при расчете пристеночных течений необходимо использовать примерно в два раза меньшее значение константы Смагоринского, чем при расчете свободных турбулентных течений. Это факт иллюстрирует рис. 2.7, на котором сравниваются результаты расчета установившегося течения в плоском канале, использованием NTS кода при C S 0.2 и 0.1. Из него видно, что в первом случае, то есть при использовании значения C S 0.2, полученного при калибровке этой константы на задаче о затухании свободной изотропной турбулентности, результаты расчета существенно отличаются от соответствующих результатов DNS [62], а при C S 0.1 практически совпадают с этими результатами.