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

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

Том XXVIII 199 7 №1

УДК 533.6.011.35

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

; Ю. П. Нуштаев, В. А. Павлов

Рассматривается алгоритм численного решения уравнения Кармана для задачи трансзвукового изоэнтропического обтекания профиля, приводятся примеры расчетов. Расчеты проведены методом малых возмущений для симметричного профиля НАСА64А00б с относительной толщиной 6%.

Уравнение Кармана выводится из точных уравнений газовой динамики путем асимптотических разложений в предположении о малости возмущений, вносимых в околозвуковой невязкий поток [1]. Решение задачи об обтекании профиля на основе уравнения Кармана помимо самостоятельного научного интереса может оказаться очень полезным для нахождения исходного приближения при решении нестационарной задачи о колебаниях тонкого профиля [2].

Уравнение Кармана для потенциала возмущения может быть записано в безразмерном виде следующим образом:

[^ - (у + 1)М^фх]фхс + = 0, (1)

1-М2

где К - ——— параметр подобия Кармана, т — относительная тол’ '

щина профиля, координата х относится к величине полухорды профиля с, координата у — к величине ст1^3, а потенциал ф — к величине

с^\ Параметр т выбирается в зависимости от числа Маха [3].

Пусть профиль расположен вдоль оси х при |х| < 1. Если форма поверхности профиля задается в виде функции F(x), то в рамках теории малых возмущений граничное условие на оси у = 0 при |х| <1 записывается в виде .

