Исследование методов сквозного счета для задач сверхзвуковой аэродинамики Текст научной статьи по специальности «Физика»

________УЧЕНЫЕ ЗАПИСКИ Ц А Г И

Т о м VII 197 6

№ 1

УДК 518:517.9:533.6.011

ИССЛЕДОВАНИЕ МЕТОДОВ СКВОЗНОГО СЧЕТА ДЛЯ ЗАДАЧ СВЕРХЗВУКОВОЙ АЭРОДИНАМИКИ

А. П. Косих, А. Н. Минайлос

На основе работ [1—2] разработан стационарный аналог метода [1] для расчета невязких сверхзвуковых течений. Его точность сопоставляется с точностью метода С. К. Годунова [2—3] на примерах простых точных решений для сверхзвуковых течений (контактного разрыва, течения Прандтля — Майера и ударной волны). Обсуждается потеря точности схем при расчетах течений разрежения.

1. Первая часть исследования во многом аналогична работе [4] и поэтому описана кратко. С целью выбора эффективной конечноразностной схемы на модельной задаче формирования и движения разрыва исследуется ряд схем [1, 2, 5—9], имеющих порядок аппроксимации от 1 до 4. Рассматривается уравнение

м^ + иих = 0, 0<£<1, 0<*< 1,

(1, О<ЛГ<0;

«(0, X) — ,2-, о<*<20;

.0, 20<*<1.

В точном решении при£ = 0 формируется разрыв. Решение имеет вид: при £ < 0

1, О<*<г+0;

26 — X

при і 0

и {і, Х) =

и(г, х) = .

в — ( 0.

¿ + 0<*<20; 20<ЛГ< 1;

1, 0 <*<-¿-(¿ + 36);

10. 1

2-(г + 30)<*<1.

Результаты расчетов, выполненных по различным схемам, рассматриваются в моменты времени ¿ = 0 и £>0 для исследования формирования и движения разрыва. На фиг. 1 для двух моментов времени ¿ = 0=0,2 и £ = 0,4 приведены результаты, полученные различными методами. Сплошной линией изображено

точное решение; результаты, полученные по схеме [1] — черными кружками, по методу [3] — белыми кружками, по схеме Мак-Кормака, описанной в [6] и по методу [5]—штриховой линией, для схем третьего порядка аппроксимации [6—7] — пунктирной линией, по схеме четвертого порядка [9] — крестиками.

Основные выводы состоят в следующем. Определяющим для схем является порядок аппроксимации. Схемы, имеющие один и тот же порядок аппроксимации, дают количественно и качественно почти одинаковые результаты. Таковы, например, схема Мак-Кормака [6] и схема [5], имеющие второй порядок аппроксимации, или схема Кутлера — Ломакса — Вармина [6] и схема [7], имеющие третий порядок (см. фиг. 1).

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

Скорость движения фронта разрыва при использовании схемы третьего порядка [8] меньше точного значения. Это единственная из рассмотренных схем, в которой фронт отстает от точного решения.

Свойством монотонности обладают из рассмотренных только схемы [1,3]. Это свойство (чрезвычайно важное при анализе результатов расчетов сложных задач), а также хорошие диссипативные и дисперсионные характеристики схемы [1] определили выбор этой схемы для дальнейших исследований.

2. Аналог схемы [3] для расчета стационарных сверхзвуковых течений предложен в работе [2]. Вместо временных распадов разрыва рассмотрена автомодельная задача о столкновении двух

Ю

сверхзвуковых потоков. Решение получено в линейном приближении относительно угла поворота потока, т. е. в приближении Аккерета. В таком виде метод применяется авторами [2] к решению широкого класса задач, в которых углы поворота скорости течения невелики. При разработке стационарного аналога схемы [1] за основу была взята работа [2]. Однако для течений с достаточно большими углами поворота потока (—20°—30°) точность расчета распада разрыва по формуле Аккерета оказалась недостаточной, несмотря на то, что в схеме [1] на монотонном решении распады слабее, а следовательно, они считаются „по Аккерету“ точнее, чем в схеме [3].

3. Расчет столкновения двух плоских сверхзвуковых потоков осуществляется по итерационной схеме. В схеме используются точные соотношения для расчета ударных волн, как и в [2], и точные соотношения для расчета течения Прандтля — Майера. Начальные значения для итераций получаются „в приближении Буземана“, т. е. с точностью второго порядка относительно угла поворота потока. Расчеты показали, что необходима высокая точность сходимости итераций. Так, относительное отличие давления в двух последующих приближениях

П Рп— 1

д =

Рп-1

(При Д =10“3 в области начальной характеристики течения разрежения ошибки в плотности достигали 20%).

4. Как известно, в методе Годунова [3] используется кусочнопостоянная аппроксимация искомых функций внутри элементарной ячейки. Полученную однородную монотонную схему первого порядка точности условно обозначим Г. Позже в работе [1] была предложена схема повышенной точности (назовем ее К1) для решения нестационарных задач газодинамики. В основу этой схемы положена кусочно-линейная аппроксимация функций с использованием специально выбираемых минимальных градиентов по пространственной переменной. В [1] показано, что таким образом построенная схема монотонна, так же как и схема Г, но неоднородна и по пространственной переменной порядка выше первого. Схема К1 в отличие от многих схем высокого порядка, применяемых для сквозного счета, не дает осцилляций и тем самым существенно уменьшает искажение истинной картины течения. Во многих случаях с помощью осциллирующих схем вообще невозможно получить приемлемое численное решение задачи. В. П. Колган рассмотрел также модификацию схемы К1 с более простым способом построения минимальных производных в ячейке (обозначим эту схему К2). Мы построили стационарный аналог этой схемы для двумерных и пространственных течений. Вариант для двумерных течений рассмотрен ниже. Запишем систему нелинейных дифференциальных уравнений плоского течения идеального газа (* = 1,4) в декартовых координатах и дивергентной форме:

~5х^Ц)+ -§у№У)--=0\ -¿(Р+^)+^(/?У1/) = 0; -¿(.^^)+Д(Р + ^^)=0.

