Расчет обтекания затупленных несимметричных конусов с большим углом раствора Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Т о м IV 197 3

№ 5

УДК 533.6.011.3/55:629.7.024.36

РАСЧЕТ ОБТЕКАНИЯ ЗАТУПЛЕННЫХ НЕСИММЕТРИЧНЫХ КОНУСОВ С БОЛЬШИМ УГЛОМ РАСТВОРА

А. П. Базжин, С. В. Пирогова

Конечноразностным методом, основанным на принципе установления, рассчитаны течения перед тупыми осесимметричными телами и пространственное течение перед телом, не обладающим осевой симметрией. Рассмотренные тела получены из конуса с углом раствора 120° путем сферического затупления носовой части и скруглення боковой кромки. Рассчитаны варианты течения, отличающиеся значением числа М сверхзвукового набегающего потока, геометрией тела, свойствами газа и величиной угла атаки. Приведены отдельные результаты расчетов.

1. Задачи о расчете смешанных течений газа перед тупыми телами неоднократно рассматривались в литературе. Двумерные смешанные течения рассчитывались многими авторами достаточно точно. Подробную библиографию и полное освещение вопроса можно найти в монографии [1]. Задача о трехмерных смешанных течениях тоже была решена различными численными методами, в том числе методами, основанными на принципе установления (см., например, [1, 2]).

Однако создание первых алгоритмов расчета пространственных смешанных течений газа означало лишь начало обширных численных исследований в этой области. Большое разнообразие таких течений, обусловленное различием в геометрии исследуемых тел и режимах полета, вызвало появление других работ, в которых предлагались иные алгоритмы, предназначенные, как правило, для расчета течений какого-либо конкретного вида. Например, в недавно опубликованной работе [3], рассматриваются смешанные течения газа, не обладающие симметрией.

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

Дано краткое описание численного алгоритма и приведены отдельные результаты расчетов осесимметричных тел рассматриваемого типа при нулевом угле атаки и одного несимметричного тела при углах атаки от —-10° до +10°.

2. Течение рассматривается в сферической системе координат гФв с подвижным центром (фиг. 1, а). Положение центра определяется величиной отсчитываемой вдоль оси Ф = 0 от

некоторой точки О. Для решения задачи используется конечноразностный метод, основанный на принципе установления. Нужно

