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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XIII 19 82

№ 5

УДК 629.735.33.015.3.024

РАСЧЕТ ОКОЛОЗВУКОВОГО ОСЕСИММЕТРИЧНОГО ОБТЕКАНИЯ ТЕЛ ВРАЩЕНИЯ

В. В. Вышинский

Предлагается основанный на применении алгоритма релаксации конечно-разностный метод расчета околозвукового осесимметричного обтекання тел вращения. Отличительными особенностями метода являются высокие скорость сходимости и устойчивость, что позволяет с небольшими затратами машинного времени вести расчет обтекаиия тел вращения в диапазоне от нулевых до звуковых скоростей набегающего потока. Используемое преобразование переменных дает возможность рассчитывать обтекание как конечных, так и незамкнутых полубесконечных тел.

В работе [1| с помощью метода релаксации рассчитывается бесциркуляционное обтекание пространственных тел околозвуковым потоком невязкого газа. В качестве исходного уравнения движения используется полное уравнение относительно потенциала, применяется преобразование координат, позволяющее вести расчет обтекания как замкнутых, так и полубесконечных тел. Повышению эффективности метода способствуют применение алгоритмов удвоения сеток [2], экстраполяции к пределу по Ричардсону [3], выбор наилучшего исходного приближения, подбор параметров сходимости и оптимальных сеток.

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

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

1. Рассматривается осесимметричное стационарное изоэнтро-пийное течение идеального газа. Решение в продольной меридиональной плоскости ищется в декартовых прямоугольных координатах х, у, где ось х связана с осью симметрии тела вращения. Уравнение относительно потенциала течения Ф в безразмерной форме имеет вид:

где а—скорость звука, и = Фх, V = Ф^,— проекции вектора скорости на оси х и у. Уравнение Бернулли

замыкает систему уравнений.

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

Уравнение (1) переписывается в полярных координатах (г, 0) внутри круга единичного радиуса. Для отображения внешности продольного контура тела на внутренность круга используется преобразование работы [6]

где а-=г-ел, г — х + — угол задней кромки продольного

контура тела вращения. Коэффициенты ряда сп = ая + 1&л находятся в результате итерационной процедуры.

Применение формул тригонометрической интерполяции совместно с быстрым преобразованием Фурье позволяет существенно сократить число элементарных операций на каждой итерации и делает процедуру отображения быстрой и эффективной.

В том случае, если

за телом моделируется след постоянной толщины.

В плоскости круга уравнение (1) имеет вид:

(а2 — и2) ? 88 — 2 г и + г2 (а3 — V2) ®,г — 2 + г (аг иг — 2гТг) ?г +

2 + х — 1 ~ 2

ііг Аз

— - а)1-*''* ехр

геэ = (1 — Е/к — с,) ехр С0 ф О,

