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

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

Том XIX 1988

№ 3

УДК 519.6

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

В. М. Нейланд

Предложен новый численный метод расчета трансзвукового обтекания плоских тел идеальным газом. Использование некоторых фундаментальных свойств течений при -*• 1 и точных решений позволило повысить быстродействие метода в среднем на два порядка по сравнению с существующими численными методами интегрирования системы уравнений Эйлера. При этом удается избежать «размазывания» скачка уплотнения, присущего всем конечно-разностным схемам. Приводятся примеры расчета обтекания профиля, дано их сравнение с экспериментом и расчетами других авторов.

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

Интегрирование системы уравнений Эйлера дает более точное описание картины течения. При этом чаще всего используется метод установления [3], [4], когда нестационарные уравнения Эйлера интегрируются при стационарных граничных условиях до получения установившегося решения. Реже применяются другие подходы: метод релаксации [5] и метод интегральных соотношений [6]. Все эти методы требуют

весьма значительных объемов оперативной памяти (приблизительно в четыре раза больше, чем методы, основанные «а использовании потенциалов) и быстродействие их в среднем на порядок меньше. Подробное изложение этих вопросов дано в [7], [8].

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

2. Рассмотрим симметричное обтекание профиля ую(х) = хЦх) относительной толщины т<1 с хордой Ь потоком идеального газа. Такое течение описывается системой уравнений Эйлера (см. например, £9]):

здесь и, V, р, р — составляющие скорости, давление и плотность соответственно, х, у — декартовы координаты.

Дополним систему (1)—(3) уравнением сохранения энтропии вдоль линии тока

где С=const, индексы ± соответствуют параметрам до и после скачка, у — показатель адиабаты, и граничным условием

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

в ряд по малому параметру є (который будет определен ниже): и = и00(1+ гщ + е2м2 + .р*~р№{1 + е/», + 22/?2 + ...),

Р = Роо(1+£Рі + Є2Р2 + •••), « = 'И*00(«1 +^2 +-)•

Индекс «оо» соответствует невозмущенному потоку. Подставляя эти разложения в систему (1) — (5) и приравнивая члены при одинаковых степенях е, получим для первого приближения однородную линейную систему уравнений:

д(рц) д(рр) = о

дх '■ ду і ’

(1)

да , да др

?ute+Pvlh=-l£r>

(2)

(3)

(5)

и y — Представим решение системы (1) — (5) в виде разложения

/>1 = ГР 1. Ul н---------------772 Р\ = 0>

7Мо°

00

dv{

1 с ь дрг

ТМ^ Н ду

"С

дх

Здесь Me» — число Маха набегающего потока.

Для невырожденности системы надо потребовать

<7>

Во втором приближении из уравнения неразрывности следует

м~~ 1 дР1 | с2 д дх дх

£ -Ь— + е2 Tr (Pi “1 + Р2 «.) +

тЬ Г dt/j

+ e-^L-(®2 4- Pl^i)

= 0. (8)

Чтобы сохранить в уравнении (8) при Моо-*-1 все три слагаемых и тем самым рассмотреть общий случай, надо потребовать

е — (М1с — 1)/М», (9)

(Ю)

На основании (7) и (10) получим систему оценок для определения е и Ь/Н по заданному т:

е — т2/3, (11)

ЫН — т» з, (12)

а из соотношения (9) — параметр трансзвукового подобия Кармана

М1— 1

-~1. (13)

Mf

Соотношение (12) представляет принципиальный интерес. Из него следует, что при Моо—->-1 поперечные размеры Я возмущенной области течения значительно превосходят соответствующие продольные размеры Ь. Это означает, что в пределах поперечной полосы шириной Ау~Ь изменение параметров течения в поперечном направлении будет много меньше соответствующих изменений в продольном направлении

-J^r) • Поэтому в пределах такой полосы течение можно рассматривать как близкое к одномерному, для которого точное решение, соответствующее одномерной газовой динамике в канале с заданным распределением площадей сечений и расходов через стенки, дает хорошее первое приближение, легко уточняющееся в дальнейшем.

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

Метод будет тем точнее, чем ближе к скорости звука будет скорость набегающего потока.

Из этой схемы выпадают области течения, где местная скорость далека от скорости звука, и все выводы, основанные на; соотношении (12), становятся неверными. При обтекании профиля — это области передней и задней критических точек. Течение в окрестности этих точек нужно рассматривать отдельно и «склеивать» с остальным решением, как это сделано, например, у Ван-Дайка [10] для тонкого затупленного профиля в несжимаемой жидкости.

3. Основными элементами предлагаемого метода являются расчет одномерного течения в канале переменного сечения с заданным расходом через границу и расчет течения в полубесконечной области (—оо<*<оо, г/>Я) при заданном граничном условии на линии у = = Н = const.

Рассмотрим одномерное течение идеального сжимаемого газа в канале переменного сечения F(x), нижняя граница которого образована осью симметрии и поверхностью тела, а верхняя — линией у = Н (рис. 1).

У

ylfcconst.

О 1 х

Рис. 1

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

X

put -f J ($v)H dx — const, (14)

—00

+ = c°nst» (15)

jr-C±, (16)

Poo /~у P—

где C+=——C_ = — --------------значение константы за скачком, кото-

pL рі

рое может быть найдено из обычных соотношений на скачке:

(ри)+ = (рй)-, (17)

(* + ї±Т *Ит + ї±ті)_. 04

(Р«2 + р)+ = (ры2 + Р)— (19)

Приведем уравнения (14) — (16) к безразмерному виду, относя все величины к параметрам во входном сечении канала Fco, Uoo, poo, р<*>. В ре-

12

зультате получим систему (для безразмерных величин сохраняем те же обозначения, что и для размерных):

ры/7 = 1,

где Г:

«2 +

р _

(7- 0 м2

1 +

(т- )) м

Р __ г

(20)

(21)

(22)

— безразмерная эффективная ширина ка-

1 — 1/^00

—со

нала с учетом расхода газа через линию у = Н.

Из системы (20)—>(22) можно получить уравнение для р:

2С.

(7 - О К

оТ-1

1.

(23)

Решение этого уравнения можно получить последовательными приближениями по формуле

1

Р*г+1

1'

(7 - 1)

1

7-1

Уравнение (23) имеет два корня, соответствующие дозвуковым и сверхзвуковым режимам течения в канале. Выбор нужного корня осуществляется логическим путем в подпрограмме расчета, устроенной следующим образом. Расчет ведется последовательно по сечениям слева направо. Если в минимальном сечении „ш канала возникла звуковая скорость, до этого сечения в решении выбирается «дозвуковой» корень. После чего из (20) — (22) можно найти и и р. Далее между критическим сечением и скачком уплотнения выбирается «сверхзвуковой» корень уравнения (23). В сечении, где расположен скачок, по известным параметрам до скачка и соотношениям (17) — (19) находятся параметры потока за скачком уплотнения. Далее расчеты ведутся по дозвуковой ветви соотношения (23), но с измененным значением константы С—

Положение скачка уплотнения в канале зависит от величины давления на выходе. При заданном давлении на выходе это положение находится итерационным путем.

Рассмотрим теперь дозвуковое течение потенциальной сжимаемой жидкости в полупространстве с граничным условием на линии у = Н:

Р2 Ухх •+• ?уу = 0.

¥* = «1. У = л>

?л

О, х2 + у2

оо.

Затухание возмущений позволяет воспользоваться для решения этой задачи преобразованием Фурье [11] и найти:

а)

У

Н=К‘const iu?i,Vuvi

S'S' ъм

цс,Рс> Ус

ffoо sOT7777y>^s

Начальное приближение

Vi

Ж

б)

