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

_____________УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА

Том 154, кн. 1 Физико-математические пауки

2012

УДК 533.6.011

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

Е.М. Котляр, Д. В. Маклаков

Аннотация

Решена задача о безотрывном обтекании профиля заданной формы дозвуковым потоком газа. Задача сведена к нелинейному интегральному уравнению типа уравнения Вилла, которое после дискретизации решено численно методом Ныотопа. Показано, что распределения скоростей, рассчитанные по модели газа Чаплыгина, по формуле Кармана Дзяпа и с помощью CFD пакета Fluent., в дозвуковом диапазоне чисел Маха очепь близки.

Ключевые слова: дозвуковые течения, уравнения Чаплыгина, газ Чаплыгина, уравнение Вилла, дискретизация, метод Ныотопа.

Введение

Один из приближенных методов учета сжимаемости базируется на линейных уравнениях газовой динамики, записанных в плоскости годографа скорости уравнениях Чаплыгина [1]:

1 дер 1 дф

А дв рдХ’

(1)

1 д<р (1 ( 1 \ дф Х~дХ ~ Ж ^) дё’

где ер - потенциал течения, ф - функция тока, А и в - модуль и угол наклона вектора приведенной скорости (скорости, отнесенной к критической скорости звука а* ), р безразмерная плотность, отнесенная к плотности ра в точке торможения потока. При этом если Ах и Ау - компоненты вектора скорости А, то

ж дх ’ у ду’ х рду ’ у рдх' ( ’

Уравнения (1) являются следствием соотношения

которое вытекает из равенств (2).

С.А. Чаплыгин [1] предложил два метода решения системы (1): точный и приближенный. Точный метод может быть применен к достаточно широкому классу задач о течениях сжимаемой жидкости, ограниченных свободными поверхностями и полигональными твердыми стенками. Обзор работ, выполненных в этом направлении, можно найти в монографиях [2, 3].

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

задачи в приближении газа Чаплыгина даны Г.Г. Тумашевым [4. 5] (см. также [6]). Вудсом ( [7. п. 6]).

Прямые задачи, когда форма профиля задается, являются более сложными, чем обратные. Здесь явного аналитического представления решения получить ие удается. либо строится приближенное аналитическое решение [8]. либо предлагается итерационная процедура, сходимость которой исследовать весьма затруднительно. В частности, такая процедура представлена в [9]. но авторы ограничились расчетом только симметричного профиля NACA 0012.

В настоящей работе, используя приближенный метод Чаплыгина, мы сводим решение задачи об обтекании произвольного профиля дозвуковым потоком газа к нелинейному уравнению типа уравнения Вилла [2]. После дискретизации это уравнение решается методом Ныотона. Численные эксперименты показали, что метод сходится для профилей любой формы всего за несколько итераций. Полученные распределения скоростей сравнивались с распределениями, найденными в CFD пакете Fluent, а также с хорошо известной формулой Кармана Дзяна [8]. Установлено, что в дозвуковом диапазоне чисел Маха все три подхода дают весьма близкие результаты.

1. Модель газа Чаплыгина

Приближенный метод С.А. Чаплыгина состоит в аппроксимации адиабатической зависимости

Р = Г (3)

в плоскости (1 /р,р) прямой линией: вместо (3) используется уравнение

_ A , ,

р = -— + В, (4)

Р

где к показатель адиабаты (для воздуха к = 1.4), р = р/ра, р = р/ро безразмерные давление и плотность, отнесенные к давлению ро и плотноети ро в точке торможения потока, A > 0 и B постоянные. Газ, в котором связь между давлением и плотностью задается уравнением (4), называют газом Чаплыгина.

Уравнение Бернулли

V2 } dp

-г- + — = 0,

2Р

Р 0

где V - скорость, с учетом равенств

2 кро а2 к +1

°о = > ~ = о

Ро а( 2

(ао - скорость звука в точке торможения, а* - критическая скорость звука) приводится к виду

А= + 1±1/1^ = 0. (5)

к .] Р ар 1

Подставив (4) в (5), после интегрирования получим

Р= /л Л с2 = о л ^ > 0, (б)

д/1 + 4с2 А2 2А (к + 1)

причем результат интегрирования не зависит от величины параметра В. Если теперь подставить (6) в систему (1) и сделать замену