фу(х, 0) = i[P-4 (2)

где а — угол атаки профиля.

Решение для профиля с острой задней кромкой должно удовлетворять условию Жуковского — Чаплыгина, которое в рамках теории малых возмущений сводится к требованию непрерывности величин ф* и фу в следе, т. е. при у = О И X > 1.

Граничное условие на большом удалении от профиля можно получить с использованием известного асимптотического выражения для потенциала возмущений

Г Су

ф = Пт — агс1ап —

Х-+СО 271 V х

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

ф(х, =

ф(дс, -оо) = ^,

4

г (3)

ф(~°°, ^) = у> •

ф(«=, у<0) = Г,

ф(оо, у > 0) = 0.

В данной работе используется преобразование координат, которое отображает бесконечную физическую область течения на внутренность квадрата. При этом неравномерная в физической плоскости ХУ сетка с концентрацией узлов по оси О У в окрестности верхней и нижней поверхностей профиля, а по оси ОХ в окрестности передней и Задней кромок профиля преобразуется в равномерную сетку в плоскости (£, ту). Преобразование имеет вид

| (еа1х _х) при 0^х^1,

| = ~7--—-----г(е“1Х-1) при -1йх<0,

И-1)

% = 1 — (1 — при х>1, ^

? = (1-?о)ев3^-1 при х<-1,

Т| = 1 — е°2У при у > 0, т| = еа7У -1 при у < 0,

здесь ~~ положение кромок профиля в новой системе координат, а\, а2> аз — параметры, влияющие на степень концентрации узлов в окрестности профиля. В новых координатах (£, ту) уравнение Кармана принимает вид

и*-(у + 1)м^**б)(е*ф*)* +(чуФп)чт1у =0-

(5)

Аналогичным образом можно записать граничные условия в новых координатах. В качестве конечно-разностного аналога данного уравнения рассматривается смешанная эллилтико-гиперболическая система уравнений, которая решается итерационным методом релаксации. Тип уравнения определяется величиной скорости ф?. Для дозвукового и сверхзвукового случаев используются различные разностные аппроксимации производных ф4, ф?5. Если в рассматриваемой точке (х, у) течение дозвуковое, то производные ф?, фк на п+1-й итерации записываются следующим образом: '

25,(т + 1)М^Л5*”+1 =(ф-Ф"_і,у)(5х),(ї + 1)М*, (6)

(і5)2(^ф5)”*‘=- »*7)(5,ц - (ф"+/■ <7>

га» (?,Ц=і((5,)/11+(5*)/), Зд“ь » -

параметр релаксации, в случае дозвукового течения 1 < со < 2. Если в точке (/, у) течение сверхзвуковое, то

2«(Ф«СЖФ”>-*7-2,у.

(8)

(ч)2^)”1 -(іи -(ф?-й -ф?-2,у)(4<),л- (9)

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

Ґ

-±^Д. (1 - Щ)С* ~ -ФГ/ -1 1 ~

(Д§) \ . 1 24 ,} ® V ®

ф” .

(5*)а(2Фї+/ - Ф?у - ф*+1/ - и)+(ф"_2,у - фйу )]'

+—-~т{чу) ■ (Ал)2 " '

= 0,

где С” = К- (у +1)М2 (|х)г (ф"+1)у а = 0 ПРИ С,П > 0

(10)

и

Ег- — 1 при С. <0, і — 1, 2, ..., іщзх, у — 1, 2, ..., Лпах- На линии т] — 0

и - Уо) вводятся два слоя расчетных точек / = ;'0+0 и У = у0 - 0, в которых производная ф,,л записывается с учетом граничного условия на теле и в следе при

(

j=jo±0 Лт1

Ф»,Уо±1 Ф»,Уо _ АЛ

(П)

при

So < |s|: Ал2((л)3,Фп) = (%)у0+1(ф/,Уо+1 " *ЛЛ+о) -

Уо=У+о

)y0-lt*/,yo+° + Г ~ ФлУо-i)’

Ат12(луФг,)

У=Уо—0

= MJ(j+l(bij0+i + Г " Ф^Уо-о) -

Выписанная система алгебраических уравнений на каждой итерации решается последовательно для каждого вертикального луча % = const. Подход осуществляется от % = -1 до 4 = 1- Решение уравнений для ^ = const осуществляется методом встречных прогонок. Для \ = const решаемая система имеет вид

Aft j_i - С уф у + 5уф j+i - -FI,

(12)

где коэффициенты А], Ср В], Р/ выражаются через решение, полученное на предыдущей итерации. Формулы метода встречных прогонок имеют вид

где

Ф j — ау+1Фу+1 + Py+i>

В; А ;В 1 + Fj

а/+1 = г--—Г> BJ+1 = г~^—а для 2 < у ^ у0 - О,

cj - ajAj cj ~a-iAj

(13)

(14)

Фу+1 - £у+1Фу + Лу+1!

А,- В.-т\1+1+р1

здесь |у = Ц] = С--г\-+1В- ДЛЯ •/'шах - О-

Из граничных условий следует, что а.2 = % утах =0, (32 = -^-Г, Г

т1утах=—. После НаХОЖДеНИЯ ПрОГОНОЧНЫХ КОЭффИЦИеНТОВ (Ху, ру,

•Пу, £у находится решение фу0+о- Для т. е. в области перед

профилем, фу0+о находится из системы уравнений

Фуо-о = аЛ+1 + Рл +1> ф/о+о = 1у„ + 1Фу0+о + 11/о+ь

Фуо+о = Фуо-о-

Отсюда следует, что

А Руо+1 + схУо+1Г|уо+1

Ф Уо+о = —;-----------------•

1~аУо+111л+1

Для области -£о ^ ^ из системы уравнений лУофУо-1 ~сл-оФуо-о = -770-0>

Фуо-о = ау0Фу0-0 + Руо>

получаем

Фуо-0 =

^/0-0 + 4/0 Рур

................. £Уо~о . 4/»аУо

а система уравнений

“^уо+офуо+о “ ^УоФуо+1 = ~-^УЬ+0> Фуо+1 = т1у0+1Фу0+0 + т1у0+1

дает решение

х /Уо+0 +

Фуо+0 =-----------

+1

(16)

(17)

(18)

(19)

(20)

^Уо+о -®УоГ1Уо+1 В области следа \ < §о решается система уравнений

вк (Фуо-о + Фуо+1)_ сУо+о Фуо+0 = -^у0+о>

Вк (ФУо-1 + ФУо+1) “ СУо-° ФУо-о = _^Уо-° >

ФУо-1 =аУоФу'о +Руо’

Фуо+1 = ^Уо+1^Уо+0 + т1Уо+1- (21)

Ее решение имеет вид

С/о+0 + -8у0 у0 (т1Уо+1 + Ру'о)+ РуоаУо (^Уо-° ~ ^Уо+°)

Фуо+о =

■ СУо-оСУо+о ~ гУо1аУоСУо+0 + ^Уо+^Уо-о)

сУо+°ФУо+о +(^Уо-о “ *Уо+о)‘ Фуо-о =---------

С,п-

Уо-о

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

Таким образом, на каждой итерации находится решение последовательно для каждого вертикального луча § = const. После расчета для последнего луча, проходящего через профиль, вычисляется величина циркуляции, которая затем используется для формирования граничных условий в следе и на границах расчетной области для следующей итерации. При расчете циркуляции используется формула Г = ф/ j0_o-

-ф,(е j+0, где 1К — номер последнего узла сетки на задней кромке профиля. Эта формула представляет собой разностную аппроксимацию условия Жуковского — Чаплыгина. В каждой точке расчетной сетки на каждой итерации производится анализ величины скорости ф^ и выбор

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

На основе изложенного алгоритма была составлена программа расчета на ЭВМ. При расчетах итерационный процесс считался сошедшимся, если погрешность решения была

В расчетах значение а\ принималось равным 0,92, а коэффициент релаксации в дозвуковой области принимался равным 1,96. При этом итерационный процесс сходился за 200—300 итераций в зависимости от числа М.

На рис. 1 приведен пример расчета стационарного трансзвукового течения около профиля МАСА 64А006. Приведено распределение коэффициента давления ср{х) на поверхности профиля при числе М.,, =0,875 и угле атаки а = 0. Коэффициент давления рассчитывался по формуле

(23)

2

NACA 6ЧА006^0,675; а ~0

•2тЗф*. (24)

~tz ijr Oj ~~ts JC

Рис. 1

Отметим, что вследствие наличия у данного профиля тупой передней кромки стационарное решение сильно зависит от расположения передней кромки относительно узлов расчетной сетки, т. е. от величины £о- Приведенные зависимости показывают, что величина §0 влияет не только на течение в окрестности передней кромки, но и на положение скач-

ка уплотнения. При = 0,62 узел расчетной сетки располагается приблизительно на расстоянии 1% хорды профиля, а при |0 ~ 0,6085 приблизительно на расстоянии 0,5%. Сравнение результатов расчетов для данного профиля, полученных данным методом, с результатами, представленными в работах [4], [5], показало хорошее их согласование.

Как показано в работах [3], [6], нестационарные аэродинамические характеристики профилей в трансзвуковом потоке претерпевают сильное изменение. Расчеты, проведенные на примере профиля КАСА 64А006, показали, что существенно нелинейными являются и стационарные характеристики [7]. Это также видно на графиках, представленных на рис. 2—4, полученных в результате решения стационарного уравнения Кармана. Эти графики показывают, что нелинейность изменения коэффициента подъемной силы по углу атаки проявляется в сравнительно небольшой зоне чисел Мш от 0,8 до 0,92, когда на профиле возникает смешанное течение. Относительная величина

приращения коэффициента су при изменении угла атаки сильно зависит от двух параметров: самого угла атаки и числа М*,. В том случае, если обтекание профиля дозвуковое (Мш = 0,8), а также, если скачки на обеих поверхностях стабилизированы, а обтекание в основном сверхзвуковое (М„ = 0,92), то графики су(а) близки к линейным. Ко-

Рис. 2

эффициент су достигает своего максимального значения при значениях числа М тем меньших, чем больше угол атаки. Соответствующее значение модифицированного параметра Кармана 1-М2

Ка = —--------:—=_ для данного экс-

(то -I- 2са)з тремального значения су, приведенное на врезке к рис. 2, при изменении угла атаки остается практически постоянным.

Вторую интересную особенность трансзвукового обтекания можно усмотреть из графиков тг(Мм) на рис. 4. Эти графики показывают, что для всего обследованного диапазона углов атаки существует область значений числа М, при которых центр давления находится позади середины профиля. Эта область зависит от угла атаки и соответствует тем значениям числа Мм, при которых наблюдается существенное возрастание несущих свойств. Кроме того, существенное влияние на коэффициент mz оказывает сильное увеличение несущих свойств задней части профиля на наветренной стороне (рис. 3, d).

Резюмируя полученные результаты, можно сделать следующие выводы:

область нелинейного изменения коэффициентов су(а) и mz(a) по

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

коэффициент су при стационарном обтекании достигает максимального значения при числах М„ тем меньших, чем больше угол атаки. Соответствующие значения модифицированного параметра Карма-1-М2

на Ка =---------2__ остаются примерно постоянными;

(то + 2са)з

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

Рассмотренный алгоритм был разработан Ю. П. Нуштаевым в 1978 г. Постановка задачи обсуждалась им с Ю. Б. Лифшицем. Численные параметрические расчеты были проведены позднее.

ЛИТЕРАТУРА

1. Karman Т. The similarity law of transonic flow // J. of Mathematics

and Physics.—1947. Vol. 26, N 3.

2. Нушхаев Ю. П. Нелинейные задачи аэродинамики и аэроупругости крыла бесконечного размаха с элероном в трансзвуковом потоке. Диссертация на соискание ученой степени кандидата физико-математических наук.—1981.

3. Ballhaus W. F. and Georgian P. M. Implicit finite difference computations of unsteady transonic flow about airfoils, including the treatment of irregular Shock wave motion // AIAA Paper 77-205.—1977.

4. Magnus R. J. and Yoshihara H. Calculation of transonic flow over an oscillating airfoil // AIAA Paper 75-98.—1975.

5. Causton M. and Angelini J. J. Solution of nonsteady two-dimensional transonic small disturbances potential flow equation // Communication presented on Symposium.—1975.

6. Лифшиц Ю. Б. К теории трансзвуковых течений около профиля // Ученые записки ЦАГИ.—1973. т. IV, №5.

7. Tijeman Н. On the motion of shock waves on an airfoil with oscillating flap in two-dimensional transonic flow // NLR TR-75038U.—1975.

Рукопись поступила 2/VIII1995