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

УДК 532+681.3

ЧИСЛЕННО-АНАЛИТИЧЕСКИЙ МЕТОД РЕШЕНИЯ ЗАДАЧ ПОТЕНЦИАЛЬНОГО ТЕЧЕНИЯ ОКОЛО ГРУППЫ ДВУМЕРНЫХ ТЕЛ

© 2004 В. А. Фролов Самарский государственный аэрокосмический университет

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

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

Модель потенциального обтекания таких объектов является наиболее простой и может быть использована в качестве одной из составляющих в более общей модели течения. Иногда при аэродинамическом проектировании летательных аппаратов (ЛА) модель потенциального течения является достаточной для получения аэродинамических характеристик на стадии предварительного проектирования. Модели на основе уравнений Эйлера или Навье-Стокса при аэродинамическом проектировании, когда выполняется оптимизация, требуют больших затрат машинного времени. Однако даже модель потенциального течения около группы тел может привести к большим временным затратам, если напрямую использовать методы особенностей (панельный метод [1, 2], МДВ [3], комплексный метод граничных элементов [4, 5]). С математической точки зрения эти методы можно отнести к группе численных методов граничных элементов, которые, в конечном счете, сводятся к решению системы линейных алгебраических

уравнений (СЛАУ) с матрицей плотной структуры. Порядок системы определяется общим количеством граничных элементов, на которые разбиваются образующие тел.