чч)_ ехР(й') м

1 — с2 ехр (25') ’ (,)

то вместо (1) получим условия Коши Римана

скр = дф_ д^ = _дф_

дв дБ’ дБ дв' { '

При этом комплексный потенциал т = ^ + гф будет аналитической функцией переменной х = Б — гв. Переход в физическую плоскость осуществляется по формуле

dz = exp (—х) dw — c2exp (х) dw (9)

(черта означает комплексное сопряжение), которая следует из (2).

Введем новую комплексную переменную zi = xi + iyi, связанную с w и х соотношением dzi = exp(—х) dw. Область изменения этой переменной называют областью фиктивного течения несжимаемой жидкости. Тогда dw/dzi = exp х -комплексно-сопряженная скорость фиктивного течения, Л = exp S - модуль этой скорости. Из формул (6) и (7) найдем

Л 2Л _ _ 1 — с2Л2

~1-с2Л2’ _ 1 + а/1 + 4с2 Л2 ’ Р~1 + С2Л2' 1 1

Предположим теперь, что каким-то способом построена область течения несжимаемой жидкости со скоростью на бесконечности

2А^

1 + л/1 + 4с2

и найден комплексный потенциал т течения в этой области. Тогда с помощью формулы (9) можно найти соответствующую область течения газа Чаплыгина и по первой и третьей из формул (10) построить поля скоростей и плотностей для

в

обоих потоках одинаков.

Выбор угла наклона прямой (4) или, что равнозначно, параметра с2 можно осуществлять различными способами (см., например, [10, с. 61]). С.А. Чаплыгин заменял адиабатическую зависимость р = р(1/р) касательной прямой, проведенной в

точке р = 1, р = 1 торможения потока. В этом случае А = к и с2 = 1/[2(к+1)]. Па-с2

зависимости

Л2 (к — 1) т1/(к-1)

Р=

1 -

к +1

(11)

между приведенной скоростью Л и безразмерной плотностью р приближенной заЛ

В частности, Г.Ю. Степанов (см. [11, п. 24]) предложил универсальное значение с2 = 0.296. Расчеты показали, что при этом относительная погрешность формулы (6) по отношению к (11) в диапазоне 0 < Л < 0.89 составляет всего 2.3%. Если с = 0

2. Постановка задачи и сведение ее к нелинейному интегральному уравнению типа уравнения Вилла

Профиль обтекается дозвуковым потоком газа в плоскости z = x + iy. На бесконечности заданы скорость Уто, давление рто, плотноеть рто и угол наклона вектора скорости (угол атаки) аа (рис. 1).

Определим число Маха Мто па бесконечности. Для адиабатического течения квадрат скорости звука а^ = крто/рто. Отсюда найдем

М2 =Yk = l p^Yk

оо о *

а5о к Р<^

Пусть У - скорость потока, а* - критическая скорость звука. Будем искать поле приведенных скоростей А = У/а* = Ax+iAy. При заданных Уто , рто и рто параметр а* является физической константой, которую легко найти по известному числу Маха Мто. В дальнейшем будем считать, что Мто задано. Из адиабатического соотношения находим

К + 1

Лоо = Моо\/2+(« —1 )М^

и отсюда определяем а* = Уто / Ато.

Форма профиля задана уравнением

в = F (sr), 0 < sr < L, (12)

где в - угол наклона касательной к профилю, sr - дуговая координата, отсчитываемая от острой кромки A против часовой стрелки, L - периметр профиля. На рис. 2 показана функция F (sr) для профиля NACA 2411.

Точка B является точкой торможения потока, еп - внешний по отношению к профилю угол в острой кромке A (1 < е < 2). Отметим, что резкие изменения

в

профиля действительно имеют место.

Перейдем в область фиктивного течения несжимаемой жидкости zi = xi + iy 1. Пусть Л = Лвгв - вектор скорости фиктивного течения, a w — комплексный потен-

dw _-а

циал в этой области. Тогда —— = Ле комплексно-сопряженная скорость фпк-

dz1

тивиого течения. Для нахождения поля приведенных скоростей реального течения будем использовать формулы (10), где параметр с2 выберем равным универсальному значению Г.Ю. Степанова с2 = 0.296.

в

б)

в)

І 1 1 1 '201і 1 1 12.02 і

Рис. 2. Функция 0 = Е (вг) для профиля ИАОА 2411 (о); поведение функции 0 = Е (вг) у задней кромки верхней части профиля (б) и пижпей части профиля (с)

Рис. 3. Область фиктивного течения и каноническая область

Будем искать конформное отображение *1 = *1 (£) внешности единичной окружности в плоскости £ = £ + на внешность профиля в области фиктивного течения с условиями *1 (то) = то, Иш7^ 0+ *1 (е'1'7) = 0 (рис. 3).