(1)

Систему замыкаем условием постоянства полной энтальпии

U2+V2 + ^lYR = const’ (2)

здесь Р, R, U, V — давление, плотность и компоненты скорости газа соответственно.

Мы будем рассматривать задачи, для которых система уравнений (1) — (2) X — гиперболична. Введем на плоскости XOY прямоугольную сетку X{_j_= (i—Y._± = (у — 4") Л Y’ i==l>

2, . . . ; j = 1, 2, . . . .Центры прямоугольных элементарных ячеек (i, у) имеют соответственно координаты Xt = iAX, Yj = уДУ. В ячейке (/, У) параметры течения представлены линейной зависимостью как функции от У. Например, для плотности R внутри ячейки вводим аппроксимацию

R(X, Г) = £,у+Щ \Y-Yj).

Если

dR

. .= ¿y min {| Rt, J - Ги, у_, I, | I, I /?Л/+1 - Rit ,■ |} ,

2

dY

то получается схема K2. Функции P(X,Y), U(X,Y), V(X, У) аппроксимируются аналогично. Очевидно, что в результате такого кусочно-линейного представления функций внутри ячеек (i, j) на границах ячеек образуются разрывы. Как уже отмечалось, эти

разрывы слабее, чем в методе [3]. Параметры на нижней^/, у---

и верхней ^г, у + границах ячейки определяем из автомодельного решения задачи о взаимодействии двух равномерных сверхзвуковых потоков.

Применяя интегральные законы сохранения к элементарной ячейке и схему К2, получим систему разностных уравнений

(RUY-' - (ВД. I/). /+^ - (R V).

(Р + RU*y>J = (Р + RlP)i,, - *¡*^RUV)U+1_- ;

(RUV)t'J-{RUV)l,J-^P + RV')i¡¡+j-(P + RV2)i .;

(£/i. /). + (vuy + = const.

Верхние индексы соответствуют слою X + АХ, а нижние — слою*.

Схема К2 на гладких функциях „почти“ второго порядка точности по координате У и первого порядка по X. Схема устойчива при условии, что возмущения с границ ячеек при переходе со слоя X на слой Х-{-АХ не достигают центра ячеек.

5. Экспериментальное исследование описанной выше схемы К2 и сопоставление ее со схемой Г [2] проведено на примерах расчета элементов плоского сверхзвукового течения: контактного разрыва, ударной волны и течения Прандтля — Майера. Все эти течения создавались при столкновении двух плоских равномерных потоков. В направлении Y в поле течения равномерно располагалось 200 (или 400) узлов расчетной сетки. Граница двух

потоков в начальном сечении *==0 располагалась между 100(200) и 101(201) узлами сетки.

Приведенные на фигурах схемы течения имеют следующие обозначения: пунктирные линии — границы течения разрежения, штриховые — контактные разрывы, двойные линии — ударные волны. На фиг. 2—4 сплошными линиями показано точное решение, штриховыми с черными кружками — решение, полученное по схеме Г, штрих-пунктирными с белыми кружками — по схеме К2. Каждый кружок соответствует узлу расчетной сетки. Результаты расчета столкновения двух потоков (кроме контактного разрыва) строились по автомодельной для точного решения переменной Y/Х. Это позволило, имея решения в различных сечениях по X, оценить точность результатов при увеличении числа узлов в 2, 4, 8, 16 раз.

Если контактный разрыв не направлен вдоль оси X, он размывается при расчетах. С увеличением значения X разрыв размывается все сильнее и охватывает все большее число узлов сетки. Схема К2 меньше размывает разрыв, чем схема Г; это отличие возрастает с ростом X. С увеличением угла между осью X и направлением скорости (т. е. разрыва) степень „размазывания“ увеличивается. На фиг. 2 показано распределение плотности R (отнесенной к Roo) и энтропийной функции S = P/R* (Р отнесено К Roo í^max) в части сложного течения состоящего из течения Прандтля — Майера, контактного разрыва и ударной волны (см. схему). Показана часть течения. Она охватывает равномерный поток после течения разрежения, контактный разрыв, равномерный поток за волной и ударную волну. Видны преимущества схемы К2.

Наилучшим образом монотонные схемы описывают ударные волны. При повороте потоков с Моо = 3 в ударной волне на 30° ошибка в энтропийной функции для схемы С. К. Годунова не превосходит 6%, а для схемы В. П. Колгана — 3%. Эта ошибка не устраняется при увеличении числа узлов, описывающих течение, а локализуется в окрестности контактного разрыва, разделяющего два

т.,

%

100

50

О

-50

Фиг. 4

потока на одном-двух узлах сетки. Вызвана она потерей точности аппроксимации схемы при расчете начальных шагов. ^

Аналогичным образом возникают погрешности значений энтропии е = 5^точн 100% в области постоянных значений параметров

¿ТОЧН

за течением Прандтля — Майера. Однако величина погрешности с увеличением угла поворота потока ¡3 начинает быстро расти и достигает при повороте потоков на 30° (Моо = 3) 66% для схемы Г и 71% для схемы К2 (фиг. 3 и 4). Такая величина ошибки

Фиг. 3

*9 /

/

/

Щ / /

У /

| У /

■20 **

<

0- 1

70° 20

V

•-1 ;г

V •и ч

к N

\\ \\ 0,5 X

V —• N

1

в функции 5 вызвана разными знаками относительных ошибок ер в давлении и sR в плотности в окрестности точки поворота потока.

Приобретая погрешность в области, где схема теряет точность, искаженная энтропийная функция в дальнейшем, как ей и положено в изэнтропическом течении, сохраняет свое значение, давление— выравнивается (ошибка в давлении уменьшается с 20—30 до

3—4%). Сохраняется погрешность в значениях плотности. Таким образом, возникает „схемный“ энтропийный слой (фиг. 4) с ошибкой в плотности еЛ, достигающей 30%.

В многочисленных работах, использующих в расчетах метод С. К. Годунова, авторы обычно не исследуют ошибки энтропийной функции или, во всяком случае, в публикациях опускают вопрос о возникновении фиктивного энтропийного слоя. При достаточно больших углах поворота течения расширения этот слой определяет точность расчета всего поля (за исключением величины давления). К сожалению, . схема К2 дает за течением расширения примерно такую же ошибку в энтропийной функции, как и схема Г. Само же течение расширения описывается при помощи схемы К2 лучше (фиг. 3). В расчетах, результаты которых представлены на фиг. 3, давление отнесено к R«з а1(а* — критическая скорость звука).

На фиг. 4 приведены кривые относительных ошибок как функций угла поворота потока р и переменной *(Р = 30’). До значения р—10° ошибок в давлении практически нет, а плотность и энтропия определены достаточно точно. С увеличением угла Р^10° кривые ошибок и г5 резко возрастают, хотя значение вр не превышает 4% (штриховые линии соответствуют схеме Г, а штрих-пунктирные — К2).

Рассмотрим, наконец, еще один, более сложный пример, выявляющий преимущество разработанного варианта схемы [1]. Исследуем взаимодействие контактного разрыва с автомодельным течением, представленным на фиг. 2. При этом в слое газа II (фиг. 5) возникает система постепенно ослабевающих отраженных волн. Первая из них обозначена цифрой 6, вторая—-цифрой 7. На фигуре приведено распределение плотности в различных сечениях X = const в области течения, выделенной на фигуре окружностью.

Цифрами в кружках на верхней части фигуры обозначены отдельные разрывы в поле течения, указанные на схеме внизу фиг. 5. Линиями без точек обозначены результаты, полученные по схеме К2. Кривая с черными точками (каждая точка соответствует расчетной ячейке) получена по схеме Г и соответствует сечению * = 2,4. Слабый скачок 6 размыт схемой Г настолько, что положение его трудно определить. Второй отраженный скачок 7, выявленный схемой К2, в результатах схемы Г уже совсем отсутствует.

Выявленные закономерности для схем К2 и Г проявляются и в результатах расчетов пространственных течений. На фиг. 6 показано распределение параметров в плоскости симметрии течения около плоской бесконечно тонкой треугольной пластины (угол стреловидности х — 44,7°, число Моо = 2,94, угол атаки а — 12°). Пластина расположена в плоскости Y ==0,4 (см. рисунок в верхней части фиг. 6). Как и на фиг. 4, штриховые линии соответствуют методу Г, а штрих-пунктирные — К2. В рассматриваемом случае ударная волна присоединена к передним кромкам.

OJ

OJ

0Л

0,2

О 0,05 OJO 0,15 i

о. . І* to is a

0,052 OflSï 0,056 S

Фиг. 6

Прежде всего отметим, что метод К2 точнее описывает ударную волну на наветренной стороне пластины (F = 0,14-f-0,18), размывая ее на меньшее количество узлов.

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

Схема К2 лучше описывает и течение расширения над пластиной. Точное значение границы возмущенной области течения над крылом (У = 0,9) показано на фиг. 6 черным треугольником. По схеме Г возмущение распространяется до 7=1,0, по схеме К2—до Y = 0,9.

Ошибки в энтропийной функции S = Р/R* в методе Г помимо причин, создающих фиктивные энтропийные слои, вызываются ростом энтропии в распадах разрывов на границах ячеек. Поскольку в методе К2 на монотонном решении распады значительно слабее по интенсивности, то и энтропийные ошибки в поле меньше.

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

ЛИТЕРАТУРА

1. Колган В. П. Применение принципа минимальных значений производной к построению конечноразностных схем для расчета разрывных решений газовой динамики. „Ученые записки ЦАГИ-, т. III, № 6, 1972.

2. Иванов М. Я., Край ко A. H., Михайлов Н. В. Метод сквозного счета для двумерных и пространственных сверхзвуковых течений. Ж. вычислит, матем. и матем. физ., т. 12, № 2, 1972.

3. Годунов С. К. Разностный метод численного расчета разрывных решений уравнений гидродинамики. Матем. сб., т. 3,

№ 47, 1959.

4. Taylor Т. D., Masson В. S., Ndefo E. A study of difference methods for solving viscous and inviscid flow problems. J. Comp.

Phys., т. 9, 99 (1972).

5. Lax P. D., W e n d г о f f В. Systems of conservation laws. Com.

Rure. Appl. Math., vol. 13, 1960.

6. Kulter P., Lomax H., Warming R. Computation of space shuttle flow fields using noncenterid finite-dlfference schemes.

AIAA Paper N 72-193, 1972.

7. Русанов В. В. Расчет и исследование многомерных течений газа методом конечных разностей. Диссертация на соискание ученой степени д. ф.-м. н. М., ИПМ, АН СССР, 1968.

8. AbarbanelS., ZwasG. Third and fourth order accurate schemes for hyperbolic equations of conservation law form. Journ. Mathem. of Comput., vol. 25, N 114, 1971.

9. Б а л а к и н В. Б. О методах типа Рунге-Кутта для газовой динамики. Ж. вычисл. матем. и матем. физ., т. 10, № 6, 1970.

Рукопись поступила 26/VII 1974 г.

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