= f p,Vjdx

1 щ\ V

I приближение-одномерное течение 5 полосе у-Н

Ус\рсРс

Ж пра^литение-РГРС'-ЦР-УС) -■

111

Метод малых Возмущений, при у

Релаксации v,

— Про &.ер к а точности

Решение с заданной точностью Рис. 2

V(X, Н)

оо

if

«1 (5) di

(24)

где интеграл берется в смысле главного значения.

При численном интегрировании (24) выбор шага интегрирования, длины интервала интегрирования и размеров малой окрестности особой точки осуществлялся путем сравнения с известным точным решением.

4. С использованием элементов, описанных выше, можно построить итерационную схему расчета. Рассмотрим ее на примере простейшей двухполосной схемы расчета обтекания симметричного профиля без угла атаки. Всю область течения разобьем линией у = Н = = const на две полосы* (рис. 2, а). В ближайшей к профилю полосе течение в первом приближении описывается уравнениями одномерной газовой динамики. Решение, полученное с помощью системы (15) — (20), припишем центральной линии полосы yc=Fc(x)/2 и обозначим

Для выбора Н воспользуемся соотношением (12).

ис, рс, Рс■ Давление на верхней границе полосы и на профиле, слабо отличающиеся от найденного решения, можно определить приближенно:

р {х, Н) - Рс + -дду {И — ус), Р(х, yJ = pc--jj-(yc — yw),

(25)

(26)

где для отыскания производной по у можно воспользоваться уравнением импульсов (3) и найденным решением для ис, рс:

Входящую в последнее соотношение скорость V находим из условия стыковки с внешней областью вдоль границы полосы у=Нл Условием совместности вдоль этой линии является равенство давлений и наклонов линий тока, которое достигается итерационным путем по следующей схеме (рис. 2, б).

Расчет начинается с некоторого нулевого приближения, полученного, например, по линейной дозвуковой теории, которое дает значение вертикальной скорости и1 и плотности р! вдоль линии у = Н. Рас-

ход через границу w(x)= ] pt г»! dx дает возможность рассчитать од-

номерное течение в полосе вблизи тела с учетом перетекания газа через нее. Найденные при этом параметры течения (первое приближение) приписываем центральной линии полосы ус> а к верхней границе переходим по приближенным формулам типа (25) (второе приближение). Выше линии у = Н лежит область линейного дозвукового течения, которое можно рассчитать известными методами (см. п. 3) и, в частности, найти новое значение скорости вдоль границы у=Н