В плоскости фиктивного течения профиль будет разомкнутым: точке А в плоскости ^соответствуют две точки А+ и А-, с которых сходят две конгруэнтные линии тока, простирающиеся до бесконечности (как в схеме Ву [2]). В самом деле, функция й*1/Л является аналитической, поэтому все интегралы от этой функции по замкнутым контурам, охватывающим единичный круг в плоскости £ будут равны между собой. Поэтому линия тока в плоскости £, отходящая от точки А, превратится в две конгруэнтные линии в плоскости *1. Так как вдоль линий тока выполняется равенство ( [13, с. 25])

Лг = (1 — с2Л2) Лг

(13)

а профиль в физической плоскости * - замкнутый, то разомкнутость в плоскости * исчезает, И конгруэнтные ЛИНИИ тока В ПЛОСКОСТИ *1 переходят в одну линию в плоскости

Изменение величины 0 по всему контуру профиля А0 = ^ (£) — ^ (0). Тогда

имеет место соотношение £ ^*1

3-^.

п

Производную —— будем искать в виде а£

І- 1

г-1

6.5

где Ф (t) — анадитическая функция, не имеющая особенностей в канонической об-

ласти, U0 =

dzi ,

— Ы]

постоянная, имеющая размерность длины.

На границе канонической области введем обозначения

Ф (е®7) = ^ (7) + гт (7), 0 < 7 < 2п ,

где 7 — полярный угол. В предлагаемом методе искомой функцией является функция т (7), для определения которой мы выведем нелинейное интегральное уравнение. Если предположить, что т (7) - известна, то тогда известна реальная часть аналитической функции —гФ (е®7) = т (7) — г^ (7). Последнюю можно восстановить с помощью формулы Шварца, а затем найти и

2п

1C ег

+1 - t

dY:

(15)

причем здесь учтено, что U0 =

dz

S(o°'

и следовательно, ReФ(то) = 0. Из (15)

найдем, что связь между функциями т (7) и ц (7) дается интегралом Гильберта:

2П

А‘(/3) = J г (7) ctg d7.

(16)

Теперь определим, как через функцию т (7) выражаются функции в (7) и вг (7). Если в выражении (14) перейти на границу канонической области, то с учетом того, что

dY dt dY

dt

получим

dzi / Y

-Ц-Ч3тг

Отсюда выводим

■-i

exp

* (е - !) + * +7) -М(7) - *т (7)

dzi п Y .

= al'g d7 = £ 2 + 2 ( “ Т (7)'

Комплексный потенциал в канонической области имеет внд

(17)

w (t) = сра lat Н—-—|- 2г sin a In t а

ной силы), не совпадающий с истинным углом атаки аа, ^0 _ постоянная, имеющая размерность потенциала скорости. Тогда

dw

dt

^0e

t-1

~

te

^(n+2a)

dw dw/dt

