Численное моделирование ламинарных и турбулентных течений в следе за телом Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 2, № 6, 1997

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЛАМИНАРНЫХ И ТУРБУЛЕНТНЫХ ТЕЧЕНИЙ В СЛЕДЕ ЗА ТЕЛОМ *

В. М. Ковеня, А. С. Лебедев Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: kovenya@adm.ict.nsc.ru

The results of numerical modeling of the flows in the near wake behind plane and axisymmetric bodies in a wide range of Mach and Reynolds numbers are presented. The problem solution was sought in the approximation of Favre-averaged Navier — Stokes equations for a compressible heat-conductive gas complemented by a semi-empirical turbulence model, which were solved with the method of splitting into physical processes and space directions. The peculiarities of the laminar, transitional and turbulent flows have been studied.

Исследование течения в следе за движущимся телом является одной из важнейших задач аэродинамики, так как оно оказывает существенное влияние на общую картину течения и, как следствие, на аэродинамические характеристики летательного аппарата и его элементов. Структура течения в следе сложна и многообразна и зависит от многих факторов, что затрудняет ее изучение экспериментальными и численными методами (см., например, [1-3]). При экспериментальном исследовании этого класса задач, как правило, удается получать распределение параметров по поверхности обтекаемого тела, а саму структуру течения в следе лишь в отдельных сечениях. Математическое исследование отрывных течений в приближении упрощенных моделей, а также интегральными и полуэм-пирическими методами не является достаточно точным и также не позволяет получить детальную картину течения в следе. По существующим предположениям в наиболее полной постановке отрывные течения описываются осредненными уравнениями Навье — Стокса сжимаемого теплопроводного газа, дополненными той или иной моделью турбулентности. Однако результаты численных экспериментов показывают, что сам выбор модели турбулентности является проблемой и, следовательно, задача изучения отрывного течения в следе становится еще более сложной и может быть решена с достаточной точностью лишь в некотором диапазоне параметров потока. Таким образом, численное моделирование служит двум целям — получению новых знаний об особенностях течения (в частности, течения в следе) и тестированию выбранной модели турбулентности с целью определения диапазона ее применимости. В настоящей работе представлены результаты численного моделирования стационарных ламинарных и турбулентных отрывных течений в следе за плоскими

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №96-01-00136.

© В. М. Ковеня, А. С. Лебедев, 1997.

и осесимметричными телами. В качестве математической модели выбраны осредненные по Фавру уравнения Навье — Стокса сжимаемого теплопроводного газа, дополненные полуэм-пирической моделью турбулентности типа д — ш (см. [4]). Для численного решения задачи исходные уравнения Навье — Стокса аппроксимировались разностной схемой расщепления по физическим процессам и пространственным направлениям [5]. Дифференциальные уравнения, описывающие д — ш модель, решались по схеме вида [5] с расщеплением операторов по пространственным переменным. Стационарное решение задачи находилось как предельное решение нестационарных уравнений со стационарными краевыми условиями. Исследованы течения в следе за конусом и торцом пластины в широком диапазоне чисел Маха и Рейнольдса, включая ламинарный, турбулентный и переходный режимы. Получены основные закономерности течений, и выявлены некоторые их особенности в определенном диапазоне параметров потока. В работе авторов [6] было обнаружено, что при сверхзвуковом обтекании осесимметричных тел в предположении ламинарности потока локальная скорость течения в рециркуляционной зоне становится сверхзвуковой Млок ^ 1.05. Эти режимы течения имеют место при диапазоне чисел И,е близких к переходному режиму. Подобный эффект был получен также в работах [7, 8]. При исследовании течений за торцом пластины в области переходного режима имеют место два стационарных решения, аналогичных явлению “гистерезиса” по числу Рейнольдса. Так, длина рециркуляционной зоны в некотором диапазоне чисел И,е различна при различных законах изменения чисел И,е. Последняя серия расчетов посвящена модификации д — ш модели турбулентности путем изменения демпфирующей функции. Расчеты течения около конуса с углом полу-раствора 7° при Мте ^ 9.16 и сравнение с имеющимся экспериментом [9] позволили сделать заключение, что модифицированная модель дает хорошее совпадение с экспериментальными данными.

1. Постановка задачи и численный алгоритм

