CC BY
53
Том XLII
УЧЕНЫЕ ЗАПИСКИ ЦАГИ 2011
№ 4
УДК 519.688, 533.6
АЛГОРИТМ ИНИЦИИРОВАНИЯ ЛАМИНАРНО-ТУРБУЛЕНТНОГО ПЕРЕХОДА ПРИ ЧИСЛЕННОМ МОДЕЛИРОВАНИИ ТЕЧЕНИЯ НА БАЗЕ УРАВНЕНИЙ РЕЙНОЛЬДСА
В. В. ВЛАСЕНКО, А. Н. МОРОЗОВ
Статья посвящена разработке «инициатора перехода» — численного алгоритма, который может быть использован в расчете течений с пограничными слоями на базе уравнений Рейнольдса. Этот алгоритм «выключает» источниковые члены модели турбулентности вверх по потоку от поверхности, на которой начинается ламинарно-турбулентный переход, и постепенно «включает» их вниз по потоку от этой критической поверхности. Описан поиск удовлетворительной версии данного алгоритма. Дано краткое объяснение идей, лежащих в основе каждой версии инициатора, и причин, по которым была забракована каждая из промежуточных версий. Также описан устойчивый алгоритм оценки толщины пограничного слоя в неоднородных течениях.
Ключевые слова: вычислительная аэродинамика, уравнения Рейнольдса, пограничный слой, ламинарно-турбулентный переход, полуэмпирическая модель турбулентности.
Во многих задачах практической аэродинамики (например, при определении характеристик крыла самолета, лопаточных устройств в турбореактивном двигателе и др.) важное значение имеет правильное предсказание положения ламинарно-турбулентного перехода (ЛТП).
При численном решении задач аэродинамики, как правило, используются уравнения Рейнольдса, замкнутые какой-либо полуэмпирической моделью. К сожалению, обычные модели турбулентности непригодны для предсказания положения ЛТП. Поэтому при численном решении задач, где положение ЛТП имеет принципиальное значение, необходимо добавить к обычной модели турбулентности полуэмпирическую модель перехода, которая должна «включить» обычную модель турбулентности в том месте, где начинается ЛТП.
Такие полуэмпирические модели ЛТП известны [1]. Наиболее продвинутой моделью данного класса является модель Лэнгтри — Мен-тера [2], которая предназначена для использования в сочетании с моделью турбулентности SST Ментера [3] и включена в коммерческий пакет ANSYS CFX. К сожалению, все модели ЛТП этого класса носят глубоко эмпирический характер. Они основаны на абстрактных аппроксимациях экспериментальных данных, а не на физических закономерностях. Поэтому авторы настоящей работы предложили новый подход к построению полуэмпирических моделей ЛТП [4, 5]. Этот подход основан на решении дифференциального уравнения для амплитуды возмущений
скорости 0 = max V , где максимум берется по различным возмущениям в данной реализа-
ВЛАСЕНКО Владимир Викторович
кандидат физикоматематических наук, начальник сектора ЦАГИ
МОРОЗОВ Александр Николаевич
кандидат физикоматематических наук, старший научный сотрудник ЦАГИ
ции процесса возмущений. Уравнение для 9 может быть выведено из линеаризованных уравнений Навье — Стокса; затем к источниковым членам этого уравнения применяется оператор максимизации, который выделяет из множества реализаций ту, что растет быстрее всего. Применение этого оператора связано с вычислением средних по реализациям характеристик возмущенного движения. Для оценки этих средних величин используются распределения возмущений в волнах неустойчивости, которые предсказываются классической теорией устойчивости пограничного слоя (параметрические аппроксимации собственных функций уравнений типа Орра — Зоммерфельда [6]). В конце концов, для амплитуды возмущений скорости в наиболее быстро растущей реализации процесса возмущений получается уравнение следующего вида:
<Эр9 д
дх.
-а~ 50
р9мг■ -ц—
дх.
= А + В-в.
(1)
^0
где коэффициенты А, В являются сложными функциями от ии —-----------,----. Уравнение (1) со-
дх]- дх]- дх]-
держит 2 + 2М эмпирических констант, где N — количество учитываемых физических механизмов перехода. В работах [4, 5] рассматривался ЛТП на поверхности стреловидного крыла, и там учитывались волны Толлмина — Шлихтинга и стоячие волны боковой неустойчивости; поэтому N было равно 2. Невязкий механизм продольной неустойчивости (неустойчивость Кельвина — Гельмгольца) учитывается в любом случае и не требует введения эмпирических констант.
Уравнение (1) следует решать численно совместно с системой уравнений Рейнольдса, замкнутой какой-либо обычной моделью турбулентности. Предполагается, что переход начинается там, где амплитуда возмущений скорости достигает некоторого критического значения 9кр (это
критическое значение определяется условием, чтобы нелинейные члены в уравнениях Навье — Стокса для возмущений стали сопоставимы с линейными). Вниз по потоку от места начала перехода «включаются» источниковые члены обычной модели турбулентности. Численный алгоритм, который «включает» обычную модель турбулентности, ниже будет называться «инициатором перехода».
Инициатор перехода имеет самостоятельное значение: в принципе, его можно использовать не только в сочетании с моделью ЛТП, разработанной авторами, но в сочетании с любой моделью, которая дает форму критической поверхности, где начинается нелинейная стадия роста возмущений. Настоящая статья посвящена созданию такого алгоритма-инициатора.
С точки зрения авторов, методический интерес могут представлять не только конечное состояние этого алгоритма, но и промежуточные версии, которые возникали по мере его создания. Будет дано краткое объяснение идей, лежащих в основе каждой версии инициатора, и причин, по которым была забракована каждая из промежуточных версий.
Нужно добавить, что описываемый алгоритм ориентирован на расчеты стационарных течений. Он также предполагает, что используются регулярные расчетные сетки с шестигранными ячейками. С точки зрения авторов, именно такие сетки предпочтительны для достижения достаточного качества описания области пограничного слоя.
1. ОБЩАЯ СТРУКТУРА ИНИЦИАТОРА ПЕРЕХОДА
В настоящей работе в качестве примера для замыкания системы уравнений Рейнольдса будет использоваться модель турбулентности SST Ментера [3]. В этой модели замыкание уравнений Рейнольдса проводится на базе гипотезы Буссинеска; кинематическая турбулентная вязкость выражается через два параметра турбулентности — кинетическую энергию турбулентности к = и\и\!2 и характерную частоту турбулентных пульсаций оз. В области небольших градиентов среднего течения V, = к/со. Параметры к и со находятся из решения дифференциальных уравнений следующего вида:
Эр к д
Э? Эх,
Эре» Э д1 Эх,
р ки1 — р
у + -
Рг;
дк
Эх,
(
р ешг - р
у + -
Рг®
гх? У
Эсо
Эх,
где ^ и Л'(|| — источниковые члены, которые связаны с производством кинетической энергии турбулентности и ее диссипацией, а также учитывают некоторые другие локальные эффекты [3].
Уравнения (2) решаются совместно с системой уравнений Рейнольдса. Одновременно решается и уравнение (1) для амплитуды возмущений скорости 0. Из решения этого уравнения в каждый момент времени известна форма критической поверхности, на которой 0 впервые достигает значения 9кр. Для инициирования ЛТП предлагается умножить источниковые члены в уравнениях (2) на инициирующую функцию п, которая равна нулю в области ламинарного течения, равна единице в области развитой турбулентности и меняется от нуля до единицы в области ЛТП:
$к ~~^/лТП $к’ $о> ^/лТгДо-
В первоначальной версии инициатора перехода функция /шп зависела от двух переменных — 9 (амплитуды возмущений скорости) и й?лтп (расстояния от данной точки пространства до ближайшей точки критической поверхности):
/лтп ^ЛТП -
о.
0<0
кр>
0...1, 0^0кр, ^лтп<^
(3)
1,
'-®кр= ^ЛТП-^
где X — длина области Л111. В настоящей работе для X используется эмпирическая формула Дэя и Нарасимхи [1, 7]:
Х =
0.411-
Яс.
К
кр V
(4)
где гКр —точка, в которой 0 = 0кр; V —кинематический коэффициент молекулярной вязкости; Ут — средняя скорость газа в направлении, касательном к поверхности обтекаемого тела;
К/5** _ **
т ; 5 — толщина потери импульса пограничного слоя; индексом е обо-
Яе5** гкр -
■кр
значены параметры на внешней границе пограничного слоя. Параметр Пр в формуле (4) определяется следующим образом [1, 7]:
I-1.453-10 3In0.2-1.6L10-3,
Ти > 0.2%,
■жр
I-1.453 • 10~31п Ти -1.61-10 3, Ти < 0.2%,
(5)
Л/ И', М-да/З
где Ти,, = —-----------100% — уровень возмущений в набегающем потоке. Предполагая, что
^ОО
возмущения скорости в набегающем потоке равнораспределены по направлениям, можно пока-1 0
зать, что Ти = —7=—^ -100%. Формула (5) должна корректироваться в случае быстрого продоль-л/З
**
ного изменения параметров течения вне пограничного слоя (см. [1, 7]). В данной работе этот эффект не учитывается.
2. ТЕСТОВАЯ ЗАДАЧА О ЛТП НА ПЛАСТИНЕ
В качестве основной тестовой задачи здесь будет рассматриваться задача о ЛТП в пограничном слое на плоской пластине. Режим течения: число Маха Ме =0.8; статическое давление ре = 1 атм; статическая температура Те = 288 К; начальный уровень турбулентности Ти., =0.1%. Длина пластины Ь — 0.5 м. Согласно классическим данным Шубауэра и Скрэмстеда
УвХ 6
[8, 9], при Тиу = 0.1 % переход на пластине начинается при Ясл. = ——«2.8-10 и за-
нач
канчивается при Ясл. кон ~3.9-106, откуда координаты начала и конца ЛТП, отсчитываемые от начала пластины, равны соответственно хнач«0.15 ми хкон«0.21 м.
Для пограничного слоя на плоской пластине область ЛТП можно идентифицировать как
область роста коэффициента трения су =------—----. На ламинарном участке су х ~ х 1 2. а по
0.5ре Vе
окончании Л111 Су х ~ х х0 ' 5. Между этими областями существует участок гладкого роста су х . Длина этого участка и будет считаться длиной области перехода X.
Расчеты проводились на сетке, которая содержала 196 ячеек вдоль пластины и была сгущена в продольном направлении к началу и к концу пластины (наименьший шаг сетки вдоль
оси х — 4-10-4 м, наибольший — 5-10-3 м) и по вертикали (наименьший шаг сетки —
2.3-10-6 м, у поверхности пластины). В конце пластины толщина пограничного слоя составляла
примерно 4-10 м, и поперек пограничного слоя помещалось около 40 ячеек.
Пластина считалась теплоизолированной поверхностью с прилипанием потока. Свободные части нижней границы расчетной области рассматривались как плоскости симметрии. На внешней границе расчетной области ставились мягкие граничные условия, основанные на анализе инвариантов Римана.
Численный метод, который использовался в расчетах, описан в [10].
3. РАБОТА ИНИЦИАТОРА ПРИ ПРЕДПИСАННОМ ПОЛОЖЕНИИ ПЕРЕХОДА
Вначале были выполнены расчеты с предписанным положением перехода. Вместо решения уравнения для 9 было принято, что 9 <9кр при х<0.15 ми 9>9кр при х>0.21 м. На рис. 1
сравниваются типичные поля продольной скорости, полученные по стандартной модели 88Т (2) и по модели, в которой источниковые члены Л'(11 были умножены на инициирующую функцию (3). Для визуализации пограничного слоя рисунок растянут вдоль вертикальной оси в 573 раза. При использовании модели SST без инициатора переход происходит слишком близко к передней кромке пластины. С инициатором переход начинается в предписанном сечении. Структура полностью развитого турбулентного пограничного слоя одинакова в обоих случаях.
На рис. 2, а сравниваются распределения с/ х , соответствующие полям течения, показанным на рис. 1, а, б. Там же показано распределение х , полученное с использованием модели турбулентности (^ — со) Коукли [11]. При использовании этой модели получается заметный участок ламинарного пограничного слоя, но положение ЛТП предсказывается некорректно.
Рис. 1. Поля продольной скорости [м/с], полученные в расчетах пограничного слоя на пластине по стандартной модели ББТ (а) и с использованием инициатора с предписанным положением перехода (б)
Рис. 2. Распределения коэффициента трения вдоль пластины: а — кривые, полученные с использованием различных моделей турбулентности; б — кривые, полученные с использованием различных инициирующих функций
На этом этапе был проведен предварительный выбор инициирующей функции. На рис. 2, б сравниваются распределения е^х , полученные с использованием в качестве /лтп различных функций. Все эти функции удовлетворяют условиям: /лтп — 0 при 9 < 9кр и /лтп — 1 при 9>9кр, б/ЛТ1[ >X. Вид каждой функции при 9>9кр, с1\\\\ \ <Х приведен в таблице. Для сравнения приведены также экспериментальные данные Шубауэра и Клебанова (1955 г.) для обтекания пластины при близком уровне турбулентности (Ти = 0.18%), взятые из статьи [2]. Экспериментальную зависимость пришлось отмасштабировать вдоль оси е^ так, чтобы последняя экспериментальная точка легла на расчетные кривые. По-видимому, отличие е^ на турбулентном участке от расчетов связано с тем, что эксперимент проводился для М ~ 0.15, а расчеты для М~0.8.
1. Функция «Peak» f 1 + COS 271 djYYYl/h — 0.5 , С^ЛТП ~ '^'ЛТП I 3 + cos 2п djijn/X-0.5 Jl, Х/2 <^ЛТП<Я,
2. Функция «Linear» /лгп = ^лтп А
3. Функция «Linear-0.5» /лгп = 1 + ^лшА /2
4. Функция «Cosinus» /лтл = 1 + COS 71- ^дТпА-1 /2
5. Функция «Arcsinus-0.3» /лШ =0.3 + 0.7. [arcsin 2^лш/А,-1 J/ti + 0.5
6. Функция «Sinus-0.3» /лш =0.3+ 0.7-sin O.SndjjjYi/X
7. Функция «Sinus» /лгп = sm 0.5ndmii/X
Распределение Cf x , которое по форме наиболее напоминает экспериментальную зависимость, было получено с использованием функции следующего вида:
0, 0<екр,
сс+ 1-сс sin 0.5лс/лтп/А, , 9>9кр, с/ЛТП/Д<1, (6)
1, ^лтп А -1 •
Параметр а был на этом этапе принят равным 0.3. Отметим, что при а —0 также получаются неплохие результаты.
4. МОДИФИКАЦИЯ ИНИЦИАТОРА ДЛЯ РАБОТЫ В СОЧЕТАНИИ С ПОЛУЭМПИРИЧЕСКОЙ МОДЕЛЬЮ ЛТП
Далее инициатор с выбранной функцией /лхп (формула (6), а = 0.3) был применен к расчету той же тестовой задачи совместно с решением уравнения для амплитуды возмущений 9. Распределение параметра IgO, которое получается в данной задаче, изображено на рис. 3. Как показывает график на врезке, вдоль выделенной линии амплитуда возмущений растет почти экспоненциально. Однако первые результаты оказались неудачными. На рис. 4, а показано поле продольной скорости. (Все изображения на рис. 4 растянуты вдоль вертикальной оси в 573 раза, как и на рис. 1.) Хорошо видно, что на турбулентном участке пограничный слой испытывает незатухающие продольные колебания. При этом его толщина растет гораздо медленнее, чем на рис. 1, б.
На рис. 5 показаны распределения Cf x , полученные в разных расчетах. Зависимость CfX , полученная с предписанным положением перехода, будет использоваться в качестве эталонной. На турбулентном участке кривая, соответствующая рис. 4, а («Вариант 1»), идет гораздо ниже эталонной кривой и испытывает колебания.
Анализ полей течения позволил установить причину этих плохих результатов. На рис. 6 представлено поле lg 9, соответствующее полю продольной скорости, показанному на рис. 4, а, но растянутое слабее вдоль вертикальной оси (в 57 раз). Там же показана критическая поверхность. Видно, что эта поверхность сильно искривлена и вытянута в продольном направлении. Если в инициирующей функции (3) рассчитывать расстояние ^дтп , как расстояние от текущей точки до ближайшей точки критической поверхности, то для большинства точек ближайшая из них помещается над текущей точкой, на расстоянии, которое намного меньше, чем толщина пограничного слоя (см. рис. 6). Но длина перехода 0.06 м» 5. Следовательно, везде
б/ди [ /Х<с 1. В такой ситуации функция /Jrn ( <sc 1, а источниковые члены модели турбулент-
ности (2), и . умножаются на это малое число. В результате турбулентность развивается намного медленнее, чем в действительности. Это приводит не только к замедленному росту толщины пограничного слоя, но и к малой турбулентной вязкости. Сглаживание, обусловленное турбулентной диффузией, оказывается недостаточным и не позволяет гасить возмущения, связанные с переходом внешней границы пограничного слоя с одного горизонтального ряда ячеек сетки на следующий.
Но в реальности возмущения развиваются практически вдоль линий тока основного течения. Следовательно, расстояние ^ЛХП должно рассчитываться не так, как предложено в разделе 1, а вдоль линий тока.
В принципе, проследить линии тока возможно; но алгоритм трассировки трудно сделать устойчивым и достаточно точным. Вместо этого была предложена идея получать значения функции из решения
дифференциального уравнения. Заметим, что форма зависимости (6) напоминает кривую релаксации к стационарному состоянию, описываемую дифферен-
Рис. 3. Поле десятичного логарифма амплитуды возмущений скорости 0 на ламинарном участке пограничного слоя (на врезке — распределение ^ 0 вдоль выделенной линии)
д/лтп , д/птп.
- Uj
циальным уравнением — = р 1 - / . Р = const >0.
ds '
Поэтому было предложено находить /лтп из решения следующего уравнения:
dt
дХ:
0,
К<
кр>
I VI
1-/лтп
0 > 0
(7)
кр>
где | V | — модуль скорости среднего течения. В набегающем потоке было предложено брать /лтп = 0.
При 0>0кр в стационарном течении уравнение (7) эквивалентно уравнению
^/лтп _ 4.1
l-/j
лтп
ds X является следующая функция:
где 5 — координата вдоль линий тока. Если X = const, то решением (8)
ЛТП
I0’
1-exp -4.1s ,
кр>
'кр-
(8)
Функция (8) близка к функции (6), если в (6) положить а = 0.3 и вычислять с/лтп вдоль линий тока. При этом функция (8) непрерывна в точке 0 = 0.
Результаты расчетов с решением уравнения (7) представлены на рис. 4, б. Теперь турбулентность развивается более интенсивно. Благодаря росту турбулентной вязкости сглаживание существенно усилилось, и колебания толщины пограничного слоя почти исчезли. Однако скорость роста пограничного слоя не увеличилась (соответствующая кривая на рис. 5 — «Вариант 2» — близка к колеблющейся кривой, полученной в Варианте 1).
Для объяснения этого результата сравним поля кинетической энергии турбулентности, полученные в расчете с предписанным положением перехода и в Варианте 2 (рис. 7,а,б). Видно, что в Варианте 2 турбулентность растет лишь в узком слое ниже уровня у = 0.0005 м. Форма
е)
Рис. 4. Поля продольной скорости [м/с], полученные в расчетах пограничного слоя на пластине с использованием различных версий инициатора перехода
Рис. 5. Распределения коэффициента трения вдоль пластины (варианты 1 —6 соответствуют рис. 4, а — е)
этого слоя примерно соответствует форме критической поверхности, поскольку, согласно уравнению (7), рост функции /лтп происходит только вдоль линий тока. За пределами этого узкого слоя по-прежнему /дХП <§с1, т. е. источники модели турбулентности занижаются. Это и препятствует расширению области турбулентного течения в вертикальном направлении.
В реальности, очевидно, после возникновения неустойчивости течения возмущения начинают развиваться уже не только вдоль линий тока среднего течения, но и в поперечном
О 0.125 0.250 0.375 0.500
Рис. 6. Поле десятичного логарифма амплитуды возмущений скорости, полученное в Варианте 1
Рис. 7. Поля кинетической энергии турбулентности Цм2/с2полученные в расчетах пограничного слоя на пластине с предписанным положением перехода (а) и в Вариантах 2—6 (б — е)
направлении — за счет турбулентной диффузии. Поэтому было решено добавить в уравнение (7) члены, описывающие молекулярную и турбулентную диффузию:
др/]
лтп
_д_
Эх,.
р/
ЛТПМг
лтп
Эх,-
0,
0 < 0
IV \
1 /лтп
|>|
кр>
кр *
(9)
Число Прандтля для параметра /лтп было принято равным Рг/ = 1. Результаты показаны на рис. 4, в; 5 (кривая «Вариант 3») и 7, в. Очевидно, что структура течения стала намного лучше. Скорость роста пограничного слоя существенно увеличилась. И колебания толщины пограничного слоя практически пропали, так как турбулентная диффузия во внешней части пограничного
слоя усилилась. Хем не менее, скорость роста пограничного слоя в турбулентной части еще занижена. Соответственно, и коэффициент трения с/ х на рис. 5 повысился примерно в 2 раза,
но он все еще заметно ниже, чем на эталонной кривой.
Причина пониженной скорости роста пограничного слоя заключается в том, что новый член в уравнении (9), описывающий турбулентную диффузию, сглаживает распределение /Лхп и, следовательно, уменьшает его максимальное значение. В результате /лп ( < 1 везде, и естественное развитие турбулентности нарушается.
Поэтому было решено отказаться от уравнения для /лхп. Вместо этого была сделана попытка решать дифференциальное уравнение для аргумента инициирующей функции — параметра ф = б/ЯП[/Х. Было предложено следующее уравнение:
Эрф д дґ дхі
рфы,- -р
Ргф гх? У
дхі
о,
рі^і 0>е <|0>
, и _ У , шах к ^
с граничным условием в набегающем потоке <рг = 0. Число Прандтля для параметра ф было принято равным Рг/ =1. Длина перехода X вычислялась по формуле (4); однако для простоты организации вычислений в эту формулу подставлялись не параметры в критическом сечении, а местные характеристики пограничного слоя. При таком упрощении формула (4) дает X —> 0, если толщина пограничного слоя 5—>0. Чтобы устранить деление на ноль, в правую часть уравнения (10) введен малый параметр г- > 0. Можно принять параметр равным наименьшему из габаритов ячеек сетки. Расчет пограничного слоя должен проводиться на сетке, которая содержит несколько сильно вытянутых ячеек на протяжении длины перехода X. Поэтому в области перехода ех <§; X.
В стационарном решении в области ЛТП, где турбулентной диффузией еще можно пренебречь, уравнение (10) практически сводится к дифференциальному уравнению = —, ф 0 =0
с1я X
(^ — координата вдоль линий тока). Поскольку в пределах области ЛТП А,«сопз1 = А, гкр , то
5 б/ ттттт
там ф і
_ ^лтп
После того как параметр ф определен из решения уравнения (10), инициирующую функцию можно вычислить по формуле (6), в которую вместо й^лтп/^ подставляется ф. Предполагалось, что рано или поздно с/лтп станет больше X. ф = й?лтп /X станет больше 1 и, соответственно, ниже по потоку всюду будет /лхп = 1.
Учитывая успех непрерывной при б/ЛТ1 ( = 0 функции (8), было принято решение в дальнейшем использовать функцию (6) с а = 0, которая также непрерывна при б/лп ( = 0.
Результаты расчетов с использованием уравнения (10) показаны на рис. 4, г; 5 (кривая «Вариант 4») и 7, г. Против ожидания, качество результатов ухудшилось. Зависимость С/ х
стала немонотонной на турбулентном участке. Но более интересно понять, почему не увеличилась скорость роста толщины пограничного слоя. Это означает, что /Лтп по-прежнему меньше 1, ЧТО ВОЗМОЖНО ЛИШЬ при ф = б/|ГГ][/л <1.
Было найдено следующее объяснение. Как уже было сказано выше, для вычисления X использовалась формула (4), в которую вместо параметров при г -гкр подставлялись местные ха-
3/2
рактеристики пограничного слоя. Но тогда формула (4) дает зависимость X х ~ 5 ' х . В турбулентной зоне при правильной скорости роста пограничного слоя 5 х ~ х4 5. В рассматривае-
мых расчетах пограничный слой рос медленнее, но все-таки более быстро, чем в ламинарной области. Следовательно, 5 х > Ах1 2 и X х > Вх* 4. Расчеты показали, что в турбулентной области параметр X рос почти линейно вдоль оси х. Но тогда на турбулентном участке некорректно интерпретировать ср, как d]vmjx /;.р , и подставлять ср вместо с/]П1[/'к в функцию Удхп- Из
уравнения = —-— < —-— следует, что ф < б/лти /х г . Турбулентная диффузия парамет-
ds X х X г
^ 'кр
ра ф только занижает максимальное значение этого параметра. Расчеты показывают, что в конце пластины параметр ф достигает лишь значений, меньших 0.5. При использовании формул (2) и (6) это приводит к тому, что источниковые члены в уравнениях для параметров турбулентности занижаются более чем в 2 раза.
Чтобы устранить этот недостаток, нужно установить X = const вниз по потоку от области перехода. Это было реализовано следующим образом.
Во-первых, вместо дифференциального уравнения для t/jrn[//- было записано дифференциальное уравнение для параметра А, который в области перехода должен совпадать с d\Y[U. а дальше вниз по потоку должен, по крайней мере, расти. Это уравнение имеет следующий вид:
<ЭрА д
8t Ох,
рАиг -р
v, дА
V н------
Рг
V rlt
дх
о, е<е
- (П)
Pin
'кр ?
с граничным условием в набегающем потоке Ау = 0. Легко видеть, что в стационарном течении в области перехода, где турбулентной диффузией можно еще пренебречь, это уравнение эквивалентно уравнению ^— = 1, откуда Д = 5=й?лтп. В функцию /дтп вместо йТщ-пА подставляет-
ds
ся А/Х. На турбулентном участке диффузионные члены в (11) приведут к более медленному
росту А в сравнении с <:/jm (: но там будет обеспечено X — const — X гкр . Поэтому там А/Х > 1 и,
соответственно, /дхп = 1 •
Во-вторых, длина перехода X теперь определяется следующим образом:
max ех;Х , А <5,
^ = <| теор , , (12) A r,t , А >8,
где Л г, t — решение дифференциального уравнения + ^Р^- ‘ = о Таким образом, в об-
ft сЦ
ласти, где Д<5 (т. е. до самого начала перехода, включительно), используется «теоретическое» значение X = л|СОр. Оно вычисляется, как и раньше, по формуле (4), в которую подставляются
местные характеристики пограничного слоя. Но при А > 8 величина X сносится вдоль линий тока без изменения.
На рис. 4, д; 5 (кривая «Вариант 5») и 7, д показаны результаты расчета, в котором использовались уравнения (11)и(12). В этом расчете было взято значение Рг(д =1. Видно, что скорость роста турбулентного пограничного слоя выросла. Соответственно, коэффициент трения в турбулентной части заметно увеличился.
И все-таки коэффициент еще ниже, чем на эталонной кривой. Это можно объяснить,
анализируя поле длины перехода X на рис. 8. Жирной линией показана поверхность, на которой достигается А = 8. Вверх по потоку от этой линии X вычисляется по формуле (4), а вниз по потоку от этой линии X сносится вдоль линий тока. Теперь очевидно, что над пограничным слоем ве-
личина /.|СОр успевает достичь довольно большого значения -- ^теор ~ 0.1 м, которое выше, чем А в этой области.
Поэтому в верхней части турбулентного пограничного слоя А/Х < 1 и /дхп < 1.
Для устранения этого эффекта было решено подставлять в функцию /дхп вместо СІЖП/Х значение
^ ^тах
Рассмотрим ряд ячеек вдоль индексной оси, пересекающей твердую стенку. Тогда Атах — максимальное значение А
в пределах этого ряда, а X Атах — значение X в ячейке, где этот максимум достигается. Таким образом, теперь /лтп = const в пределах каждого вертикального ряда ячеек, а величина /ЛТП определяется ячейкой, которая находится на максимальном расстоянии от критической поверхности среди ячеек этого ряда (т. е. ячейкой, где процесс перехода развит более, чем во всех остальных ячейках этого ряда). Из-за этой модификации инициатор перехода становится более близким к инициатору с предписанным переходом, который был рассмотрен в разделе 3 и дал наиболее естественные результаты.
Результаты расчета с этой модификацией показаны на рис. 4, е; 5 (кривая «Вариант 6»)
и 7, е. Толщина пограничного слоя выросла, и коэффициент трения достиг уровня, который был
получен в «эталонном» расчете с предписанным переходом.
Поэтому модификация инициатора перехода, основанная на решении уравнений (11) и (12) А„,„
Рис. 8. Поле длины области ЛТП X [м], полученное в расчетах пограничного слоя на пластине (Вариант 5)
и на использовании
в качестве аргумента функции fmп, была принята в качестве
окончательной версии алгоритма инициирования перехода.
Был произведен ряд дополнительных численных экспериментов, которые показали, что коэффициент РггА можно повысить до значения 100. Это снижает интенсивность «турбулентной диффузии» параметра А в 100 раз. Благодаря этому параметр А оказывается близким к с1ти не только в области ЛТП, но и в области развитого турбулентного течения. Тем не менее, этой малой интенсивности искусственной «турбулентной диффузии» оказывается достаточно для сглаживания возмущений, порождаемых сеткой.
После этого была произведена окончательная настройка эмпирических констант в уравнении для 9, благодаря которой переход стал начинаться в сечении А'нач «0.15 м, как и в экспериментах [8]. Соответствующее распределение с^ х помечено на рис. 5 как «Вариант 7». Эта кривая достаточно близка к эталонной.
5. ПРИМЕНЕНИЕ ИНИЦИАТОРА ПЕРЕХОДА К ДРУГИМ ТЕСТОВЫМ ЗАДАЧАМ
Кроме задачи о ЛТП в пограничном слое на плоской пластине, были решены и другие тестовые задачи — задача о плоском течении около профиля крыла и задача о трехмерном обтекании стреловидного крыла. В обеих задачах рассматривался один и тот же режим обтекания (число Маха набегающего потока =0.78, число Рейнольдса, посчитанное по хорде профиля, —
н
Яес -2-10 ). Угол атаки был равен а = 0.31 °
Здесь необходимо напомнить, что согласование расчетов с экспериментом определяется не столько алгоритмом инициирования перехода, сколько правильностью определения критической поверхности (поверхности, на которой начинается ЛТП), а также правильностью задания длины Л111 X. Эти факторы являются внешними по отношению к алгоритму инициирования
перехода и в настоящей статье не рассматриваются. Поэтому тестовые расчеты профиля крыла и полного крыла не будут описываться здесь подробно. В рамках данной статьи эти расчеты представляют интерес лишь как проверка устойчивости алгоритма при его применении к более сложным течениям, чем течение в пограничном слое на плоской пластине.
И в самом деле, данные расчеты выявили много новых эффектов, которые не встречались в задаче на плоской пластине. Среди этих эффектов стоит указать на сложность определения местной толщины пограничного слоя § (что существенно в обеих задачах) и на необходимость учитывать возможность развития возмущений по механизму боковой неустойчивости (crossflow instability) при обтекании стреловидного крыла.
Трудность определения местной толщины пограничного слоя 8 связана с такими факторами, как сложность идентификации внешней границы пограничного слоя в течении с непараллельными линиями тока, возможность наличия точки перегиба в профиле пограничного слоя и отрыва пограничного слоя, а также наличие особенностей типа точки торможения и острой кромки профиля крыла. При этом алгоритм определения 8 должен работать и для ламинарного, и для турбулентного пограничного слоя. После долгих экспериментов был выбран следующий алгоритм определения 8.
Рассмотрим опять ряд ячеек вдоль индексной оси, пересекающей твердую стенку. Будем нумеровать эти ячейки индексом / (ячейка / = 1 прилегает к твердой стенке). Двигаясь от стенки,
найдем в этом ряду ячеек ячейку г =/тах, в которой величина
dyw
Ji
впервые достигает мак-
симума (У — касательная к стенке компонента скорости, yw — расстояние до стенки). В ячей-
ке i=i„
вычисляются касательное напряжение х = pv
dVr
w
и параметр ут ~v/x. После этого
в каждой ячейке ряда вычисляется координата >’(+ = . /ут. Внешние характеристики погра-
ничного слоя, обозначенные индексом е в формуле (4), вычисляются в точке, где у+ = 3000. Далее вычисляется толщина потери импульса пограничного
слоя
yw
8** —
Р V,
Р%е
-max
( V ^ 1--Ч0
yw . В качестве оценки толщины пограничного слоя вычисляет-
ся величина 8 Коэффициент берется из теории ламинарного пограничного слоя
на плоской пластине (см. [12]). Расчеты показывают, что этот способ определения толщины пограничного слоя обычно дает наиболее гладкое продольное распределение толщины пограничного слоя. В турбулентной области следовало бы использовать иное значение этого коэффициента, но для рассмотрения ЛТП это не существенно.
Однако 8т имеет нефизичное значение около точки торможения, где скорость внешнего
потока У^ —>0. Около такой точки используется иной способ. Пусть у®-05 — касательная скорость в ближайшей к стенке ячейке, где впервые выполняется условие
(
Жх.
дУі¥
\
П
<0.05
дУі¥
. Тогда в качестве альтернативной оценки толщины пограничного
слоя выбирается расстояние от стенки до точки, где Ух = 0.99УТ°'°5. Это расстояние обозна-
5-0.05 чим о0 99.
Близость К точке торможения определяется ПО величине р0 . ! IРц (р() — давление торможения). Если р0 . ! I [\} < 0.9, данный ряд ячеек расположен далеко от точки торможения, и в качестве окончательного значения толщины пограничного слоя выбирается 8 = 5т. Если
о
р р
-----— > 0.95. то точка торможения близко, принимается 5 = §099- Если ------— е 0.9; 0.95 , то
Ро ' Ро
о V? о0.05
значение 6 определяется линеинои интерполяциеи между От И О0 99.
Еще один важный эффект, который имеет прямое отношение к алгоритму инициирования перехода, был обнаружен в ходе расчетов обтекания профиля крыла. Этот эффект проявляется в глубине пограничного слоя, на очень малом расстоянии от стенки, и его неудобно иллюстрировать на рисунках из-за наклона твердой поверхности и большой толщины профиля по сравнению с размерами интересующей области течения. Однако подобное явление можно обнаружить и в расчетах пограничного слоя на пластине (см. рис. 8), где показано поле X (длины области ЛТП).
У самой твердой стенки можно заметить чрезмерный рост параметра X. Это связано с тем, что расчет начинается с однородного поля, где 9 = 0да <к: 9кр и Д = 0. Поэтому во всем поле течения параметр X рассчитывается по формуле (4), в которую подставляются местные характеристики пограничного слоя. Из-за этого в начале счета X достигает высоких значений в районе конца пластины. Когда достигается 9 = 0кр, то значения X начинают сноситься с критической поверхности вдоль линий тока, постепенно заменяя большие значения л|Сор.
Расчет проводился по явной схеме с локальным шагом по времени [10]. У твердой стенки размер ячеек сетки очень мал, и потому величина шага по времени также очень мала. Поэтому в районе стенки очень долго сохраняется область больших значений X (она показана стрелкой на рис. 8). Этот эффект не проявляет себя заметным образом в задаче о пластине, но в расчетах профиля крыла он может изменить значение X Дтах . Дело в том, что в районе конца профиля
крыла линии тока отклоняются от поверхности профиля, и информация из окрестности твердой стенки с прилипанием проявляется достаточно далеко от поверхности. Чтобы исключить этот численный эффект, была предложена следующая модификация инициатора перехода:
1. Если в данном ряде ячеек Дтах > §, то /.|Сор не вычисляется по формуле (4), а берется
равным значению X Дтах на предыдущем шаге по времени.
2. Атах и X Атах выбираются только среди ячеек с уп- >0.055.
3. Формула (12) для вычисления X модифицируется следующим образом:
~ дрХ
где Л г,? —решение уравнения-----------
а
р Хи1 -р
А <8, А >5, ( \ г_
А
у +
Рг
V
дХ
сЬс,-
(13)
= 0, в которое введена искусст-
венная «турбулентная диффузия». Число Прандтля в (13) выбирается равным Рг[ = РГ|Х =100. Слабая искусственная диффузия сглаживает поле X и приносит малые значения X от средней части пограничного слоя к стенке.
Следует отметить, что первая из перечисленных модификаций приводит к тому, что результат зависит от предыстории развития численного решения. Поэтому рекомендуется начинать расчет с использованием ступенчатой инициирующей функции:
Го, 0<екр,
/лтп 9>^лтп „ (14)
[1, 0 > 9кр.
Функция (14) соответствует нулевой длине области ЛТП. Она обеспечивает хорошее начальное приближение к решению. На конечных стадиях установления решения можно включать инициатор перехода с плавной инициирующей функцией и с перечисленными выше модификациями.
Описанные модификации позволили достичь нормальной работы инициатора перехода во всех рассмотренных тестовых задачах.
ЗАКЛЮЧЕНИЕ
Основным результатом описанной работы является построение устойчивого численного алгоритма постепенного «включения» полуэмпирической модели турбулентности в области ЛТП в пограничном слое на поверхности обтекаемого тела.
Авторы надеются, что описанный в данной статье поиск алгоритма инициирования перехода имеет методическое значение и может оказаться полезным при решении других задач. В частности, алгоритм определения толщины пограничного слоя может иметь более широкое применение. Логика описанных поисков, которая в конце концов привела к идее решения транспортных уравнений для расстояния от критической поверхности А и для длины области Л 111 X, позволяет понять, почему в полуэмпирической модели ЛТП Лэнгтри — Ментера [2] было предложено решать транспортное уравнение для такого, казалось бы, локального параметра, как Reg** rhp .
Что касается согласования расчетов ЛТП с экспериментальными данными, то оно в основном определяется точностью предсказания положения критической поверхности и длины области перехода X. Возможно, что формула Дэя и Нарасимхи (4) может дать неудовлетворительные результаты в более сложных задачах (в особенности при наличии отрывов, при влиянии шероховатости поверхности на ЛТП и т. д.). Для таких случаев, возможно, следует разработать специальную модель нелинейной стадии ЛТП.
Требуются дополнительные работы для обобщения описанного алгоритма на нестационарные задачи, на нерегулярные сетки. Следует также искать способы устранения зависимости алгоритма от предыстории развития решения.
Большая часть расчетов, описанных в настоящей статье, была выполнена в рамках европейского проекта TELFONA (6-я Рамочная программа). Продолжение работ поддержано Российским фондом фундаментальных исследований (грант 10-08-00274-а). Авторы также благодарны
Е. В. Кажану (НИО-1 ЦАГИ), который предложил идею использования параметра р0 ,_XJ
для идентификации точки торможения.
ЛИТЕРАТУРА
1. Singer B. A. Modeling the transition region // NASA CR-4492. 1993, p. 1 — 90.
2. Langtry R. B., Menter F. R. Correlation-based transition modeling for unstructured parallelized Computational Fluid Dynamics codes // AIAA J. 2009. V. 47, N 12, p. 2894—2906.
3. Menter F. R. Two-equation eddy-viscosity turbulence models for engineering applications // AIAA J. 1994. V. 32, N 8, p. 1598—1605.
4. Vlasenko V. V. New semi-empirical model for the prediction of laminar turbulent transition // Proceedings of 7th ONERA-TsAGI Seminar. 2009.
5. Власенко В. В., Морозов А. Н. Новая полуэмпирическая модель для предсказания ламинарно-турбулентного перехода // Материалы XX школы-семинара «Аэродинамика летательных аппаратов». 2009, с. 41 — 42.
6. Жигулев В. Н., Тумин А. М. Возникновение турбулентности. — Новосибирск:
Наука, 1987, с. 1 —282.
7. D e y J., N a r a s i m h a R. Integral method for the calculation of incompressible two dimensional transitional boundary layers // J. Aircraft. 1990. V. 27, N 10, p. 859—865.
8. Schubauer G. B., Skramstad H. K. Laminar boundary layer oscillations and stability of laminar flow // JAS. 1947, V. 14, p. 69—78.
9. Шлихтинг Г. Теория пограничного слоя. — М.: Наука, 1974, с. 1 — 711.
10. Власенко В. В. О математическом подходе и принципах построения численных методологий для пакета прикладных программ EWT ЦАГИ. — Сб.: «Практические аспекты решения задач внешней аэродинамики двигателей летательных аппаратов в рамках осред-ненных по времени уравнений Навье — Стокса» // Труды ЦАГИ. 2007, вып. 2671, с. 20—85.
11. C o a k l e y T., H s i e h T. Comparison between implicit and hybrid methods for the calculation of steady and unsteady inlet flows // AIAA-85-1123. 1985.
12. Абрамович Г. Н. Прикладная газовая динамика. — М.: Наука, 1991. Т. 1, с. 1 —600.
Рукопись поступила 24/IV 2010 г.
