Фиг. 1.
Фиг. 2.
Конфигурация расчетной области близка к использовавшейся в [2], однако там ось симметрии и выходная граница задавались на одной координатной линии, что не слишком удобно. В данном случае при построении сетки применялся параболический генератор в варианте из [8]. Используемый подход является альтернативой блочному. В ряде случаев его применение может оказаться предпочтительным.
Пример разностной сетки, использованной в расчетах, приведен на фиг. 2. Отдельно показан фрагмент сетки в окрестности задней кромки сопла. В качестве характерного размера здесь принимался радиус критического сечения сопла, радиус затупления задней кромки составлял от него 0.5%. Сетка имела необходимое сгущение к поверхности обтекаемого тела для расчета турбулентного пограничного слоя.
Для решения задачи использовалась цилиндрическая система координат, ось х которой совпадала с осью сопла. Решение определялось путем численного интегрирования нестационарных двумерных уравнений Навье-Стокса осредненного турбулентного течения. Обезразмеренные
традиционным способом по параметрам набегающего потока на бесконечности и характерному линейному размеру, они имеют следующий вид в цилиндрической системе координат (х, у):
ЭШ + ЭТ + ЭС + Су = яе-1 № + дд^
Э t Эх ду у V Эх Эу
у
р
и - ри , г -
рv
е
р и
2 ,
ри + р + р(
р и V (е + р + р() и
С =
рv
ри V
р v 2 + р+р,
(е + р + р() V
С у =
рv
р и V
2
рv
(е + р + рг) V
0 0 0
г - А V О хх Т ух 5 Сvy - Тух
Т ху Оуу 2ц( Vу - Vу"1)
Г х Гу Гу
где р1 = 2 р^ Охх = ц(2 их - 2 Уи) , тху = тух = ц(и + Vx), Оуу = ц^2 Vy - 2 V ^,
гх = иОхх + УГху + ц Рг1 Их, Гу = иТух + УОуу + ц Рг-1 Ну,
р = (У - 1)Т-1рИ, е = р[у-1И + 0.5(и2 + V2) + к], Уи = их + Vy + vy-\
ц = ц + ц, ц Рг1 = ц Рг-1 + ц Рг-1, Яе = р^и^Н ц-1.
Здесь , - время, р - плотность, и и V - компоненты вектора скорости в направлениях х и у, р -давление, И - энтальпия, е - полная энергия, к - кинетическая энергия турбулентности, у - отношение удельных теплоемкостей, ц и ц, - коэффициенты молекулярной и вихревой вязкостей, Рг; = 0.72 и Рг, = 0.9 - ламинарное и турбулентное число Прандтля, Яе - число Рейнольдса. Система дополнялась зависимостью молекулярной вязкости от энтальпии, в данном случае - в виде формулы Сатерланда (см. [9]). В том случае, когда газ струи имеет отличный от внешнего потока показатель адиабаты, необходимо использовать также уравнение концентрации.
Для замыкания системы уравнений использовалась дифференциальная двухпараметрическая (д-у)-модель (см. [10]). Параметрами модели являются q = к1/2 - скорость турбулентности и V = = ^¿^У'2, где Ь - масштаб турбулентности. Ее уравнения имеют вид
Эи, ЭЭС, С, „ -1 ЭС^ С^Ч тт ' + - + - + -! - Яе 1 ( —^ + —£ + I + И,
о, Эх Эу у V Эх Эу у
Л =
рq 1ру]
Г =
рuq р uv
С, =
р Щ
р VV
Г, V =
ЦqPГq1qx 5 СtV - qy , И, - Htq
ц Рг-1 v х Vy_
где И, - вектор источниковых членов, компоненты которого имеют вид
Нщ = Ш (^-2 ^ - V2
Н V - ~
а CV#-3 -
-1 2 2 2 —1 2 ^2 D = ux + vy + vy , S = 2ux + uy + vx + 2( Vy + vy ) - 3D ,
MqPr-1 = Ml + Pr-1 Mt, Mv PrV1 = Ml + PrV1 M,
a = a/ + «2, / = 1 - exp (MEB^T),
где n - расстояние по нормали к поверхности, а CM = 0.09, a1 = 0.56, a2 = 0.065, Р = 0.92, РМ = 0.0075, Prq = 1 и PrV = 1.3 - константы. Коэффициент турбулентной вязкости определяется по формуле
М = Re C/Mp(q/V)2.
На границе BD задавались условия прилипания для компонент скорости и температура торможения, разная для внутреннего и наружного участков поверхности. Значение плотности определялось или путем решения уравнения неразрывности, или из положения о равенстве нулю производной от давления по нормали к стенке. В качестве граничных условий использовались также тождество q = 0 и аппроксимация односторонними разностями уравнения для V в узлах сетки, совпадающих с положением твердой поверхности.
На линии AH задавались условия симметрии течения. На входной-выходной границе DH фиксировались внешнее давление, остальные параметры рассчитывались на основе характеристических соотношений и предположения о сохранении на границе вектора скорости. Параметры на линии AB определялись исходя из полных давления и температуры газа в форкамере, а также из уравнения для характеристики, направленной против потока. Вертикальная компонента вектора скорости полагалась равной нулю.
Расчеты проводились после перехода к произвольной криволинейной системе (см. [11]) в квадратной области с координатами ^ и п (см. фиг. 1). Линии ^ = const подходили перпендикулярно как к твердой стенке, так и к оси симметрии. В свою очередь, линии п = const в зоне пограничного слоя были параллельны обтекаемой поверхности.
2. РАЗНОСТНАЯ СХЕМА И МЕТОД РЕШЕНИЯ
Численное интегрирование уравнений проводилось методом установления с использованием разностной схемы пятого порядка (см. [12]). На сетке wx = {х, = ih, h = const} оператор для первой
производной /х выглядит следующим образом:
6 E + ^)Ао- 5 Л2
+ 5 [ E - Л52]Л2 /,
где А/ = / + 1 - / _ 1, А/ = / + 1 - 2/, + / _ 1, а 5 - параметр, определяющий его ориентацию; его значение в расчетах составляло 2.5. Первая компонента данной формулы является известной симметричной разностью шестого порядка. Вторая представляет собой монотонизирующий оператор для устранения нефизичных высокочастотных осцилляций. Неявные части одинаковы и записаны на трех точках, что позволяет использовать при решении метод прогонки. Данная разностная формула использовалась для представления гиперболических членов уравнений. Симметричная ее часть применялась для описания смешанных производных с вязкостью и вычисления метрических коэффициентов преобразования системы координат.
Для описания вторых производных с вязкостью вида (м /х )Х использовался оператор
(М/г) г = ¡Ц [ + 1/2/г+1/2- М-г - И2/1 - 1/2 ) - Мг + 1/г+1 + Мг-\/г-\]
Мг ± 1/2 = 1 [9(Мг + Мг ± 1) - Мг- 1 + Мг + 2 К
/ ± 12 = 215 [27(±/ ± 1) ± /г + 1 - /г ± 2 ],
f± i = ^ [ Ff + 2 ± 6/ +1 + 18/, ± 10/±! ± 3/±2].
Полученная в результате разностной аппроксимации уравнений Навье-Стокса и модели турбулентности система линейных алгебраических уравнений решалась итерационным методом Гаусса по линиям ^ = const. Для интегрирования по времени использовалась схема второго порядка. Турбулентным течениям с отрывами пограничного слоя, вдувом струй и т.д. свойственна нестационарность. Решение считалось сошедшимся, когда после формирования течения и снижения невязки разностных уравнений на 3-4 порядка по сравнению со своим максимальным значением она совершала относительно небольшие периодические колебания.
3. РЕЗУЛЬТАТЫ РАСЧЕТОВ 3.1. Расчеты обтекания кормовой части осесимметричного тела внешним потоком
Рассматривалось обтекание хвостовых частей, исследованных экспериментально в работах [13], [14] и численно в [2]. Число Рейнольдса, посчитанное по параметрам набегающего потока и характерному размеру, составляло от 107 до 1.2 х 107. Значение кинетической энергии к как в набегающем потоке, так и на входе в сопло равнялось 10-6, что давало отношение турбулентной и ламинарной вязкостей ~10-2. Температура торможения струи совпадала с температурой торможения набегающего потока. Нулевое значение координаты х соответствовало началу сужения внешней поверхности тела. Во всех расчетах различных конфигураций кормовых частей с соплами использовалась сетка с 141 х 81 узлами. Минимальное расстояние между узлами равнялось 10-5.
На левой границе внешнего течения задавался турбулентный пограничный слой толщиной 0.01. Профиль скорости внутри него соответствовал закону "1/7", а распределение температуры - интегралу Крокко (см. [9]). Ниже по течению пограничный слой утолщался, в нем формировались профили турбулентных параметров. В процессе расчета в районе сужения кормовой части, как правило, возникала область сверхзвукового течения с замыкающим ее скачком уплотнения, а за срезом сопла развивался слой смешения внешнего потока и вытекающей струи.
Обтекание хвостовой части тела, исследованного в работе [13] и имеющего сужающуюся часть, описываемую дугой окружности, рассчитывалось при значении числа Маха 0.8. Отношение полного давления струи к давлению в набегающем потоке = 5. Помимо этого рассматривались варианты обтекания кормовой части тела с жестким Цилиндрическим имитатором, а также при отсутствии струи. В последнем случае расчет проводился в области, одной из границ которой являлась координатная линия, проходящая через критическое сечение сопла. На ней задавались условия, соответствующие твердой поверхности.
На фиг. 3 представлены полученные (линии 1) распределения коэффициента давления при вдуве струи (график (а)), ее отсутствии (график (б)) и жестком цилиндрическом имитаторе (график (в)). Там же приведены экспериментальные данные из [13] (линии 2). Видно, что наличие газовой струи или ее жесткого имитатора ведет к повышению давления в донной области. Резкое снижение значения Ср на фиг. 3а приходится на заднюю кромку сопла, за которой формируется донное течение.
Фиг. 3.
Поле линий M = const представлено на фиг. 4. Минимальное значение изолинии 0.055, шаг 0.05. Внешнее течение полностью дозвуковое. Перед струей возникает зона отрыва потока. Сопло в данном случае не имеет сверхзвукового участка. Однако давление в выходном сечении выше, чем во внешней среде, поэтому за кромкой сопла струя продолжает расширяться и становится сверхзвуковой. Видны скачки уплотнения в струе и сдвиговый слой.
Другая форма кормовой части рассматривалась при значении числа М» = 0.9. Она имела коническую сужающуюся часть с углом наклона 7.9° и коническое же сопло с углом полураствора 11.77°. За характерный размер принимался радиус критического сечения сопла rc. Отношение полного давления струи к давлению в набегающем потоке p/p» = 7.82. Данный случай является более сложным для расчета, поскольку на боковой поверхности образуется зона сверхзвукового течения, замыкающаяся скачком уплотнения. Происходит взаимодействие скачка с турбулентным пограничным слоем.
Полученное распределение коэффициента давления по внешнему контуру тела при вдуве струи (графи
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.