Исследуется стационарное течение вязкого сжимаемого теплопроводного газа в ближнем следе за телом. Типичная картина течения около обтекаемого тела представлена на рис. 1. Течение предполагается плоским или осесимметричным. Расчетная область начинается на некотором расстоянии вверх по потоку от донного среза и включает зону ближнего следа (см. рис. 1, область ЛБОПЕЕ). На поверхности тела задаются условия прилипания для скорости и режим тепловой изоляции или охлаждения для температуры. Прямая ЕП является линией симметрии. На входной границе задается развитый пограничный слой, полученный либо из решения задачи об обтекании головной части тела, либо, как при обтекании торца пластины, из приближенного решения. Внешняя граница БО либо является головным скачком уплотнения, либо выбирается достаточно далеко от тела с постановкой на ней условий невозмущенного потока. На выходной границе задаются “мягкие” условия.

В качестве математической модели для описания течения в рассматриваемой области используются осредненные по Фавру уравнения Навье — Стокса, которые в обезразмерен-ном виде в декартовых координатах могут быть записаны в виде

ди+^=<*• р=(1.1)

Рис. 1.

и =

3 = 1, 2.

где

/ Р \ / Риз \

ри ; . = ри и. + а.

ри2 ’ . ри2 и. + ^2.

V рЕ / \ рЕи. + и 6. — qj )

Здесь и ниже по повторным индексам предполагается суммирование, р — плотность; и1, и2 — компоненты скорости; Т — температура; Е = Т + и.и./2 — полная удельная энергия; р = (7 — 1)рТ — давление; 7 = ср/с^ = 1.4. Тензор напряжений а. и тепловой поток qj представляются в виде

ач =

6ц(р + 2 кр) — ( £ + №)

' 2 Л Зик

8x3 3 3 дхк

Яз = 7

р + ^ ЯеРг Рг^

дТ дх3 ’

где віз = дих/дхз + диз/дхх; — турбулентная вязкость; к — турбулентная кинетическая

энергия; Рг= 0.72; Рг^ = 0, 9, р = р(Т) — коэффициент динамической вязкости.

Турбулентная вязкость определяется через масштаб длины I и скорость турбулентности я = -/к: Р = Рі(Я, I). Величины I и я находятся из дифференциальных уравнений

д д д — (р£п) + — (рЯ из) =

д£

дх3

дх3

р +

Яе + Рг,

рг \ д$г

дх3

з

(1.2)

Здесь Sn = (1, q) (п = 1, 2); Нп = Нп(р, 51, $2, ди/дх.) — источниковые члены; Рг1 = 1;

Рг2 = 1, 3. Для замыкания уравнений (1.1) в данной работе выбрана полуэмпирическая модель [4] с двумя дифференциальными уравнениями для S1 = q = у/к и S2 = и = q/1. Источниковые члены Нп в (1.2) имеют вид

Н1 = 2 (сд!^я/и — 3— Р? Н2

Р?

где

О

дик; дхк ’

т___ 2 \ дих

7 = — 3 3 дхз;

РЯ

сі = 0.405# + 0.045; С2 = 0.92;

сд = 0, 09; # = 1 — ехр(—вЯт); Ят = Яе—; в = 0, 0018; р£ = сД#ря1. На поверхности тела

р^

ЛЕЕ заданы условия

дТ

иі = и2 = 0, — =0, (1.3)

дп

где п — внешняя нормаль к поверхности, а на выходной границе СД “мягкие” условия

д/

= 0, / =(р,Иі,И2,Т)

т

(1.4)

Внешняя граница ВС выбирается в области невозмущенного потока, и на ней задаются условия

и1 1; и2 р рж 1; Т ТСХ> /_ 1 \ л * (1*5)

7(7 - 1)м^

2. Численный алгоритм

Численный алгоритм решения уравнений Навье — Стокса (1.1) основан на неявной схеме расщепления по физическим процессам и пространственным направлениям [5]. Для построения алгоритма представим уравнения (1.1) в недивергентном виде

3=1

5х3

к

(2.1)

где матричные операторы В, содержат газодинамические члены уравнений и вязкие члены (без смешанных производных) по направлению Xj, а вектор правых частей ^ включает смешанные производные в уравнениях движения и диссипативные члены в уравнении энергии. Представим одномерные операторы Bj• в форме расщепления по физическим процессам

д

Bj дХ" = С + С?'+2’ (2*2)

где С, содержит конвективные и вязкие члены по направлению х, , а — члены с

давлением в уравнениях движения и члены вида ^уи в уравнениях неразрывности и энергии. С учетом (2.2) исходные уравнения Навье — Стокса могут быть представлены в форме

3 = 1

(2.3)

где, например, для вектора / = (р,и1 , и2,Т)т матричные операторы С, имеют вид (^ 1, 2),

С,

щ 1 - РдХ"Ь

О'.Х'і

_д_

5ж,

С,-

3+2

/ 0 Р ^2, р 0 \

2 0 0 Ь2 5

2 2 0 0 ^ Ь2 5х3

V 0 6ц С2 2 С 2 0 /

1 др 1 др

Здесь I — единичная матрица, Ь, — диагональные матрицы, а2 = — —-, Ь2 = — ——, с2

р др р дТ

б, — символ Кронекера. Очевидно, (1.1) и (2.3) связаны соотношениями

(2.4)

Р 7“, Р

/ = - £ с, / + *

3=1

-А

-і

(2.5)

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

£ "+1 £

£ £ I \ А /~'Ц ГЛ, Г"+1 | /1 Л,\ £ПЛ __ Т?П

т

+ £ С3>/ ”+‘ + (1 - а)/п] = к,