по известному давлению Р\—РСЛ- ~ Ус)- Для этого восполь-

зуемся соотношением (24), предварительно определив по уравнению Бернулли величину возмущенной скорости ии соответствующей давлению Р1.

Новое значение скорости о, и плотности р! = рс + (Я— ус),

где дает новое значение расхода газа через гра-

°У 7 Рс оу

ницу полосы у = Н, и итерации повторяются до получения заданной точности. Для улучшения сходимости итераций вводится стандартная процедура релаксации. Обычно сходимость до 4-го знака после запятой достигается за 20—25 итераций.

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

Пример расчета по описанному способу представлен на рис. 3. Здесь рассмотрено обтекание профиля из дуг окружности с относительной толщиной 10%! под углом атаки а=0. При докритических скоростях наблюдается очень хорошее совпадение с экспериментом. Экспериментальные данные были получены в условиях, когда влияние сте-

X

•со

Про роль из дуг окружное та

Профвль вз дуг окружности *С*0,1; 14^=0, В

■ данный метод, Л =2 --— точное решение

• эксперимент Неъ= 1,1 щ6

Рис. 3

Параболический просриль %=0'06 0,303 ^

0,5 х/Ь

г /----данный метод Н =1,4

I------метод Мёрмана-Коула

-------- Джеймсона

• 4 эксперимент Ннехтеля В.еь=210в

Рис. 6

т

нок незначительно искажает обтекание профиля. (Затенение потока 3,75% при 20% перфорации стенок).

Сравнение с расчетом по методу установления Годунова (расчетная сетка 99x20 ячеек) при наличии на профиле обширной сверхзвуковой зоны и скачка уплотнения показано на рис. 4. Видно, что положение скачка совпало точно, а небольшие расхождения вблизи передней и задней критических точек обусловлены погрешностями примене-

ния одномерной теории в окрестности этих точек. Специальное рассмотрение этих областей или измельчение ширины полос позволит улучшить описание течения вблизи этих точек. Особо следует отметить, что решение по предложенной двухполосной схеме получено в 170 раз быстрее, чем по методу установления. Здесь же приводятся данные из работы [12], полученные при интегрировании уравнений Эйлера.

Рис. 5 демонстрирует сравнение расчетов по данному методу с экспериментом на толстом профиле [13], когда имело сильное вязкое взаимодействие в районе скачка уплотнения. Совпадение с учетом этого факта можно считать удовлетворительным.

Сравнение с экспериментом на тонком профиле в условиях слабых вязких эффектов (рис. 6) демонстрирует хорошее совпадение, особенно по положению скачка, и значительные преимущества предложенного метода перед существующими аналогичными методами.

ЛИТЕРАТУРА

1. Murman Е. М., С а 1 е Т. D. Calculation of plane steady transonic flows. — AIAA J„ 1971, N 9.

2. G о о r d j i a n P. М., Van Brisk R. Implicit calculations of transonic flows using monotone methods. — AIAA, p. 81-0331, 1981.

3. Годунов С. К., Забродин А. В., Прокопов Г. П. Разностная схема для двумерных нестационарных задач газовой динамики и расчет обтекания с отошедшей ударной волной. — Ж. вычисл. матем. и матем. физ., 1961, т. 1, № 6.

4. Tameson A., Schmidt W., Т u г k е 1 Е. Numerical solutions of the Euler equations by finite volume method using Runge—Kutta time-stepping schemes. — AIAA, p. 81-1259, 1981.

5. Atkins H. L., Hassan H. A. A new stream function formulation for the steady Euler equations. — AIAA J., 1985, N 5.

6. T a i Т. C. Application of the metrod of integral relations to transonic airfoil problem.—AIAA, p. 71-98, 1971.

7. Аэродинамика летательных аппаратов при трансзвуковых скоростях. Ч. I и II. — Обзор ОНТИ ЦАГИ, 1974, № 441, 442.

8. Методы расчета обтекания элементов летательных аппаратов при трансзвуковых скоростях. — Обзор ОНТИ ЦАГИ, 1980, № 585.

9. Кочин Н. Е., К и бель И. А., РозеН. В. Теоретическая гидромеханика. Т. I.—М.: Физматгиз, 1963.

10. В ан-Дайк М. Методы возмущений в механике жидкостей.—

М.: Мир, 1967.

11. Диткин В. А., Прудников А. П. Интегральные преобразования и операционное исчисление.—М.: Физматгиз, 1961.

12. Klopfer G. A. A Lagrangian method for the numerical solution of the Euler equations for transonic flows. — Proc. 6-th Int. Conf. of Numerical Meth. Fluid Dyn., 1978.

13. D e v i 11 Me. Т. B., Levy L. L. jr. De'iwert G. S. Transonic flow about a thick circular-arc airfoil. — AIAA J., 1976, N 5.

14. Численные методы в динамике жидкостей. — Сб. переводов с англ./Под ред. Белоцерковского С. М. и Шидловского В. П. — М.: Мир,

1981.

Рукопись поступила 23/1 1987 г.