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

УДК 532.5

М. Е. Макарова

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

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

E-mail: masha-mak@mail.ru

Ключевые слова: вихревой слой, профиль, обтекание, комплексный потенциал, циркуляция, метод вихревых элементов.

Постановка задачи. Рассмотрим задачу моделирования плоского стационарного безотрывного обтекания неподвижного профиля, занимающего область С С R2, потоком идеальной несжимаемой среды. Н.Е.Жуковским в его фундаментальной работе "О присоединенных вихрях" (1906) было показано, что данная задача может быть сведена к поиску интенсивности вихревого слоя, расположенного на границе профиля дС, оказывающего воздействие на среду, эквивалентное действию профиля. Вихревой слой, в свою очередь, можно рассматривать как поверхность разрыва касательной компоненты поля скоростей. Величина разрыва равна интенсивности вихревого слоя, поэтому при аналитическом решении задачи ограничимся вычислением предельного значения скорости потока на границе профиля.

Стационарное поле скоростей среды v в области R2 \ С удовлетворяет условиям

'V ■ v(r) = 0,

v(r) ■ n\dC = Vn|dC = 0, (1)

v(r) ^ v.. при \r| ^ TO,

где П — внешняя единичная нормаль к границе профиля дС.

Решение задачи методом конформных отображений. Если обтекаемый профиль является кругом С' единичного радиуса с центром в начале координат, то комплексный потенциал течения в области R2 \ С' имеет вид [1]