Рассмотрим частный случай, когда каждый контур имеет приблизительно одинаковую длину образующей. Разобьем образующую каждого тела на N граничных элементов. Если число контуров принять за к, то общее количество элементов равно к^. Известно, что машинное время при решении СЛАУ с матрицей плотной структуры методом исключения пропорционально (к-^3, а при применении итерационных методов -т(к-Щ2, где т - количество необходимых итераций для достижения заданной точности. Ясно, что при т <k■N итерационные методы оказываются более экономичными по сравнению с методами исключения. Для удов -летворения высокой точности требуется большое число элементов, что независимо от выбранного численного метода решения СЛАУ неизбежно будет приводить к значительным затратам компьютерного времени. Поэтому поиск путей повышения экономичности методов решения потенциальных задач обтекания тел является актуальной задачей вычислительной гидроаэродинамики.

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

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

Суть предлагаемого ниже численноаналитического метода состоит в совместном применении конформного отображения ТФКП и МДВ. Использование инверсии дискретных вихрей относительно окружности, на которую преобразуется центральное тело, позволяет обеспечить точное выполнение условий непротекания поверхности центрального тела. Условие непротекания тел, расположенных около центрального тела, обеспечивается решением СЛАУ Подобная идея использования конформного отображения ТФКП и присоединенных дискретных вихрей использовалась в работе [10] для расчета присоединенных масс эллиптического контура, двигающегося вблизи неплоской бесконечной границы.

Использование симметрии задачи позволяет рассматривать одновременно обтекание двух или четырех одинаковых тел, расположенных около центрального тела симметрично. При этом сохраняется порядок СЛАУ как для одного тела. Симметрия течения используется при вычислении матрицы аэродинамического влияния. Основной трудностью задачи является проблема отыскания преобразующей функции конформного отображения внешности центрального тела в физической плоскости на внешность окружности во вспомогательной плоскости. Идея метода демонстрируется на примерах известных функций преобразования двуугольника [11] и эллипса [12] в физической плоскости на окружность во вспомогательной плоскости. Геометрия тел, расположенных около центрального тела, может быть произвольной. Возможны следующие варианты геометрии задачи:

- обтекание изолированного центрального тела;

- обтекание центрального тела в присутствии одного внешнего тела;

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

- обтекание комбинации тела с четырьмя симметрично расположенными внешними телами.

Рис. 1. Физическая и вспомогательная плоскости комплексного переменного при конформном преобразовании двуугольника в окружность

Модель позволяет решать задачи при отсутствии и наличии циркуляции. Наложение циркуляции может быть осуществлено независимо как на центральное тело, так и на внешние цилиндры.

2. Модель потенциального течения около группы тел. Рассматривается течение идеальной несжимаемой жидкости около группы двумерных тел. Пусть физическая область течения О ограничена кривой Ь, являющейся образующей центрального тела, и одинаковыми симметрично расположенными произвольными контурами £>1, ^, £ и

Д4. Пусть контуры £ и £ симметричны контурам £ и П3 относительно оси ОХ , а контуры П3 и симметричны контурам £

и £ относительно оси ОУ (рис. 1 и рис. 2).

Как известно, потенциальное течение описывается уравнением Лапласа

на бесконечности

Р V, 00В(У,, г),

дх

р V, 00$(Г',,у),

(3)

Ар = 0,

(1)

где р - потенциал скорости.

Для неподвижного контура граничные условия записываются в виде: на границе контура

V, = р 0;

П ^ ’

дп

(2)

где п - внешняя нормаль к образующей кон -

тура; х, у - декартовы координаты; V,, V, -

нормальная составляющая скорости на границе тела и скорость течения на бесконечности, соответственно.

Фундаментальное решение задачи (1)-(3) можно получить методом ТФКП, используя принцип суперпозиции частных фундаментальных решений для диполя; N присоединенных вихрей, распределенных по внешнему контуру П и линейной функции для плоскопараллельного потока. Интенсивности дискретных точечных вихрей должны быть выбраны так, чтобы выполнялись условия непротекания на всех телах одновременно. Непротекание на контуре Ь центрального тела обеспечивается аналитически применением инверсии вихрей относительно окружности во вспомогательной плоскости комплексного переменного. Непротекание на контуре П1 выполняется удовлетворением равенства нулю нормального компонента скорости в контрольных точках этого контура.

Рис. 2. Физическая и вспомогательная плоскости комплексного переменного при конформном преобразовании эллипса в окружность

Разобьем произвольный контур В на N элементов, как показано на рис. 1 и рис. 2. Расположим на каждом элементе точечный дискретный вихрь и контрольную точку. Координаты этих точек расположим на расстоянии 0,25 и 0,75 длины элемента соответственно. Комплексные координаты таких точек в физической плоскости равны

г . = х + іу ., г ■ = х + іу ■

У У •'У}’ С Сі У Сі 5

где і2 = -1 - мнимая единица;

координаты ,-го точечно-

(, уч), (ха, уС1) -

го вихря и координаты /-ой контрольной точки, соответственно.

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

(4)

где д = % + іп - комплексная переменная во вспомогательной плоскости; I - полухорда двуугольника (рис. 1); г0 - радиус окружности во вспомогательной плоскости; параметр определяется формулой

к =

п

2{ж-8У

где 5 - угол при пересечении окружностей двуугольника (рис. 1).

Конформное преобразование внешности эллипса в физической плоскости на внешность окружности выполняется при помощи функции преобразования Н. Е. Жуковского [12]

д = г + л/ г2 -с2

(5)

где с2 = а2 - Ь2 - квадрат межфокусного расстояния эллипса; а, Ь - большая и малая полуоси эллипса.

В зависимости от геометрии контура Ь комплексные координаты во вспомогательной плоскости у-го точечного вихря дч- и

/-ой контрольной точки дс/ определяются по формулам (4) или (5).

Комплексный потенциал течения во вспомогательной плоскости может быть представлен в виде суммы пяти составляющих: комплексного потенциала обтекания цилиндра; комплексного потенциала, определяющего циркуляцию скорости по контуру Ь; комплексного потенциала от суммы точечных вихрей, расположенных на контуре В1 , и их инверсий; комплексного потенциала от суммы точечных вихрей, расположенных на кон -туре В2 с соответствующей инверсией вихрей; комплексного потенциала от суммы точечных вихрей, расположенных на контуре В3 и В4 с соответствующей инверсией вихрей.

Таким образом, суммарный комплекс -ный потенциал течения можно представить в виде следующей суммы:

^ = '№1 + А^2 + А 2 ^3 + А 3 ^4 + ^,

(

V д-

ад 7

Vад Г

2

ад 0

д у

2пі

}=1

ІП

(д-д)

2пі

ІГ

,=1

V

(

ІП

V ду у

(д-ду

(6)

2пі

(д

ІП

’У У

( 2 'А

д + — д

0_ у у

(д + ду,)

г0

2пі

Іп д,

где А1 = А2 = А3 = 0 для случая обтекания

изолированного контура Ь; А1 = А 2 = А 3 = 1 для случая обтекания контура Ь и симметрично расположенных контуров Б и Б2;

1

2

Г

0

2

Г

0

1

W3 =

2

0

1

Г0 - величина циркуляции скорости на контуре Ь; - действительный параметр, ха-

рактеризующий отношение скоростей на бесконечности в физической плоскости г и во вспомогательной плоскости д :

V V

' оо _ оо

V ~V

Voэ,Voэ - комплексная и сопряженная скорости течения на бесконечности во вспомогательной плоскости, соответственно. Величина равна

4/Ь;

для эллипса; для двуугольника.

Комплексно-сопряженная скорость в любой точке области □ и на ее границе находится формулой

— . dw dд

V = и - /V =----------------------,

dд dz

(7)

где и, V - компоненты скорости вдоль оси ОХ и ОУ, соответственно.

Найдем неизвестные интенсивности присоединенных точечных вихрей из решения СЛАУ, составленной из условий непро-текания в контрольных точках внешнего кон -тура:

№Н4 ()

где [А] - матрица аэродинамического влияния для произвольного внешнего контура с учетом его зеркальных отражений относительно осей ОХ и ОУ Элементы матрицы Ау/ вычисляются по формулам (6) и (7) для ком -плексных потенциалов ^2, ^3 и ^4; [г] -вектор-столбец неизвестных интенсивностей присоединенных вихрей Гу , входящих в формулы (6); [^] - правые части, в которые входят нормальные компоненты скорости в кон -трольных точках контура Б Нормальные компоненты рассчитываются по формулам (6)

и (7) для комплексных потенциалов w1 и ^5.

СЛАУ (8) является переопределенной, поскольку должно быть выполнено условие заданной циркуляции для внешнего контура. Случай бесциркуляционного обтекания означает, что значение циркуляции равно нулю.

Рис. 3. Линии тока потенциального течения

около двуугольника (5 = 165°; / = 0,5) и двух цилиндров (Я = 1,5)

Рис. 4. Линии тока потенциального течения

около двуугольника (5 = 1650; / = 0,5 ) и двух эллипсов (Ь = 1,5)

=<

1

у=1

где Ыз,4

нормальный компонент скорости, индуцированный у-м вихрем в /-й кон -

Рис. 5. Линии тока потенциального течения около пяти цилиндров

Чтобы выполнить условие по заданному значению циркуляции, необходимо исключить одну из контрольных точек из решения. Тогда (8) принимает вид

('"Л»,=«„,. '=1. * -1.

у=1

N

= г.>

трольной точке и вычисляемый по комплексным потенциалам ^2, ^3 и ^4; )15 -

нормальный компонент скорости в /-й контрольной точке, рассчитываемый по комплексным потенциалам w1 и ^5; Г - заданная величина циркуляции вокруг контура .

После нахождения распределения интенсивностей присоединенных точечных

вихрей гу компоненты скорости, обусловленные внесением группы тел в однородный поток, определяются по формулам (6) и (7). Компоненты скорости в любой точке облас-

Рис. 6. Линии тока потенциального течения около пяти цилиндров при наличии циркуляции около малых цилиндров (Г1 < 0)

Рис. 7. Линии тока потенциального течения около пяти цилиндров при наличии циркуляции около малых цилиндров (Г1 > 0)

ти Q и на ее границах L, D D D3 и D4 находятся по формулам

u=

Re(V ) v = - Im(V ).

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

dx dy

— = — = dt,

4. Погрешность метода. Граничные условия непротекания на контуре L выполняются точно, поэтому на основании теоремы о максимуме основная погрешность метода сосредоточена на границе D1 внешнего контура. Для оценки данной погрешности вычислим среднеквадратичное отклонение RMS нормального компонента скорости на контуре D1 от нулевого значения по формуле

u v

где dt - произвольная константа, задающая шаг интегрирования.

Интегрирование уравнения (9) проводилось по схеме Эйлера первого порядка

хм = x +Tu {i Ум = У +Tv {i),

где xi+1, yi+1 и xi, yi координаты в (/+1)-ой и i-ой расчетных точках, соответственно; u(i), v(i) - компоненты скорости в i-ой расчетной точке; т « dt = const.

На рис. 3...8 показаны линии тока обтекания различных комбинаций двумерных тел. Результаты, представленные на рис. 5 и на рис. 6, могут быть получены двумя способами: преобразованием (4) двуугольника

и преобразованием (5) эллипса ( ^ b) в окружность.

(9) RMS = ■

—X-n/V7

N -1 “Ю m ■'

(10)

где Vni - нормальный компонент скорости, вычисленный на самом контуре D1 в точке, полученной от пересечения прямой, соединяющей центр тела с i-ой контрольной точкой на i-ом элементе и образующей контура

D1.

В общем случае RMS зависит от взаимного расположения тел, от их геометрии и от заданных значений циркуляций Г0 и Г Максимальная погрешность наблюдается, когда тела касаются друг друга. На рис. 9 показаны типичные зависимости среднеквадратичного отклонения от количества вихрей и величины циркуляциина контуре Dj. В качестве примера выбраны варианты расчета для комбинации тел, изображенных на рис. 4 и рис. 5, когда внешние цилиндры касаются двуугольника. Для варианта комбинации тел, показанной на рис. 5, расчеты величины RMS выполнены для двух значений циркуляции

5 8

rj = 0 и Г = 10. Из рис. 9 следует, что при одинаковом количестве вихрей значение RMS снижается с уменьшением величины циркуляции и количества тел.

5. Заключение. Предложенный численно-аналитический метод расчета потенциальных течений около группы тел обладает высокой точностью и экономичностью по сравнению с прямым применением методов граничного элемента.Экономичность метода и повышенная точность расчета обеспечивается аналитическим выполнением граничных условий на контуре L .

Минимальный выигрыш во времени

счета составляет приблизительно 8 раз ( к3 ,

где к = 2 ) по сравнению с численным методом граничных элементов при одинаковых количествах элементов и дискретных вихрей N на равных по длине образующих двух тел.

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

Список литературы

1. Андерсон Д., Таннехилл Дж., Плет-чер Р. Вычислительная гидромеханика и теплообмен: В 2-х т.: Т. 1: Пер. с англ. - М.: Мир, 1990. - 384 с.

2. Флетчер К. Вычислительные методы в динамике жидкостей: В 2-х т.: Т. 2: Пер. с англ. - М.: Мир, 1991. - 552 с.

RMS

0,04

0,03

0,02

0,01

0

№

I

ш № ч о' II комбинаци я рис. 5

к-. -ч- ГтЧ и 1 ация рис. 5

ы ь- —

=0, комбинация рис.4

0

500

1000

1500

2000

Рис. 9. Зависимость среднеквадратичного отклонения от количества вихрей и величины циркуляции

3. Белоцерковский С. М., Котовский В. Н., Ништ М. И., Федоров Р. М. Математическое моделирование плоскопараллельного отрывного обтекания тел. - М.: Наука, 1988.

- 232 с.

4. Громадка II Т., Лей Ч. Комплексный метод граничных элементов в инженерных задачах: Пер. с англ. - М.: Мир, 1990. - 303 с.

5. Афанасьев К. Е., Стуколов С. В. КМГЭ для решения плоских задач гидродинамики и его реализация на параллельных компьютерах: Учебное пособие. - Кемерово: КемГУ, 2001. - 208 с.

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

7. Фролов В. А. Модель потенциального обтекания комбинации двух круговых кон -туров в присутствии пары дискретных стационарных вихрей //Ракетно-космическая техника. Научно-технический сборник, серия XII, выпуск 1. - Самара: Волжское конструкторское бюро РКК «Энергия», Самарский го -сударственный аэрокосмический университет, 2001. - С. 194-201.

8. Frolov V. A. High-speed flows of the compressible fluid around two circle contours

with a pair of symmetric vortices, The International Summer Scientific School “High Speed Hydrodynamics”, June 16-23, 2002, Cheboksary, Russia, 2002. - P. 331-338.

9. Фролов В. А. Модель потенциального обтекания произвольного контура в присутствии двух круговых контуров // Управле -ние движением и навигация летательных аппаратов. Сборник трудов X Всероссийского научно-технического семинара по управлению движением и навигации летательных аппаратов (Самара, 26-27 июня 2001 г.). - Самара: Самарский государственный аэрокосмический университет, 2002. - С. 315-322.

10. Горбань И. Н. Гидродинамические характеристики эллиптического контура, двигающегося вблизи неплоской границы // Гидродинамика, Республиканский межведомственный сборник научных трудов, вып. 59.

- Киев: Наукова думка, 1989. - С. 64-68.

11. Фабрикант Н. Я. Аэродинамика. Часть первая. - Л., М.: Государственное издательство технико-теоретической литературы, 1949. - 624 с.

12. Лойцянский Л. Г. Механика жидкости и газа. - М.: Наука, 1978. - 736 с.

NUMERICALLY - ANALYTICAL METHOD OF SOLVING OF PROBLEMS OF POTENTIAL FLOW ABOUT GROUP OF TWO-DIMENSIONAL

© 2004 V. А. Frolov Samara State Aerospace University

A new numerically-analytical method of solving problems of potential flow about group of two-dimensional bodies is developed. The method is based on the use of conformal mapping of the theory of complex variable function (TCVF) and the method of discrete vortices (MDV). The efficiency of the method is achieved by using the numerical procedure for one body only. The accuracy of the developed method is higher due to analytical fulfillment of zero normal velocity for the central body contour.