рассчитать предельное состояние (при ^ оо) течения в зоне влияния, ограниченной поверхностями тела [г= в (Ф, 0,Щ, головной ударной волны [г = (Ф, 0, /)] и некоторой замыкающей по-

верхностью пространственного типа. Газ считается совершенным или находящимся в состоянии термодинамического равновесия.

Разностная схема и решение строятся на равномерной сетке в пространстве нормированных независимых переменных га, <р, 0, Т. Переход к этим переменным производится по формулам:

„ = р fa г г — а (Ф. е, t) .

П Г J уп), п w (ф> Q t)_Q (ф> 0) () -

Т = ^2(Ф,в); =^з(в); T=t.

(1)

Выбором функций F3, F2, F: задается желаемое расположение меридиональных плоскостей (0 = const), лучей (0=const, ® = const) и узловых счетных точек на луче в физическом пространстве гФ0.

Ю

В пространстве га<р& расчетная сетка остается всегда равномерной по всем направлениям. В проведенных расчетах были использованы функции Fl(ti)=n, /73(0)==0, а в качестве функции сгущения лучей р\(Ф, 0) было использовано несколько функций, не зависящих от 0. Вместе с выбором функции 50(®, определяющей положение центра сферической системы координат, это дает возможность расположить счетные лучи нужным образом в каждой меридиональной плоскости, т. е. так, чтобы они сгущались в тех областях течения, где можно ожидать резкого изменения газодинамических величин. Этот прием был ранее предложен в работе [1].

Система газодинамических уравнений в нормированных независимых переменных была использована в виде

df df df df w + Mirn + Ni>t + L-m + H = 0’

где компонентами вектор-функции / являются составляющие скорости и, v, w, логарифм давления Р—\пр и логарифм плотности /? = 1пр. Все газодинамические параметры соответствующим образом приведены к безразмерной форме. Матрицы М, N, L размером 5X5 зависят от компонентов вектор-функции /, от координат г и Ф и, кроме того, от функций С(Ф, 0, t) и ХР(Ф, 0, t), представляющих поверхности тела и ударной волны, их производных по Ф, 0, t, а также от £оф и £О0 — производных функции, определяющей движение центра системы координат. Эти матрицы не выписаны здесь в явном виде из-за их громоздкости. Вектор Н выглядит следующим образом:

rr f V2 W2 , VW f. UV W2 , UW f.

+ — Soe, — — ctg Ф Ft7

^ + J~ctg<D, o, ^- + ^ + ictg®}.

Здесь rj = r -f- !оф sin Ф. При неподвижном центре, когда £0(Ф, 0) = = const, вектор Н принимает обычный для сферической системы координат вид.

Система конечноразностных уравнений, соответствующая системе (2), имеет следующий вид:

хгп-г1 'Ут хт ftn xm _ хгп

ijk J ijk , ямтп J i+\Jk Ji-\jk I \тШ J ij+\k J ij-\k I

-----If------ + Mi*--------2lh------' + "Ч*-------2~Д<р------Г

ftn __ ftn

+ L?jkIiik+\Jijk-x -f HTjk = o , (3>

т. e. переход на следующий временной слой во внутренних точках поля осуществлялся по явной одношаговой схеме. Величина fTjk* входящая в (3),

/"* = (1 -2 & - 2 р, - 2 рз)/"* + Pi ifT-ijk +f?+w) +

+ Рз (fiJ+1 k +■ f?j—\ k) + (ffjk+i +/“*-1 )• (4)

Индекс m обозначает m-й временной слой, индексы i, j, k относятся к направлениям п, f, ft соответственно (фиг. 1, б).

На основе метода Фурье было проведено исследование устойчивости этой разностной схемы. Выделив из системы скалярное

п

уравнение для логарифма плотности можно получить, что в трехмерном случае множитель перехода |М<1, если

А2

Щ

(5)

В этом неравенстве (в случае неподвижного центра, т. е. при ?0 = const)

где индексами t, Ф, 0 обозначены соответствующие производные. Неравенство (5) дает качественную зависимость шага по времени от весовых множителей р2. Рв> входящих в схему. Эта зависимость учитывалась в расчетах. Более точно область устойчивости схемы определялась численно.

3. В нормированной системе координат тело является координатной поверхностью п = 0. Условие непротекания на поверхности тела (в случае неподвижного центра, ?0 = const) имеет вид

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

хцие производной условие совместности, выполняющееся на

волновой характеристической поверхности, и уравнение сохранения энтропии в частице газа. На фиг. 1, в показаны все узловые точки в окрестности луча (0 = сопз1, <р = соп5^, использованные в разностной схеме у поверхности тела, за исключением точек в соседних плоскостях & + Д0 на временном слое Т, участвующих в образовании производных по 0.

Опуская здесь подробное описание схемы и порядка вычислений, скажем лишь об определении значения энтропии на поверхности тела. Если уравнение сохранения энтропии записать в разностной форме для ячейки с центром Споказанной на фиг. 1, в, то, вообще говоря, можно не задавать значения энтропии на поверхности тела. При достаточно малых, но конечных размерах ячейки расчетной сетки можно надеяться на получение почти точного ее значения в процессе установления. Известно, что вопрос о точном значении энтропии на поверхности тела в трехмерном потоке газа до сих пор остается открытым. Численные эксперименты [1], проведенные при помощи метода установления, не могут дать исчерпывающего ответа на этот вопрос, однако их результаты подтверждают широко используемую рабочую гипотезу, согласно которой энтропия на поверхности тела, обтекаемого трехмерным установившимся потоком газа, по крайней мере, очень близка к энтропии в струйке тока, проходящей через прямой скачок уплотнения.

Использование разностного уравнения сохранения энтропии, записанного относительно центра в точке (3 (см. фиг. \,в), приво-

В

w, в= И7(Ф, 0, t)-G(09,t),

w

(G 4- п&) и — (Go 4" nsФ) v — (G© 4- nse)-^ф — 0.

(6)

дило на практике к очень медленному изменению плотности (и энтропии) на поверхности тела в процессе установления. Установление течения на поверхности тела сильно запаздывало по сравнению с установлением течения в других областях. Поэтому в дальнейшем ради экономии машинного времени при расчетах задавалось упомянутое выше „точное* значение энтропии на поверхности тела, а разностное уравнение сохранения энтропии было записано для точки на самой поверхности.

4. Течение на ударной волне рассчитывалось по схеме, обычно используемой в методе характеристик. Шесть неизвестных величин —v, w, Р, R и D (где D — скорость перемещения волны вдоль внешней нормали) — определяются из пяти соотношений для разрыва на волне:

рз =: Роо Ооо ;

Р “Ь Р°2 := Рса ~Г роо °оо j

, . 1 2 , 1 2 (7)

"2~а — 2~~ а°°;

и условия совместности, выполняющегося на волновой характеристической поверхности. В соотношениях (7), кроме обычных, использованы обозначения: a—V^ — D, ^ и V-— нормальная и касательная к поверхности волны составляющие скорости.

После определения скорости перемещения волны D можно найти величину связанной с ней производной ds/dT и вычислить положение волны на новом временном слое, т. е. найти величину е(Т-(-Д7') в каждой точке волны. Однако в разработанном алгоритме для определения нового положения волны применялся прием, описанный в работе [4]. Для этой цели была использована связь Зоо=1Лоо — D, которая при заданной величине а*, представляет собой нелинейное дифференциальное уравнение в частных производных относительно е. Это уравнение линеаризируется, и полученное линейное уравнение относительно приращения Де решается методом прогонки.

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

Обычно расчетная сетка в одной меридиональной плоскости состояла из сорока лучей и одиннадцати узловых точек на луче. Время расчета одного двумерного варианта зависит от того, какое решение принято за исходное, и может составлять от часа до трех часов на машине БЭСМ-6. В используемой схеме сумма весовых множителей Pi + р2 + Рз может изменяться от 0,5 (схема Лакса) до нуля. При расчетах эта сумма не превышала 0,1, в наиболее точных двумерных вариантах она была равна 0,02 (р, = р2 = 0,01, рз = 0). Иначе говоря, искусственная вязкость вводилась лишь в размерах, необходимых для обеспечения устойчивости счета.

Точность решения оценивалась по расходу газа через контрольную поверхность (в двумерном варианте — через предпоследний луч счетной сетки) и по величине интеграла Бернулли в установившемся решении. В лучших вариантах (при Р, + р2 + h — 0,02)

погрешность в расходе газа составляла 2—3%. Погрешность в интеграле Бернулли в тех же вариантах лишь в отдельных точках превышала I —2% (в точках разрыва кривизны тела при скругле-нии кромки и близких к ним), а в остальных точках не превосходила долей процента.

5. С использованием описанного алгоритма были проведены расчеты тел типа затупленных конусов с большими углами раствора, одно из которых показано на фиг. 2. Все тела рассмотренного

Фиг. 2

семейства получены из одного исходного тела — конуса с углом раствора 120°, с высотой, равной единице. Этот конус имеет переднее сферическое затупление радиусом Д! и скругление боковой кромки. Радиус скругления Д2 может быть постоянным, в этом случае скругление представляет собой часть тороидальной поверхности, а тело является осесимметричным. Если радиус Д2 зависит от меридионального угла 0, то в этом случае получаем некое несимметричное тело, аэродинамические характеристики которого могут заметно отличаться от аэродинамических характеристик осесимметричного тела.

Некоторые из полученных результатов расчета представлены на фиг. 3—7, на которых указаны значения параметров, определяющих течение. Кроме геометрических параметров Аг и Д2 к ним относятся число М набегающего потока Мао, угол атаки а, а также показатель адиабаты если газ был совершенным. В качестве реального газа, находящегося в состоянии термодинамического равновесия, был взят воздух. Словом „Воздух" на графиках обозначены результаты, относящиеся к случаям течения несовершенного газа.

і

f

На фиг. 3 показано распределение относительного давления по осесимметричному телу. В качестве абсциссы используется относительная длина дуги 5/1,45, отсчитываемая от оси системы координат вдоль контура тела. Здесь же нанесены экспериментальные данные, взятые из работы [5]. Видно, что численные результаты хорошо соответствуют экспериментальным*.

На фиг. 4 показано распределение давления по поверхности симметричных тел с различными радиусами переднего сферического затупления Aj при Моо = 4,63, а = 0, Д2 — 0,2 и 7 = 1,4 Уменьшение радиуса переднего затупления приводит к понижению давления в центральной части тела. В окрестности боковой кромки происходит обратное явление — заметное повышение давления. Можно поэтому ожидать, что аэродинамический момент, действующий на подобные тела, будет более чувствителен к форме боковой кромки у тел с меньшим радиусом переднего затупления. Тело с радиусом затупления Дь равным 2,8, уже является, по существу, сферическим сегментом с затупленной боковой кромкой. На этой же фигуре показано влияние отличия воздуха от совершенного газа на распределение давления по поверхности симметричного тела с радиусом затупления Aj = l,0. Этот вариант течения соответствует значениям Моо=10 и высоте Н — 0. Заметное уменьшение толщины ударного слоя в реальном газе по сравнению с его толщиной в совершенном газе приводит к усилению локальных свойств течения: на сферическом участке тела происходит более сильный разгон потока, после которого наблюдается очень медленное изменение давления на коническом участке.

Результаты, представленные на фиг. 5—7, бтносятся к расчету обтекания несимметричного тела с Д 4 = 1,0 при Моо — 4,63, 7=1,4. Несимметричность тела определяется зависимостью закругления кромки от меридионального угла 9. Радиус закругления был задан

формулой Д2 = 0,2 (l + sin2 . На фиг. 5 показано распределение

давления в трех меридиональных плоскостях 0—180°, 45° — 225° и 90" — 270°. Оказалось, что при нулевом угле атаки несимметрия тела, несмотря на эллиптический характер течения, проявилась в распределении давления практически лишь вне сферической носовой части тела. На этом же графике черными точками нанесено распределение давления по осесимметричному телу, контур которого совпадает с контуром несимметричного тела в плоскости 0 = 0. Различие в давлениях на осесимметричном теле и несимметричном в плоскости 0 = 0 объясняется целиком влиянием изменения формы тела. При углах атаки, отличных от нуля, распределение давления по переднему сферическому затуплению становится различным в разных меридиональных плоскостях.

На фиг. 6 показано распределение скорости по несимметричному телу в плоскости симметрии течения при различных углах атаки. Графики дают возможность проследить за перемещением критической точки и обеих звуковых точек в плоскости симметрии при изменении угла атаки. (При т = 1,4 скорость звука, отнесенная

* Это сравнение не является абсолютно строгим. В работе [5] были испытаны затупленные конусы с острой боковой кромкой. В нашем случае звуковая точка, согласно расчету, имеет ординату приблизительно 1,450 (к ней и отнесены длины дуг; ордината точки сопряжения скругления с конической поверхностью равна 1,432), а зона влияния скругления на течение в дозвуковой области отсутствует.

к максимальной скорости, равна 0,408.) Четырьмя тонкими вертикальными линиями здесь отмечено положение четырех точек сопряжения переднего сферического затупления и поверхности скруг-ления кромки с коническим участком поверхности тела (5сопр = 0,525 и 1,600 при 0 = 0; 5сопр = 0,525 и 1,225 при 0=180°). В рассмотренном диапазоне углов атаки градиент скорости в критической точке изменяется очень мало.

Наконец, на фиг. 7 представлена зависимость коэффициента продольного момента тг и аэродинамического качества К несимметричного тела от угла атаки. Коэффициент продольного момента вычислен относительно точки, расположенной на оси на расстоянии 0,7 от вершины исходного конуса (см. фиг. 2), и отнесен к характерной длине, равной единице, и характерной площади, равной я. При вычислении момента и сил, действующих на это тело, учитывалось только давление в рассчитанной области течения. Расчет был проведен до угла атаки а = —15°. При этом угле атаки АГ= — 0,296. Из зависимости т2(а) следует, что при а ^ — 20° это несимметричное тело будет самосбалансировано с аэродинамическим качеством, близким к —0,4. Следовательно, введение небольшой несимметрии в форму тела можно эффективно использовать при проектировании летательных аппаратов с заданными аэродинамическими характеристиками на режиме самобалансировки. Летательные аппараты подобного рода были недавно предложены и экспериментально исследованы в работе [6].

Авторы выражают благодарность Ю. Я. Михайлову за полезные советы, высказанные им при подготовке статьи к печати.

ЛИТЕРАТУРА

1. Л ю б и м о в А. Н., Русанов В. В. Течения газа около тупых тел. М., .Наука”, 1970.

2. Moretti G. and Bleich Q. Three-dimensional flow around blunt bodies. A1AA J., vol. 5, No 9, 1967.

3. Липницкий Ю. М., Михайлов Ю. Я., Савинков К. Г. Расчет пространственных течений идеального газа без плоскости симметрии. „Изв. АН СССР —МЖГ“, 1972, № 3.

4. Б а б е н к о К. И. и др. Нестационарное обтекание головной

части затупленного тела идеальным газом. Препринт ИПМ АН СССР, 1969. ,

5. Stallings R. L. and Tudor D. H. Experimental pressure distributiops on a 120° cone at Mach numbers from 2.96 to 4.69 and angles of attack from 0 to 20°. NASA TN D—5054, 1969.

6. A p т о н к и н В. Г., Петров К. П. Влияние несимметричного уменьшения площади миделя тела вращения на его аэродинамические характеристики. „Ученые записки ЦАГИ, т. 4, № 4, 1973.

Рукопись поступила 16jXlI 1972 г.

2—Ученые записки ЦАГИ № 5