Найдем (po и aa. Для этого запишем равенство —— = ————

dzi dzi /dt

Ле

(18)

на

бесконечности, в выражениях (14), (15), (18) также устремим і к бесконечности. В результате получим

2п

(20)

При переходе на границу канонической области выражение (18) принимает вид

Воспользуемся тем, что (14), (19) и (21), получаем

(в*7) =— 4^о*є

іиі Іш/іИ

І2і І2і/ІІ

'Т \ 2

1

2'

(21)

= Л (7), тогда, учитывая соотношения

2-е

сое ( — — а

Из выражения (13) следует, что

7

и0 [ (2вш^У е-^с1г{-- 4с2Л2

3-е

I — — о. ] е

.М(7)

^7

. (22)

Подставив (17) и (22) в основное уравнение (12), задающее форму профиля, найдем нелинейное краевое условие, которое должно выполняться для аналитической функции Ф (£) на параметрической окружности:

ио

2 вігі %

г-1

,-М(7)

^7 —

- 4с2Л^ У (2эт ^ соэ2 - а) ем(7)

.

Учитывая, что периметр данного профиля известен и равен Ь, получим условие для определения

Ь = ио

2п

2 йіп %

-1

-М(7)

^7 —

22

2п

- 4с2Л

2 єіп -2

3-е

I — — а ) с

.М(7)

^7

.

Краевое условие (23) превращается в нелинейное интегральное уравнение типа уравнения Вилла [2], если учесть, что связь меду функциями т (7) и /и (7) дается линейным сингулярным оператором (16) с ядром Гильберта. Тогда (24) нелинейное функциональное соотношение для определения параметра . Кроме того, мы имеем интегральное соотношение (20) для отыскания теоретического угла атаки а. Таким образом, задача сводится к решению системы (16), (20), (23), (24) относительно функции т (7) и параметров Щ, а. Для дискретизации этой системы применялся метод, предложенный в [13] для расчетов кавитационного обтекания профилей.

а

3. Дискретизация системы уравнений

На отрезок [0, 2п] нанесем сетку узлов 7* , * = 1,..., п; 71 =0; 7„ = 2п. Нам

потребуется дискретные аналоги линейного функционала

2п

Пусть Б* (7), г = 1,..., п, - система фундаментальных кубических сплайнов, построенная на сетке 7* с граничными условиями типа “по^а-кпо^’ [14]. Фундаментальный кубический сплайн Б* (7) - это сплайн, значения которого во всех узлах равны 0, кроме узла с номером г, в котором Б* (7*) = 1. Системой фундаментальных сплайнов будем пользоваться для приближенного вычисления интегралов. Имеем

Численно коэффициенты матрицы Аг, вычисляются следующим образом. Возьмем вектор 6, = (0, 0,0,..., 1,..., 0,0) (единица находится на ,7-том месте). На значениях этого вектора строим кубический сплайн Б, (7) [14]. Интегрируем этот сплайн аналитически по рекуррентной формуле

Тем самым находим 7-й столбец матрицы Аг, . Пусть а, = Ап,, 7 = 1,..., п. Тогда

Кроме того, с помощью фундаментальных сплайнов аналогичным образом построим матрицу для дифференцирования

о

и оператора Гильберта

о

Тогда

где

о

А у = 0, Лі+у = Л*, + Б, (7) ^7, * = 1,...,п - 1.

(25)

П

У=1

Оператор Гильберта дискретизируется формулой

П

где Н, - коэффициенты квадратной матрицы, которые находятся следующим образом.

Представим оператор Гильберта в виде

и учтем, что

А* (/?) = I [А (7) - А (/?)] (]П

о

7 - в4

Иш [А (7) - А (/3)] сіш

V 2

2А' (в).

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

Если г =1, г = п, то

^2 а,с^ у'=1>у'=г

Ъ - 7г

у=1

у'=1>у'=г

7, - 7г

Если * = 1,то

У = 2

—1

у = 1

у = 1

В последней формуле учтено, что А' (0) = А' (2п), а значит,

п п

А' (0) йуА,-, А' (2п) ^А,.

5=1 ,=1

*=п

в силу условия м(0) = м(2п) ■

М (7п) = М (71)

Собираем коэффициенты при А, в каждой строке с номером г, в итоге для Н

получим соотношения:

Яи =-------------

2тг

П—1

2 (аігіи Н- ^п^пі) ^ ^ ^ —-

3 = 2

і н1п — (^І^іп + (іп(іпп),

Яу - “і

, 7 = 2, ...,п - 1,

Ни =--------------

у 2тг

2агйго + а, ctg

7, - 7г

, * = 2, ??. — 1, 7 = 1,..., п, і 7,

Н.и =-----------------

2тг

* = 2, п - 1, 7 = 1,..., ?г, і = 7,

НП, Н1, , 7 1, . . . , п.

В дискретной форме уравнения (20), (23), (24) примут вид:

1 ~

р. ^ ^ (1-і Ті

2п А—*

П

(27)

(28)

Получили систему (26)—(28) из п+2 уравнений относительно п+2 неизвестных: а, Щ, , ] = 1, Эта система решалась численно методом Ньютона, причем квадратурные коэффициенты вычисляются до начала итераций. Вид

системы (26) (28) позволяет на каждом шаге итерационного процесса заполнять якобиан системы аналитически, что является основным ресурсом экономии машинного времени. Неравномерная сетка узлов 7* строилась следующим образом. Исходный профиль задается в виде дискретного множества точек ж*, у*, г = 1,..., п (последняя точка совпадает с первой). С помощью численного конформного отображения внешности профиля на внешность круга единичного радиуса эти точки переводятся в точки 7* параметрической окружности. Тем самым определяется сетка 7* и находится нулевое приближение

где 0* — известные углы наклона касательной в точках ж*, у*. Это нулевое приближение в случае с = 0 (несжимаемой жидкости) будет удовлетворять системе

(26) (28) с очень высокой точностью, определяемой ошибками дискретизации и численного конформного отображения.

На рис. 4 показаны распределения приведенных скоростей А по профилю NACA 2411. рассчитанные при различных числах Маха на бесконечности и различных углах атаки аа . Расчеты были проведены по модели газа Чаплыгина (сплошные линии) с помощью CFD пакета FLUENT (пунктирная линия) и по формуле Кармана Дзяна (штриховая линия). Из приведенных графиков видно, что в докрнтнческом диапазоне чисел Маха все три подхода дают очень близкие распределения. Отличия наблюдаются лишь для критических значений чисел Маха на бесконечности. Совершенно аналогичные результаты были получены для профилей NACA 0012 и

Отметим, что формула Кармана Дзяна также выводится в предположениях газа Чаплыгина, однако при ее выводе не учитываются различия между формами профилей в фиктивном потоке н физической плоскости. Наши численные эксперименты показывают, что для реальных профилей эти различия крайне незначительны, поэтому неудивительно, что формула Кармана Дзяна дает результат, близкий к полной модели газа Чаплыгина.

4. Результаты числовых расчетов

KLARK К.

Рис. 4. Распределения приведенной скорости по профилю NACA 2411

Укажем также, что расчет во FLUENT на уже построенной сетке одного распределения приведенной скорости занимает около 3 мин, в то время как расчет по полной модели газа Чаплыгина и по формуле Кармана Дзяна занимает не более 3 с.

Работа выполнена при финансовой поддержки РФФИ (проект X- 12-01-00996, 12-01-00333).

Summary

Е.М. Kutlyar, D.V. Maklakov. Flow past a Profile of a Given Shape in the Chaplygin Gas.

In the paper, the problem of a subsonic continuous flow past a profile of a given shape is investigated. The problem is reduced to a nonlinear Villa-type integral equation, which after discretization is solved by Newton’s method. It is shown that the velocity distributions

computed using the Chaplygin gas model, the Karman Tsien formula and the CFD Fluent package are very close to each other in the subsonic range of Macli numbers.

Key words: subsonic flows, Chaplygin equations, Chaplygin gas. Villa’s equation, discretization, Newton’s method.

Литература

1. Чаплыгин С.А. О газовых струях // Чаплыгин С.А. Поли. собр. соч. Л.: АН СССР, 1933. Т. 2. С. 3 90.

2. Гуревич М.И. Теория струй идеальной жидкости. М.: Наука, 1979. 536 с.

3. Суигурцев Ю.В. Плоские струйные течения газа. М.: Изд-во Моск. ун-та, 1989. 256 с.

4. Тумалиев Г.Г. Нахождение формы профиля по заданному распределению скорости

с учетом сжимаемости жидкости // Изв. Казап. физ.-матем. о-ва. Сер. 2. 1945.

Т. 13. С. 127 132.

5. Тумалиев Г.Г. Определение формы границ потока жидкости по заданному распределению скорости или давления // Учен. зап. Казап. уп-та. 1952. Т. 112, кп. 3.

С. 3 41.

6. Тумалиев Г.Г., Нужии М.Т. Обратные краевые задачи и их приложения. Казань:

Изд-во Казап. уп-та, 1965. 333 с.

7. Woods L.C. The theory of subsonic plane flow. Cambridge: Cambridge Univ. Press,

1961. 594 p.

8. Седое Л.И. Плоские задачи гидродинамики и аэродинамики. М.: Наука, 1980. 448 с.

9. Daripa Р.К., Sirovich L. Exact and approximate gas dynamics using the tangent gas // J. Comp. Phys. 1986. V. 62, No 2. P. 400 413.

10. Елизаров A.M., Ильинский Н.Б., Пота,шее А.В. Обратные краевые задачи аэро-

гид родипамики: теория и методы проектирования и оптимизации формы крыловых профилей. М.: Физматлит, 1994. 436 с.

11. Степа,нов Г.Ю. Гидродинамика решеток турбомашип. М.: Физматгиз, 1962. 512 с.

12. Елизаров А.М., Касимов А.Р., Маклаков Д.В. Задачи оптимизации формы в аэрогид родипамике. М.: Физматлит, 2008. С. 18 29.

13. Маклаков Д.В. Нелинейные задачи гидродинамики потенциальных течений со свободными границами. М.: Япус-К, 1997. 281 с.

14. de Boor С. A Practical Guide to Splines. Berlin: Heidelberg: N. Y.: Springer-Verlag,

1978. 346 p.

Поступила в редакцию 25.11.11

Котляр Евгения Михайловна аспирант кафедры аэрогидромехапики Казанского (Приволжского) федерального университета.

E-mail: Kem-vgiMimail.ru

Маклаков Дмитрий Владимирович доктор физико-математических паук, профессор кафедры аэрогидромехапики Казанского (Приволжского) федерального университета.

E-mail: Dinitri.MaklakovQksu.ru