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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XXIV 1993

№ 1

УДК 629.735.33.015.3.025.73 : 532.526

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

А. В. Волков, С. В. Ляпунов

Описан метод расчета вязкого трансзвукового обтекания профиля при больших числах Рейнольдса. Все поле течения разделяется на внутреннюю и внешнюю области. Во внешней части решаются уравнения, аппроксимирующие уравнения Эйлера, которые позволяют учесть изменение энтропии на скачках уплотнения и возникновение завихренности за ними. При решении внутренней задачи (расчет пограничного слоя) как на профиле, так и в следе используются интегральные методы. Учитывается кривизна вязкого следа и поперечный градиент давления в пограничном слое. Решение определяется с помощью итерационной процедуры. Учет вязко-невязкого взаимодействия осуществляется путем изменения граничных условий для внешней задачи распределением эквивалентных источников по профилю и в следе за ним. Показано, что учет завихренности за скачками уплотнения в настоящем методе позволяет улучшить согласование с экспериментальными данными как по распределению давления, так и по значению профильного сопротивления при незначительном увеличении времени расчета.

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

1. Решение внешней задачи. При решении внешней задачи расчета течения идеального сжимаемого газа необходимо решать систему уравнений Эйлера. Однако необходимое расчетное время в этом случае велико, что осложняет проведение массовых расчетов. Значительным

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

В данной работе, при решении внешней задачи использован подход, предложенный в [!]• Подход основан на применении преобразования Клебша. Вектор скорости потока V представляется в виде суммы потенциальной части уФ и части, обеспечивающей завихренность Ху;х:

V = уФ -|- Ауи. (1)

Функции Яиц могут быть выбраны так, чтобы выполнялись уравнения Эйлера:

V (рК) = 0 , (2)

— уравнение неразрывности;

VXV = -*VS (3)

— уравнение Крокко;

— уравнение Бернулли.

Здесь р — плотность газа, u) = rot V — вектор завихренности, а — скорость звука, х — отношение удельных теплоемкостей, S=^-~-у In ~ —

энтропия, р — давление, М» — число Маха набегающего потока. Все величины, имеющие размерность скорости, отнесены к скорости набегающего потока, давление и плотность к соответствующим величинам в набегающем потоке. Уравнение (3) является следствием уравнения Бернулли и уравнения сохранения импульса и заменяет последнее в системе уравнений Эйлера.

Однако, если использовать вместо скаляра К величину энтропии S,

а вместо вектора yjj. вектор V, т. е.

У = 7Ф-^,У, (5)

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

Так после подстановки преобразования (5) в уравнение неразрывности (2) имеем:

у(руф) = 0, (6)

где вводится модифицированная плотность

Л 1 р=р—V '

1 + Ш*

Величина плотности в каждой точке в соответствии с определе-

1

нием энтропии равна р = (М^а!)‘_ 1 е~\

Данное уравнение решается совместно с соотношением для скорости

К’( 1 - V 5) - VIV® 1 + | (5^ + !=!) = 10;, (8)

получаемого при подстановке уравнения Бернулли (4) ,в уравнение (5). Значения энтропии непосредственно за скачком уплотнения могут быть получены из соотношения:

5 =