+ 4- ("2 +^2) (иГСЪ + ^ Фг) + а2 Ц7 (3_ и + = 0, (2)

где ® — регулярная при г -+ 0 часть полного потенциала, Ф = <р4-+ cos 6/г, W = r2\dzjda\ — модуль якобиана отображающей функции, /г = 1 [/"срв — sin 0J, х? = W''—1 [/^2 — c°s 6J — компоненты ско-

рости в плоскости круга.

Контур тела является поверхностью уровня r=const=l, условие непротекания записывается как

?,|г-1 = cos 6.

Условие на бесконечности имеет вид

?|г=0 = О.

Решение задачи ведется конечно-разностным методом в переменных г, 0, выбор которых обеспечивает естественное сгущение расчетной сетки в области больших градиентов параметров.

Уравнение (1) на оси у — 0 имеет особенность, которая раскрывается посредством соотношения

lira Ф /у = Ф .

у^О

и последний член в уравнении (2) при у-»-О заменяется на

г3 . о /cos 6 \ , г!

у- РИ-------_j + _ye>y r6V

Уравнение (2) аппроксимируется конечно-разностным оператором смешанного эллиптико-гиперболического типа, обеспечивающим второй порядок аппроксимации в дозвуковой области поля течения и первый порядок в сверхзвуковой зоне [6]. С целью повышения устойчивости метода релаксации используется предложенная в работе [5] следящая схема.

В уравнении выделяются группы членов, соответствующих вторым производным в направлении s вектора скорости К=(и, v) и перпендикулярном к нему направлении п

(a*-V'*)<Pw + a«v+... = 01 (3)

где .

= («2 Yee -г 2гищг е + г2 v'1 <¥n)l V2, (4)

= (®2 - Zruvfrtt + rs и2 <prr)f V2. (5)

В узлах разностной сетки, где уравнение (2) является гиперболическим, для аппроксимации всех вторых производных в (4) используются односторонние разности, а для тех же производных в (5) — центральные разности. Центрально-разностные операторы используются также для аппроксимации всех первых производных [в (3) члены с первыми производными опущены], независимо от типа уравнения в расчетном узле, и для аппроксимации вторых производных в узлах разностной сетки, где уравнение (2) является эллиптическим. Односторонние разности берутся в направлении, противоположном вектору V, и в зависимости от знака v при этом используются различные узлы расчетной сетки.

Подстановка разностного оператора в уравнение (2) приводит к нелинейной системе алгебраических уравнений, которая решается с помощью алгоритма релаксации по лучам 9( = const. Система уравнений записывается на луче i, а ее решение находится в итерационной процедуре, где на каждом шаге обращается матрица

соответствующей линейной системы, и производится последующий пересчет коэффициентов. Релаксационный процесс идет шаг за шагом для последовательных лучей I = const в направлении от головной части тела к корме. Параметр релаксации со из соображений устойчивости выбирается в диапазоне 1<;ш<2.

Система уравнений на луче имеет четырехдиагональную матрицу (как следствие применения следящей схемы). Для сохранения трехдиагональности матрицы и диагонального преобладания вводятся значения ? на луче i на предыдущей и текущей итерациях [5], что позволяет рассматривать итерационный процесс как процесс установления по некоторой искусственной временной переменной t, а разностные формулы — как аппроксимацию членов с производными по времени. Обращение трехдиагональной матрицы легко осуществляется с помощью устойчивого алгоритма прогоики.

В некоторых случаях с целью ускорения сходимости и повышения устойчивости при больших местных числах М в разностное уравнение вводится дополнительный демпфирующий член:

= p-55-(«?#( +гг?(6)

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

Координаты контура исходного тела задаются в виде таблицы и интерполируются с помощью кубического сплайна. Используется раздельное параметрическое представление переменных л (и.) и 3^ (н-)-Параметр [0, и определяется как

(1= arccos ^1 —

где s—текущая, a s0—полная длины дуги. Для устранения ошибки аппроксимации и, в частности, неточности исходных данных используется алгоритм сглаживания [7].

После выполнения отображения внешности контура на внутренность круга для всех узлов расчетной сетки находятся якобиан преобразования W и производные уг ув или уг, угъ (для узлов сетки на оси у = 0). При реализации метода на ЭВМ эти величины наряду со значениями ф во всей расчетной области хранятся в оперативной памяти.

В качестве исходного приближения для итерационного процесса метода релаксации могут использоваться невозмущенный набегающий поток <р = 0, предыдущий (по числу М«>) расчет для данного тела, расчет при том же Моо для другого тела или точное решение для обтекания потоком несжимаемого газа эллипсоида вращения [8].

Применение алгоритма удвоения, как и в методе [1], позволяет существенно сократить время вычисления. Максимальное число итераций на сетке после удвоения задается равным 200, до удвоения — 100. Используются одинаковые значения параметра релаксации ш для расчетов на сетках до и после удвоения. Исходя из численных экспериментов, параметр <и выбран для точек разностной сеткн, соответствующих эллиптичности уравнения, равным 1,6, для точек разностной сетки, соответствующих гиперболичности уравнения, равным 1. Значение параметра J3 в формуле (6) обычно задается равным 0,2 и несколько большим для расчетов при Моо ->■ 1.

При решении систем алгебраических уравнений на лучах * = const можно ограничиться однократным обращением матриц линейных систем.

В качестве стандартной сетки для массовых расчетов используется сетка с числом узлов 41X31. В тех случаях, когда не требуется высокой точности, можно ограничиться применением более редкой сетки с числом узлов 41 X 17* Применение последней сетки оправдано также в случае расчета обтекания тел с затупленными иосовой и кормовой частями. В частности, она используется в последнем из приводимых в п. 2 примеров.

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

2. С целью демонстрации особенностей вычислительного метода и возможностей расчетной программы приводится ряд примеров построения поля течения при осесимметричном обтекании тел вращения.

В качестве первого примера рассматривается околозвуковое обтекание тела вращения с относительной толщиной 14,4%, контур которого в хвостовой части оканчивается точкой возврата. На рис. 1 представлены контур продольного сечения тела и распределение чисел М на теле для последовательных режимов от Мес = 0,3 до Мж =1,01. Высокая устойчивость метода позволила провести расчет обтекания вплоть до Мм=1,01, однако используемое преобразование координат, по-видимому, не позволит (в отличие от преобразования работы [5]) использовать метод для расчета обтекания при числах Мж, существенно больших единицы. На этом же рисунке приведено поле линий М = const для решения с Моо = 0,975.

Интегрирование распределения давления на теле позволяет определить величину волнового сопротивления Сх. Зависимость С*(Мос) представлена на рнс. 2. Как видно, волновое сопротивление становится существенным лишь при -+■ I, когда на эпюрах М(х) появляются заметные скачки уплотнения.

При увеличении числа М дозвукового набегающего потока Мсо 1 распределение давления вдоль обтекаемой поверхности и значения других распределенных параметров до скачка уплотнения меняются слабо с изменением скорости набегающего потока, хотя сам скачок довольно быстро перемещается в кормовую область. Асимптотическая теория предсказывает, что при вариации Моо вблизи Моо=1 параметры газа как в поле течения, так и на поверхности тела будут отклоняться от своих предельных значений для звукового потока пропорционально |Моо—1|5/3.

Пте М-сбт1

Рис. 2

В частности, для волнового сопротивления Сх в работе [9] получено разложение

ДС, = С, (1)- Сл (Ми) = С, (1 - Моо)1'3 + С, (1 - Мо;,)2'3 + . . .,

справедливое для тел малой и умеренной относительной толщины. На рис. 2 с целью его проверки приведена зависимость Сх от (1--МС0)2а. Линейный характер зависимости говорит о том, что

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

в данном случае собственное значение ^ = 0 и разложение начинается со второго члена.

Во втором примере рассматривается обтекание полубесконеч-ного тела, контур которого получен посредством наращивания толщины вытеснения по нормали на контур первого тела и осевую линию следа при числе Ие = 22,5* 10е. Обтекание полученного полутела рассчитывается при тех же значениях Мс*,. На рис. 2 приведены зависимости Сл(Мю) и Сх [(1 — Мго)2/3[, Для полутела получаются меньшие значения волнового сопротивления при тех же числах М набегающего потока. Скачки уплотнения одновременно с уменьшением интенсивности смещаются вперед. В этом случае при Моо>0,95 также наблюдается линейный характер зависимости Сх от (1 — Мот)ад. Коэффициент пропорциональности С2 зависит от формы тела и принимает в последнем случае меиьшее значение. Аналогичные выводы получены в работе [10] для профиля

в околозвуковом потоке. Контур полутела представлен на рис. 3, где приведены линии М = const и линии тока'для режима с Мм = = 0,975,

Возможность расчета обтекания „тупохвостых“ тел с большой относительной толщиной демонстрируются на примере расчета безотрывного обтекания сферы. На рис. 4 представлено распределение чисел М на сфере при числах М набегающего потока Mo, = 0,3-s-0,9. На рис. 5 результаты расчета при Мж = 0,8 [распределение давления Ср(х) и поле линий М = const] сравниваются с результатами, представленными в работах [5] и [11]. В методе [И] применяются уравнения Эйлера, расчет ведется по методу установления и выполнен на достаточно редкой сетке, чем, по-видимому, объясняется расхождение последних результатов с результатами работы [5] и полученными предлагаемым методом.

Описанная вычислительная процедура реализована на ЭВМ БЭСМ-6 в виде программы расчета, написанной на языке

tp

-2

Пане M^const

схч,т

M=8,B

Рис, 4

Рис, 5

ФОРТРАН. Потребное время решения на стандартной сетке с числом узлов 41X31 занимает около трех минут. Использование графопостроителя позволяет оперативно и с высокой точностью получать результаты в графическом виде. Все расчеты могут сопровождаться построением поля линий тока и линий уровня М = const, Ср = const И Р = const.

Автор приносит благодарность С. В. Ляпунову за полезное обсуждение результатов работы.

ЛИТЕРАТУРА

1. Вышинский В. В. Расчет околозвукового обтекания трехосного эллипсоида. .Ученые записки ЦАГИ“, т. XI. № 6, 1980.

2. Федоренко Р. П. Релаксационный метод решения разностных эллиптических уравнений. „Ж. вычисл. матем. и матем. физ,“, т. 1, № 5, 1961.

3. М а р ч у к Г. И., Шайдуров В. В. Повышение точности решений разностных схем. М., „Наука", 1979.

4. W h 11 с о ш b R. Т. A study of the zero-lift drag rise characteristics of wing body combination near the speed of sound. NACA TR 1273, 1956.

5. Jameson A., South J. C. Relaxation solutions for inviscid axisymmetric transonic flow over blunt or pointed bodies. A1AA Computational Fluid Dynamics Conference, Palm Springs, California, 1973.

6. Bauer E., Garabedian P., Korn D„ Jameson A. Supercritical wing sections II. Springer-Verlag, Berlin—Heidelberg—New York, 1975. (Lecture Notes in Economics and Math. Syst. 108).

7. Вышинский В. В, Применение метода градиента к задаче сглаживания поверхностей. „Ученые записки ЦАГИ‘, т. IX, №6, 1978.

8. К о ч и н Н. Е., К и 6 е л ь И. А., Розе Н. В. Теоретическая гидродинамика. Ч. I. М., Физматгиз, 1963.

9. Днесперов В. Н„ Л и ф ш и ц Ю. Б. О сопротивлении тел вращения при трансзвуковых скоростях потока. ПММ, т. 39, вып. 2, 1975.

10. Л и ф ш и ц Ю. Б., Р ы ж о в О. С. О трансзвуковых течениях около несущего крылового профиля. „Изв. АН СССР, МЖГ‘,

1978, № 1.

11. Липницкий Ю. М., Лифшиц Ю. Б. О расчете обтекания тел вращения трансзвуковым потоком. ПММ, т. 34, вып. 3, 1970.

Рукопись поступила 25ЦП 1981 г.