F (С) = W. С + + 2G Ln(°, (2)

где £ — точка на комплексной плоскости, Wx — комплексная скорость невозмущенного потока, звездочка означает комплексное сопряжение, G — циркуляция поля скоростей v по любому замкнутому контуру, охватывающему круг C'.

Через комплексный потенциал течения определяется комплексная скорость потока w(£) = (F')*(£) в любой точке течения, в том числе и на границе контура.

Для случая произвольного профиля C задача расчета обтекания будет решена, если известно конформное отображение £ = ((z), удовлетворяющее следующим требованиям: 1) ((z) переводит внешность C на внешность единичного круга; 2) ((ж) = ж; 3) ('(ж) = P = const. При решении практических задач требуется также знание обратного отображения z = Z(£), переводящего внешность единичного круга на внешность профиля C.

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

здесь а, Л, Н — константы, зависящие от вида профиля.

В случае профиля Жуковского а, Л, Н определяются следующим образом:

Ф = р(а, Ь) = аге^(Ь/а), Н = Н(а, й, Ь) = %Ь — йв

Произвольные параметры а, й, Ь задают длину, толщину и искривление соответственно [2].

Для эллиптического профиля

а = a(ai,bi) = J а\ — Ь\, R = R(a\,Ь \) = а,\ + Ь\, H = 0.

Произвольные параметры а 1, Ь 1 задают большую и малую полуоси эллипса соответственно.

Зная функции, осуществляющие необходимые конформные отображения, можно записать выражение для комплексного потенциала течения [1]:

R = R(a, d, h) = у/(а + d cos(^))2 + (h + d sin(^))2,

(3)

где w00 — комплексная скорость набегающего потока при обтекании профиля, которая известна из постановки задачи.

Для вычисления величины , входящей в (2), отметим, что

dF

W * = — W ж d?

df (Z)

dZ

dZ (?)

Z ^ж

d?

ж

=

dZ

d?

R *

= 2 ^;

(4)

поскольку для рассматриваемых отображений эллипсов и профилей

^ Л

Жуковского — ^ — при 4

С учетом (2) и (4) комплексный потенциал течения окончательно

имеет вид

R

G

(5)

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

V = < дфдФ >, (6)

dx dy

где Ф(х,y) = Re(f (x + iy)).

Для определения скорости на профиле введем параметризацию его границы дС и будем считать, что г = г(р), р Е [0, 2п). Наиболее простой способ введения параметра р — это использование известного отображения £ (4), тогда г(р) = С(ег(р-^), где в случае профиля Жуковского вычисляется по формуле аг^(Л,/а), а для эллипса ^ = 0.

Тогда

/ (г(р)) = ^ (¿(г(р))) = ^ (¿(С (е^))) = ^ (е^) =

W

G

= W* ei(p-v) + —^г + — Ln(ei(p-v)).

ei(p-v) 2ni

Положим Wi» = |Woo |eie, следовательно, WJ* = |Woo|e" f (z(p)) = |W00|ei(p-v-e) + |W00|e-i(p-v-e) + - =

2П

-iß

и

= 72|w ж |(e

i(p-v-ß) + e-i(p-v-ß)) + G(p - ^

2n

= R|Voo| cos(p - ^ - в) +

G(p - У)

2n

Для вычисления скорости на границе профиля отметим, что

df = df dp dz

df

откуда w*(z) = — dz

• dz = R\v^\sin(^ + ß - p) + de dp 2n

dC

dz

dp

dz

-1

G

R\v^\sin(^ + ß - p) + — ) •

Вычисляя производную —, получаем

df

dz

dp

G

R\v^\ sin(^ + ß - p) + —

2П

de

iRei(p-v)

2

1 -

(Rei(P-v)) + H)2

(7)

dz

Можно убедиться, что при моделировании обтекания эллипса

— = 0, в то время как для профиля Жуковского на острой кромке

dp

производная обращается в ноль. Согласно постулату Чаплыгина-Жуковского-Кутты, в этом случае для обеспечения конечности скорости на острой кромке числитель дроби также должен обращаться в ноль. Из этого условия определяется циркуляция скорости по профилю Жуковского, которая пропорциональна скорости набегающего потока, зависит от формы профиля и угла атаки в:

G(ß) = -2n\v^\sin(e + <p)(V h2 + а2 + d).

(8)

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

Известно, что скорость безвихревого потока несжимаемой среды в точке (х; у) области течения можно определить по интенсивности вихревого слоя 7(р0) = 7(хо,уо), моделирующего профиль, используя закон Био—Савара:

v(x,y) = v^ +

de

Y(Po) (-(У - yo)i + (x - xo)j) 2n ((x - xo)2 + (y - yo)2)

d1(xo; yo),

(9)

(х,у) е М2 \ С.

Для упрощения выкладок рассмотрим трехмерное поле скоростей V = (ух,Уу, 0)т, и будем считать, что точка вычисления скорости задается радиусом-вектором г = (х,у, 0)т, а точка на границе профиля, по

2

а

координатам которой ведется интегрирование, — радиусом-вектором

Го = (xo,yo, 0)т. Пусть Y(r) = (0, 0,7(г))т, тогда

¿(Г) = V» + / Y^.X (ГГ-2Г0) «и™, Г € R2 \ С.

J 2n|r - Го|2

дС

Предельные значения скорости потока v±(r) на границе профиля дС определяются по формуле [3]

ад = / Y^—^ ± (^ X n(r)) , r € ЭС, (10)

дС

где n(r) — единичная внешняя нормаль к профилю в точке г. Величина v+(r) и знак плюс в правой части равенства (10) соответствуют предельному значению скорости со стороны потока, а V_(r) и знак минус — предельному значению со стороны профиля. В выражениях (9) и (10) вектор V» = (v»x, v»y, 0)т — заданная скорость невозмущенного потока на бесконечности.

Для нахождения неизвестной интенсивности вихревого слоя 7(Г) используется условие равенства нулю предельного значения скорости V_(r) со стороны профиля. Можно доказать [3, 4], что для выполнения этого условия в методе вихревых элементов достаточно потребовать равенства нулю либо нормальной, либо касательной компоненты вектора V_(r).

Первый метод, соответствующий обнулению нормальной компоненты скорости и получивший широкое развитие и обоснование в трудах отечественных исследователей [3, 5], предполагает скалярное умножение уравнения (10) на орт нормали к профилю и приводит к интегральному уравнению вида

/ _ 7Ы X (r- го) = _ V», r € dC, (11)

J 2n|r -ro|2

дС

Такой подход будем сокращенно называть НМВЭ (метод вихревых элементов с нормальными компонентами скоростей).

Ядро уравнения (11) неограничено и имеет неинтегрируемую особенность при |г — Го| —У 0 , интеграл понимается в смысле главного значения по Коши, а само уравнение представляет собой сингулярное интегральное уравнение первого рода. Его решение существует в силу вида правой части, но не единственно; одним из способов выделения единственного решения является задание значения интеграла

(r)d/r = G, (12)

дС

определяющего циркуляцию поля скоростей вдоль любого контура, охватывающего обтекаемый профиль [3].

Если граница профиля дС — гладкая кривая (в частности, эллипс),

то единственное решение можно выделить равенством = 0.

дС

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

то 7 (г)дХг = О (в), где О(в) определяется по формуле (8).

дС

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

Выделение главного значения по Коши в сингулярном интеграле уравнения (11) обеспечивается расчетной схемой на профиле, которая предполагает чередование расчетных точек, в которых располагаются вихревые элементы, и контрольных точек, в которых обеспечивается выполнение дискретного аналога уравнения (11). При этом для получения верных результатов контрольные точки должны располагаться строго посередине между расчетными. Данная система уравнений дополняется дискретным аналогом уравнения (12) и становится в общем случае несовместной. Вопросам ее регуляризации посвящен ряд работ, в частности [3].

Второй метод решения интегрального уравнения (10) состоит в обеспечении равенства нулю касательной компоненты предельного значения скорости у-(г) на профиле. Такой подход описан в работе [4] и предполагает скалярное умножение уравнения (10) на орт касательной к профилю. Направление касательной т(г) выберем так, чтобы п(г) х т(г) = к, где к — орт оси Ог. В результате приходим к интегральному уравнению

/Г(г) • ?^Х-У ** - ^ = -т« • г е дС, (13)

дС

или

£<<(г - Го)!(?о)(йг0 - 12) = -Т(г) • Уж, г е дС, (14)

дС

- \ (r - Го) ■ n(r)

где Q(r — го) =-гц—ц-тг--ядро интегрального уравнения.

2п|г — г0|2

Этот подход будем сокращенно называть КМВЭ (метод вихревых элементов с касательными компонентами скоростей).

Анализ ядра Q уравнения (14) позволяет сделать вывод, что при

гладкости контура dC оно является ограниченным, Q ^ — при

4п

|r — r0| ^ 0, где к — кривизна границы профиля в точке г. Уравнение (14) является интегральным уравнением Фредгольма 2-го рода. Для его решения также будем разбивать контур на панели, однако расчетная схема может быть достаточно произвольной: в силу ограниченности ядра нет жесткого требования к расположению расчетных и контрольных точек. Будем полагать интенсивность вихревого слоя постоянной на каждой панели. При этом интеграл в (14) заменяется суммой интегралов по отдельным панелям, каждый из которых вычисляется аналитически. Равенство нулю касательной компоненты скорости будем обеспечивать не в отдельных контрольных точках, а в среднем на каждой панели; соответствующие интегралы вычисляются численно с использованием квадратурных формул Гаусса.

Поскольку решение уравнения также не является единственным, вводится дополнительное условие вида (12), а получаемая система регуляризируется так же, как и в НМВЭ.

Результаты расчетов. Рассмотрим несколько тестовых примеров, для которых методами ТФКП удается построить отображение и вычислить аналитически интенсивность вихревого слоя на контуре, пользуясь формулой (7). Далее приведены вычисленные НМВЭ и КМВЭ значения интенсивностей вихревых слоев на профиле в сравнении с точными решениями для различного числа n панелей, а в

n

таблицах — погрешности численных решений ЦДЦ^ = lYi — 70|1¿,

i=1

|| Д |оо = max ^ъ — Y0|, где Yi — найденное в расчетах среднее значение

ii

интенсивности вихревого слоя на i-й панели, yÍ0 — среднее значение интенсивности на панели, 1ъ — длина панели.

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

Таблица 1

Расчет обтекания кругового профиля

Число 50 200 500

панелей,n НМВЭ КМВЭ НМВЭ КМВЭ НМВЭ КМВЭ

CondA 14 14 32 32 80 80

IIAIk 0, 0045 0,0074 0,0003 0,0008 0, 0001 0, 0002

IIAIU 0, 0011 0,0018 0,0001 0,0002 0, 0000 0, 0001

а 6

Рис.1. Интенсивность вихревого слоя при обтекании эллипса (ах/Ъх = 10), рассчитанная методами НМВЭ (а) и КМВЭ (б) в сравнении с точным решением

На рис. 1 показаны результаты расчета интенсивности вихревого слоя при расчете обтекания эллиптического профиля с соотношением полуосей а\/Ь\ = 10, установленного под углом атаки в = п/6, по сравнению с точным решением при дискретизации профиля на п = 200 панелей.

Значения числа обусловленности соответствующих матриц, а также погрешности численных решений для этого случая приведены в табл. 2.

Таблица 2

Расчет обтекания эллиптического профиля (а^/Ъ^ = 10)

Число 50 200 500

панелей, n НМВЭ КМВЭ НМВЭ КМВЭ НМВЭ КМВЭ

CondA 43 115 90 242 144 384

IIAIk 1, 4505 0, 621 0,4193 0,0058 0,1728 0, 0015

IIAIU 1, 7864 0, 3894 0, 6455 0,0444 0, 2717 0, 0083

Число обусловленности матрицы системы в КМВЭ несколько больше, чем в НМВЭ, однако КМВЭ в данном случае позволяет получить

Рис.2. Интенсивность вихревого слоя при обтекании эллипса (0,1/61 = 20), рассчитанная методами НМВЭ (а) и КМВЭ (б) в сравнении с точным решением

на 1-2 порядка более точное решение по сравнению с НМВЭ. Отметим, что в расчетах НМВЭ основная погрешность сосредоточена у концов большой оси профиля (см. рис. 1).

При рассмотрении более тонкого профиля (эллипса с полуосями 01/61 = 20 под углом атаки в = 0) погрешность НМВЭ даже при использовании п = 200 панелей оказывается достаточно большой, в то время как расчет КМВЭ обеспечивает намного большую точность (рис. 2)

При увеличении числа панелей погрешность обоих методов уменьшается, при этом точность КМВЭ существенно выше. Числа обусловленности матриц КМВЭ, как и в предыдущем случае, несколько выше, чем в НМВЭ (табл. 3).

Таблица 3

Расчет обтекания эллиптического профиля (ai/bi = 20)

Число 50 200 500

панелей,n НМВЭ КМВЭ НМВЭ КМВЭ НМВЭ КМВЭ

CondA 77 264 167 632 375 1011

IIAIk 1, 3164 0,0009 0,4439 0,0001 0,1876 0, 0000

IIAIL 0, 4244 0,0082 0,1461 0,0061 0, 0613 0, 0012

При расчете обтекания симметричного профиля Жуковского с параметрами а = 1,0, Н = 0, ( = 0,2, установленного под углом атаки в = п/6, наблюдается следующий эффект: при использовании НМВЭ численное решение существенно расходится с аналитическим вблизи острой кромки и на носке профиля, чего не наблюдается в КМВЭ (рис. 3).

При расчетах методом НМВЭ с увеличением числа панелей ошибка численного решения в норме || • ||оо растет, а в норме || • Ц^ уменьша-

а 6

Рис. 3. Интенсивность вихревого слоя при обтекании симметричного профиля Жуковского, рассчитанная методами НМВЭ (а) и КМВЭ (б) в сравнении с точным решением

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

Таблица 4

Расчет обтекания симметричного профиля Жуковского

Число 50 200 500

панелей, n НМВЭ КМВЭ НМВЭ КМВЭ НМВЭ КМВЭ

CondA 2.3 ■ 103 3.0 ■ 102 7.0 ■ 104 2.0 ■ 103 6.9 ■ 105 7.0 ■ 103

ЦЛЦг1 0, 560 0,018 0, 309 0, 002 0,161 0, 000

IIAIL 4,179 0,065 21, 218 0, 034 57,138 0, 024

Отмеченные для симметричного профиля Жуковского эффекты проявляются в еще большей мере для несимметричного профиля Жуковского (а = 3,5, Н = 0,3, ( = 0,4), установленного под тем же углом атаки в = п/6. При использовании НМВЭ ошибка численного решения при приближении к острой кромке растет катастрофически и становится неприемлемо большой даже при небольших значениях п (рис. 4). Резкий рост числа обусловленности матрицы осложняет решение системы линейных уравнений (табл. 5), что также может повлечь значительную погрешность.

Норма погрешности || ■ для НМВЭ резко растет, а || ■ Ц^ снижается крайне медленно при увеличении п; при расчетах методом КМ-ВЭ ошибка получается весьма малой и уменьшается в обеих нормах

а б

Рис. 4. Интенсивность вихревого слоя при обтекании несимметричного профиля Жуковского, рассчитанная методами НМВЭ (а) и КМВЭ (б) в сравнении с точным решением

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

Таблица 5

Расчет обтекания несимметричного профиля Жуковского

Число 50 200 500

панелей, n НМВЭ КМВЭ НМВЭ КМВЭ НМВЭ КМВЭ

CondA 1, 7 ■ 104 4, 5 ■ 102 5, 2 ■ 105 3,1 ■ 103 5,1 ■ 106 1,1 ■ 104

ЦЛЦг1 14, 903 0,100 11, 835 0,010 8, 870 0,003

IIAIL 118, 028 0,138 1638, 860 0,072 7933, 360 0,051

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

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Свешников А. Г., Тихонов А. Н. Теория функций комплексной переменной. - М.: ФИЗМАТЛИТ, 2004. - 336 с.

2. Лаврентьев М. А., Ш а б а т Б. В. Методы теории функций комплексного переменного. - М.: Наука, 1965. - 716 с.

3. Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент. - М.: Янус, 1995. - 520 с.

4. Kempka S. N., Glass M. W., Peery J. S., Strickland J. H. Accuracy Considerations for Implementing Velocity Boundary Conditions in Vorticity Formulations // SANDIA REPORT SAND96-0583 UC-700, March, 1996.

5. Андронов П. Р., Гувернюк С. В., Дынникова Г. Я. Вихревые методы расчета нестационарных гидродинамических нагрузок. - М.: Изд-во МГУ им. М.В. Ломоносова, 2006. - 184 с.

Статья поступила в редакцию 25.10.2011