аппроксимирующую уравнения (2.3) с порядком 0(т + к1). Здесь С3., — разностная

аппроксимация дифференциальных операторов С и ^ с порядком 0(к-). Приближенно факторизуя оператор I + та ^4= С". по формуле

4 4 4

I + ™ Е 3 ю П(I + Т“СЙ) = / + та £ Сд + 0(т2),

3 = 1 3 = 1 3=1

рассмотрим разностную схему

4 /п+1 /га 4

П(1 + таС'„)/------------^ = - £ С/п + ¥,

3=1 3=1

или, с учетом соотношений (2.5), разностную схему

4 рп+Л пп 2

I —1\п \ А л кт

£ "+1 £ "

Д(1 + таС-'Л)7 т 7 = -(А-1)" £ Л‘И<7*. (2.6)

3=1 Т 1=1

аппроксимирующую уравнения Навье — Стокса (1.1) с порядком 0(т+тк1 + кк), а при установлении — стационарные уравнения Навье — Стокса в консервативной форме с порядком 0(кк). Здесь Л* = 5/ <9ж3- + О(к^), к = тах | к, |, 1 < ^ < 2. В линейном приближении

разностная схема (2.6) безусловно устойчива при а > 0.5 для I = к и при а =1 для I < к.

Для численной реализации схемы (2.6) рассмотрим эквивалентную ей разностную схему в дробных шагах:

еп = -£ л и,к

1=1

(I + таСУ Єп+1/4 =

2

(I + таС-у £"+■ = е"+3/4

/"+1 = /" + т^"+1. (2.7)

В соответствии с расщеплением операторов по физическим процессам и пространственным переменным разностная схема (2.7) на каждом дробном шаге реализуется скалярными прогонками. Для упрощения реализации операторы С3- аппроксимируются с первым порядком (/ = 1), а операторы правой части со вторым (к = 2), что позволяет свести решение разностных уравнений к скалярным трехточечным прогонкам. Уравнения (1.2) для определения турбулентных параметров решались независимо от (1.1) после каждого шага по времени т с помощью аналогичной (2.7) разностной схемы

Сга+1 сга 2

(р" + таЯ1,,,) (р-1)" (р" + таЯ2„)--------;-----= - £ Я”. + Н (2.8)

3 = 1

или эквивалентной ей схемы в дробных шагах

е = - Е 3 +

3=1

(рп + тад;,) еп+1/2 = еп, (рп + тад2,) еп+1 = рпеп+1/2,

£п+1 = £п + теп.

(2.9)

д

Здесь Я3. = Л3 р"и" — Л3- Л3 — аппроксимация дифференциальных операторов —— ри3- —

0X3

д д

——^7 ——, с порядком 0(к3-), Н. — аппроксимация источниковых членов Нп. Подобно (2.7) 0x3 0x3 ^

разностная схема (2.9) реализуется скалярными прогонками и безусловно устойчива. Типичное число расчетных узлов в продольном и поперечном направлении изменялось от 100 до 200 по каждому направлению. Разностная сетка сгущалась в окрестности пограничного слоя и области ближнего следа. Как правило, при нахождении решения в качестве начального приближения задавались параметры потока, полученные при решении задачи при близких параметрах (числах Маха и Рейнольдса).

2

3. Результаты расчетов

Большинство расчетов, в которых течение в ближнем следе предполагалось ламинарным, выполнены для короткого затупленного конуса с углом полураствора 25°. Длина обтекаемого тела вдоль оси симметрии равнялась двум радиусам затупления. Числа Маха (Мте ) и Рейнольдса (И,ете ), подсчитанные по параметрам набегающего потока и радиусу затупления, изменялись в диапазонах: 2 < Мте < 6, 103 < И.ете < 105.

Проведенные расчеты, позволив для каждого набора определяющих параметров получить детальную картину течения и воспроизведя, например, такие известные из эксперимента и соответствующие теоретическим представлениям тенденции, как приближение точки отрыва к угловой точке и увеличение длины отрывной зоны с ростом числа Рейнольдса или уменьшением числа Маха, выявили также и ряд интересных особенностей, среди которых выделяются формирование в окрестности точки отрыва мелкомасштабных вихрей и возникновение локальной сверхзвуковой зоны в области возвратного течения. Следует отметить, что изменение структуры течения в окрестности точки отрыва, заключающееся в формировании системы локальных вихрей, обнаруживалось в расчетах всякий раз при изменениях определяющих параметров, приводящих к смещению точки отрыва к угловой точке, то есть при увеличении числа Рейнольдса или уменьшении числа Маха, независимо от того осесимметричное рассматривалось течение или плоское. Иначе обстояли дела с возникновением сверхзвуковой (Мте ~ 1.05 — 1.08) зоны в возвратном течении, положение и размеры которой см. на рис. 2, где приведены изолинии числа Маха для параметров набегающего потока Мте = 2.5, И.ете = 5000. И это кажется особенно уместным подчеркнуть в связи с тем, что, в отличие от довольно обширного экспериментального материала, подкрепляющего возможность изменения картины течения возле точки отрыва, экспериментальные данные, свидетельствующие о существовании при некоторых режимах обтекания локальной сверхзвуковой зоны в рециркуляционной области ближнего

1.0 0.8 0.6 0.4 0.2

Рис. 2. Рис. 3.

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

При исследовании стационарного турбулентного течения в окрестности задней вертикальной кромки пластины конечной толщины числа Маха и Рейнольдса, вычисленные по значениям газодинамических параметров в однородном изэнтропическом потоке над пластиной и полутолщине пластины К, принимали значения: М = 2.8; 4.0 и 6.0; Ие = 2.5 • 103 — 106 . Задаваемые профили величин во входном сечении АВ соответствовали развитому турбулентному пограничному слою, причем процедура построения этих профилей была такова, что в рассматриваемом диапазоне определяющих параметров рост числа Ие приводил к увеличению турбулентой вязкости ^. Влияние числа Рейнольдса на длину отрывной зоны х0/К показано на рис. 3.

Обратим внимание на то, что при М = 2.8 для Ие приблизительно из интервала от 105 до 4 • 105 на графике показаны два значения длины отрывной зоны, это соответствует существованию двух различных стационарных решений задачи при совпадающих определяющих параметрах. Отличие расчетов состояло лишь в выборе начального поля: верхняя ветвь графика получена при выборе в качестве начального поля уже найденного стационарного решения для небольшого значения Ие, а результаты, соответствующие нижней ветви графика, получены, когда за начальное поле принималось стационарное решение для большого значения Ие. Имеет место, таким образом, “гистерезисное” по числу Рейнольдса явление. Для М = 4.0, как видно из рис. 3, тоже существует интервал двойственности

0,3 ....... I I I I ■ 11 I ■ I ■ ■ -...г11 | I I | |

О 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 О 0,1 0,2 0,3 0 0,1 0,2 Х/к

Рис. 4.

решения, и при больших Ие происходит резкое сокращение длины отрывной зоны. Следует подчеркнуть подтвержденное многочисленными расчетами обстоятельство: какие бы данные ни выбирались в качестве начальных, получавшееся стационарное течение всегда соответствовало изображенным на рис. 3 кривым и никаких других решений не возникало. Направление скорости в окрестности угловой точки для М = 2.8 показано на рис. 4. Появление пары вращающихся в противоположных направлениях мелкомасштабных вихрей в окрестности точки отрыва на вертикальной кромке пластины происходит при увеличении Ие аналогично тому, как это имело место при расчетах ламинарного течения, однако эти вихри исчезают, когда Ие достигает значения 6 • 105. Точка отрыва основного потока на вертикальной кромке с увеличением Ие монотонно приближается к угловой точке, а затем несколько удаляется от нее.

Ре=300000

Рис. 5.

На рис. 5 при М = 2.8 представлены профили турбулентной вязкости в нескольких поперечных сечениях. Для Ие = 3 • 105 приведено то их двух решений, которое соответствует “ламинарному” режиму с протяженной отрывной зоной, а для Ие = 6 • 105 — единствен-

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

В заключение приведем результаты расчетов, о качестве которых можно судить путем сравнения с имеющимися экспериментальными данными, описанными в [9]. Обтекаемое тело представляет собой острый конус длиной 0.578 т с углом полураствора 7°, закрепленный на цилиндрической державке диаметром 0.102 т. Параметры набегающего потока:

= 59.8 К, Ието/т = 5.5 • 107, Мто = 9.16. Температура поверхности конуса поддерживалась равной 290 К.

Профили величин во входном сечении АВ получались следующим образом. Сначала рассчитывалось течение возле боковой поверхности конуса в предположении ламинарно-сти потока возле носка. На расстоянии 5-10 мм вниз по течению от носка турбулентным параметрам д и I приписывались некоторые, вообще говоря, произвольные значения, подчиненные, однако, условию, выполнение которого достигалось методом проб и ошибок. Целью было внести уровень турбулизации потока возле носка, который достаточно высок для того, чтобы течение далее вниз по потоку развилось в турбулентное, не ламинари-зовалось, но в то же время как можно более низок, чтобы получающееся турбулентное течение не зависело от произвольно заданных д и I возле носка. Полученное таким образом течение над боковой поверхностью конуса очень хорошо согласуется с найденными экспериментально профилями числа Маха и давления в нескольких нормальных к конусу сечениях, а величины на АВ (см. рис. 1) — с корректно заданными для последующего расчета и сравнения в области ближнего следа.

Рис. 6.

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

Pt/Pinf

Рис. 7.

плавно отключить действие этой функции вдали от твердой поверхности (не меняя ее поведения вблизи) и предотвратить тем самым избыточную ламинаризацию потока в волне разрежения, мера удаленности от поверхности при этом оценивалась с помощью естественного поведения демпфирующей функции в сечении, расположенном несколько выше по потоку от донного среза. На рис. 6 показаны полученные в расчете изолинии плотности рядом с фотографией течения, взятой из [9], на которой различимы пограничный слой у боковой поверхности конуса, начало волны разрежения и хвостовой скачок. Длина отрывной зоны согласно эксперименту составляет 40.3 мм, расчет дал очень близкое значение 39 мм. Вычисленные и экспериментальные профили давления в нескольких нормальных к поверхности державки сечениях показаны на рис. 7. Приведенные сравнения свидетельствуют о хорошем совпадении результатов расчета с экспериментальными данными, особенно для рассматриваемого класса течений.

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

[1] Чжен П. Отрывные течения. Мир, М., 1972-1973, Т. 1-3.

[2] Швец А. И., Швец И. Т. Газодинамика ближнего следа. Наукова думка, Киев, 1976.

[3] Ковеня В. М., Лебедев А. С. Численное моделирование отрывного турбулентного течения в ближнем следе за пластиной. ПМТФ, 37, 1996, 106-113.

[4] CoAKLEY T. J. A Compressible Navier —Stokes Code for turbulent flow modelling. N.Y., 1985 (NASA TM 85-899).

[5] Ковеня В. М., Яненко Н. Н. Метод расщепления в задачах газовой динамики. Наука, Новосибирск, 1981.

[6] Ковеня В. М., Лебедев А. С. Численное моделирование отрывного турбулентного течения в ближнем следе за пластиной. ПМТФ, 30, №5, 1989, 53-59.

[7] Hollanders H. Devezeaux de Lavergue D. High speed laminar near wake flow calculation by implicit Navier —Stokes solver. In “AIAA 8th CFD Conf.”, 1987, 598-607.

[8] Егоров И. В., Зайцев О. Л. Об одном подходе к численному решению двумерных уравнений Навье — Стокса методом сквозного счета. Журн. вычисл. матем. и матем. физики., 31, №2, 1991, 286-299.

[9] Denman P. A., Harvey J. K., Hiller R. Hypersonic Boundary Layer and Base Flow. In “Proc. Workshop on Hypersonic Flow for Reentry Problems”, Antibes, France, 1990, 2, 43-62.

Поступила в редакцию 16 декабря 1997 г.