1п(гттмг.-;-тт) + *1п(Ы

(9)

Далее, путем расчета линий тока и использования условия сохранения энтропии вдоль них, определяется поле энтропии за скачком уплотнения.

Подстановка величины скорости, определяемой уравнением (5), в уравнение Крокко (3), дает выражение:

УХ» = -|?5-;£-,1'Х»-т»'ХтнгХ у- (1°)

Два последних слагаемых являются дополнительными по сравнению с уравнением Крокко. Поскольку величины 5 и ш имеют третий порядок малости относительно интенсивности скачка уплотнения, то первое из дополнительных слагаемых имеет шестой порядок малости и им можно пренебречь. Порядок второго из дополнительных слагаемых формально совпадает с порядком основных членов уравнения. Однако в работе [1] показано, что он пропорционален произведению Бк, где ^ — кривизна линий тока, и поскольку величина кривизны обычно мала, данным слагаемым также можно пренебречь. Следовательно, представление скорости в виде (5) позволяет аппроксимировать систему уравнений Эйлера с точностью до малых слагаемых.

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

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

Алгоритм расчета во внешней области представляет следующее:

1. Во всем поле течения полагается 5 = 0.

2. Осуществляется несколько итераций решения уравнения (6) с помощью консервативной разностной схемы [2].

3. Определяются положения и интенсивности скачков уплотнения и вычисляется изменение энтропии вдоль них по формуле (9).

4. Рассчитываются формы линий тока путем численного интегрирования уравнений линий тока.

5. Определяются значения энтропии в узлах разностной сетки с помощью линейной интерполяции по известным значениям на линиях тока.

6. Вычисляются значения модифицированной плотности в центрах ячеек разностной сетки по формуле (7). Значения скорости при этом определяются из уравнения (8).

7. Пункты 2—6 повторяются до достижения сходимости.

Скорость сходимости итерационного процесса и время расчета при

этом являются практически такими же, как и при решении уравнения для потенциала скорости.

2. Решение внутренней задачи и расчет вязко-невязкого взаимодействия. Расчет пограничного слоя ведется интегральными методами. По заданной скорости, получаемой из решения внешней задачи, рассчитывается толщина вытеснения, толщина потери импульса и коэффициент поверхностного трения. Ламинарный участок течения рассчитывается по методике [3]. Расчет может вестись в режиме с фиксированной (заданной) точкой перехода и в режиме, когда точка перехода рассчитывается. Турбулентный участок потока рассчитывается по методу запаздывающей эжекции Грина [41. Интегрирование системы дифференциальных уравнений, описывающей развитие турбулентного пограничного слоя как ,на профиле так и в следе, проводится методом Рунге—Кутта—Мерсона 5-го порядка точности с автоматическим выбором шага.

Согласование решений внутренней и внешней областей осуществляется путем формулировки и решения эквивалентной невязкой задачи (ЭНЗ) в соответствии с работой [5]. Моделирование влияния вязкости на течение вне пограничного слоя осуществляется в. ЭНЗ путем изменения граничных условий на поверхности профиля и на осевой линии вязкого следа, которые ставятся следующим образом. На поверхности профиля задается величина нормальной скорости

что обеспечивает выполнение условия непротекания на линии, отстоящей на б* от поверхности профиля с точностью порядка (б*)2. Здесь и — тангенциальная скорость и р — плотность на поверхности профиля.

Вязкий след моделируется линией, на которой терпит разрыв как нормальная Ую, так и тангенциальная составляющая скорости.

Разрыв нормальной скорости обусловлен влиянием толщины вытеснения а*т следа и моделируется аналогичным способом:

(П)

где ри, — значение плотности газа на осевой линии следа.

Разрыв тангенциальной скорости обусловлен влиянием кривизны следа. Этот разрыв может быть выражен через интегральные характеристики пограничного слоя следующим образом [5]:

у(с/л = -^в/и*;+вв), (13)

где /Сш — кривизна следа, 0№ — толщина потери импульса в следе.

С целью уточнения получаемых результатов при решении ЭНЗ проводится учет градиента давления по нормали к поверхности в пограничном слое. В классической теории пограничного слоя этот градиент давления др/дп считается равным нулю. В настоящей статье, в уравнении для поперечной составляющей импульса сохраняются члены, пропорциональные кривизне поверхности, т. е. учитываются слагаемые более высокого порядка по сравнению с сохраняемыми в классической теории пограничного слоя. Учет этих слагаемых существенен в области задней кромки профиля. Уравнение сохранения нормальной к поверхности составляющей импульса имеет вид:

(14)

где ко — кривизна поверхности профиля. Интегрирование уравнения (14), в случае вязкой задачи, дает:

Рь-Ре = Ьо?.Щ(~-Ъ*-Ъ + Ъ), (15)

где ре и рь — соответственно значения давления на внешней границе пограничного слоя и на поверхности профиля в реальном вязком течении, 6 — толщина пограничного слоя на профиле, 0 — толщина потери импульса.

В ЭНЗ при интегрировании (14), с учетом малого изменения продольной компоненты скорости, получаем:

Р — Ре = Ь й?еиЧ. (16)

Приравнивание значений давления на внешней границе пограничного слоя в (15) и (16) позволяет получить выражение для разности давлений на поверхности профиля в ЭНЗ и реальном вязком течении в сле-

дующем виде:

Р-рь = коРеи1(Ъ* + в). (17)

Как показано в [5], детальный анализ членов более высокого порядка в (14), приводит к замене величины кривизны поверхности профиля ко на величину кривизны контура вытеснения к*. Таким образом, получаемые при решении ЭНЗ значения давления на поверхности профиля и в следе отличаются от действительных значений давления в вязком потоке на величину, прямо пропорциональную кривизне поверхности вытеснения и сумме толщин вытеснения и потери импульса. В реализованном методе, для учета указанной поправки, после решения ЭНЗ осуществлялась коррекция локальной скорости по формуле:

и„ = и[\+11*(Ъ* + Ь)\. (18)

Для учета взаимодействия сильного скачка уплотнения с пограничным слоем применялась процедура, предложенная в работе [6], в

которой по существу моделируется внезапное утолщение пограничного слоя за ударной волной.

Расчет профильного сопротивления производился путем интегрирования по поверхности профиля коэффициентов давления и трения. Для оценки величины волнового сопротивления использовалась формула работы [7].

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

3. Результаты расчетов. Существенным элементом настоящей методики расчета обтекания профиля является учет завихренности за скачками уплотнения. Исследования данного способа учета завихренности на примерах течений идеального газа проведены в работе [1]. На рис. 1 демонстрируется сравнение результатов расчета обтекания профиля 1\АСА-0012 идеальным газом при М(х, = 0,8, а=1,25. Здесь приведены результаты, полученные путем решения уравнения для потенциала, с помощью настоящей методики, а также путем решения уравнений Эйлера [8]. Видно, что расчет по данной методике хорошо согласуется с решением уравнений Эйлера как по распределенным, так и по интегральным характеристикам. Результат, полученный в рамках потенциальной модели, является существенно неверным.

Далее приводятся расчеты, полученные с учетом влияния вязкости. На рис. 2—5 представлены результаты расчетов обтекания различных профилей по настоящей методике, а также полученные с использованием решения уравнения для потенциала скорости с помощью неконсервативной разностной схемы. Эти результаты сравниваются с экспериментальными данными и с решением уравнений Навье—Стокса (уравнений Рейнольдса, рис. 2).

Рис. 2 представляет результаты расчетов обтекания сверхкритичес-кого профиля 05МА-523 с острой задней кромкой, с = 11%, на режиме трансзвукового обтекания Мс = 0,765, су = 0,7, Ке = 2 млн. Сравниваются результаты, полученные по данной методике, с результатами решения уравнений Рейнольдса с алгебраической моделью турбулентности

Рис. I

Прифаль В5МА-523;М=0,7£5, Я*=2Ю‘\с,*0,71

о решение уравнении Набье-Стокса С?]

0,70 0,0159

удк------решение с использованием 1,4-2 0,0217

' неконсерВатибной схемы

— настоящий расчет 0,85 0,0147

Рис. 2

Профиль ЯАЕ-2822,М-0,723-6,5-10 ;су^0,74

„оОООО

ко

ОС Сх

эксперимент [10] 2,32 0,0127

решение с использованием 2,75 0,0141 неконсерВатиВной схемы

настоящий расчет ' 23В 0,0133

Рис. 3

а сх

эксперимент [//7] 3,86(геом.) 0,0127

■ решение с использованием 3,03 0,0151

неконсербати&ной. схемы

■ настоящий расчет 2,88 0,0140

Рис. 4

О,В

ОМ

0,2

Профиль N АСА-0012 М-0,7 ,В.е=9-10в

Г/'

(/ о эксперимент []0\

'/ — — решение с использованием" неконсерВа-т ив но и схемы настоящий расчет

0,02

ОМ сх

Рис. 5

[9]. Интегрирование уравнений Рейнольдса проводилось на криволинейных сетках (204x55 узлов) двумя методами: явным методом Маккормака и неявным методом Бима—Уорминга, причем расчеты по обеим методикам близки и на рисунке нанесены крестиками. При расчетах положение перехода ламинарного пограничного слоя в турбулентный фиксировалось соответственно на верхней и нижней поверхностях при х = 35 и 18%. Видно, что расчет с учетом изменения энтропии за скачком уплотнения более близок к результатам работы [9] как по эпюре давления, так и по величине угла атаки и профильного сопротивления.

На рис. 3 даны результаты расчетов обтекания сверхкритического профиля РАЕ-2822 при Моо = 0,729, ы = 2,92, Не = 6,5-10в. Экспериментальные результаты, соответствующие этому случаю, были использованы в качестве одного из контрольных примеров в [10]. В эксперименте

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

На рис. 4 и 5 приводятся сравнения результатов расчетов с экспериментальными данными для профиля NACA-0012 при М<х> = 0,7, Re = = 9-10б [10]. Переход был фиксирован при х = 5%. Расчет по предлагаемой методике улучшает согласование с экспериментальными данными в области скачка уплотнения, а рассчитанная поляра хорошо согласуется с экспериментальными данными (рис. 5). Более высокое сопротивление, получаемое при использовании неконсервативной схемы, объясняется присутствием схемных источников на скачке уплотнения.

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

ЛИТЕРАТУРА

1. Ляпунов С. В. Ускоренный метод решения уравнений Эйлера в задаче о трансзвуковом обтекании профиля//Математическое моделирование.— 1991, № 3.

2. Jameson A. Acceleration of transonic potential flow calculations on arbitrary meshes by the multiple grid method//AIAA Paper. — 1979,

N 1458.

3. Cohen С. B., Reshotko E. Similar solution for the compressible laminar boundary layer...//NACA Rep. — 1294, (1956).

4. Green J. E„ Weeks D. J„ Brooman J. W. F. Prediction of turbulent boundary layers and wakes in compressible flow by a lag — entrapment method//RAE Teh. Rep. —72231 (ARC—RM 3791) (1977).

5. Lock R. C., Firm in М. C. Survey of techniques for estimating viscous effects in external aerodynamics//Proceedings of IMA conference on numerical methods in aeronautical fluid dynamics, 30 march—1 april 1981, edited by P. Roe, Academic Press (1983).

6 Desai S. S., Rangarajan R. Viscous transonic flow over aerofoils using transonic full potential equation in system of cartesian coordinates//AIAA Paper.— 1987, N 411.

7. Боксер В. Д., Серебрийский Я. М. Приближенный метод определения волнового сопротивления при наличии местной сверхзвуковой зоны//Ученые записки ЦАГИ. — 1978. Т. 9, № 5.

8. AGARD-AR-211. Test cases for inviscid flow field methods.— 1985.

9. Visbal M. R., S h a n g J. S. Comparative study between two navier-stokes algorithms for transonic airfoils//AIAA Paper.— 1985, N 480.

10. Holst T. L. Viscous transonic airfoil workshop compendium of results//AIAA Paper. — 1987, N 1460.

Рукопись поступила 27/V 1991 г.