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

Вычислительные технологии

Том 14, № 1, 2009

Метод криволинейных панелей для решения сингулярных интегральных уравнений с ядром Коши по замкнутому контуру*

Д. Н. ГОРЕЛОВ, Д. Г. РЕДРЕЕВ Омский филиал Института математики им. С.Л. Соболева СО РАН e-mail: gorelov@ofim.oscsbras.ru

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

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

Введение

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

*Работа выполнялась в рамках проекта СО РАН № 117.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2009.

некоторая "мертвая зона" по геометрическим характеристикам замкнутых контуров, в которой известные методы решения СИУ практически не работают.

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

1. Постановка задачи

Рассмотрим задачу обтекания крылового профиля (замкнутого контура L с одной угловой точкой) потенциальным потоком идеальной несжимаемой жидкости, имеющим в бесконечном удалении перед профилем скорость v— (рис. 1). Соответствующую краевую задачу для комплексной скорости v(z) в плоскости z = x + iy можно сформулировать следующим образом. Найти вне контура L аналитическую функцию v(z), удовлетворяющую граничным условиям непротекания жидкости через контур L и затухания возмущенной скорости в бесконечно удаленной точке:

Im {w(z)ei0(z)} = 0, z G L, lim v(z) = = |v-|e-ia. (1.1)

Здесь 0 — угол наклона касательной к контуру L в точке z G L. К граничным условиям (1.1) следует добавить постулат Кутты—Жуковского—Чаплыгина о конечности скорости жидкости в задней кромке профиля, позволяющий определять циркуляцию скорости вокруг контура L.

Контур L будем моделировать вихревым слоем с интенсивностью y(s), равной разрыву предельных значений касательной составляющей скорости жидкости при подходе к контуру L с внутренней и внешней сторон контура (s — дуговая координата). Комплексную скорость вне контура L представим интегралом Коши:

e(z> = «г- + i / z-dS) • z (s> GL. a.2)

L

Рис. 1. Задача обтекания крылового профиля

Требуя выполнения граничных условий (1.1) для г>(г), придем к двум интегральным уравнениям относительно искомой вещественной функции 7(в) [2-4]:

2Y(s) + Re <( eiö(z)

v^ +

1

Y(a) da

2пг J z — Z(a)

L

Im < e"

(z)

v™ +

Y (a)da

2тл / z — Z (a)

L

0, z,Z G L;

0, z,Z G L.

(1.3)

(1.4)

Уравнения (1.3), (1.4) являются соответственно сингулярными интегральными уравнениями второго и первого рода. Отметим две основные особенности этих уравнений. Первая из них связана с тем, что на гладких участках контура Ь ядро интегрального уравнения второго рода является ограниченной функцией, тогда как ядро уравнения первого рода имеет полюс первого порядка. Вторая особенность связана с толщиной профиля и проявляется при сближении верхней и нижней его сторон, когда начинают действовать формулы Сохоцкого—Племели для предельных значений интеграла Коши. Эту параметрическую особенность можно исключить, если перейти от уравнений (1.3), (1.4) к комбинированной системе уравнений, предложенной Д.Н. Гореловым [3]. С учетом этих замечаний можно ожидать, что решение краевых задач обтекания телесных профилей целесообразно проводить с помощью СИУ второго рода либо комбинированной системы СИУ.

Перейдем к непосредственной постановке задачи. Нужно построить численный метод решения СИУ (1.3), (1.4) на основе метода панелей, удовлетворяющий следующим условиям:

1) уравнения панелей и контура Ь имеют одинаковую асимптотику в окрестности передней кромки;

2) искомое решение на каждой панели представляется в виде, учитывающем асимптотику решения в предельном случае бесконечно тонкого профиля;

3) уравнения криволинейных панелей и представление искомого решения на них позволяют точно вычислить соответствующие интегралы в методе панелей;

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

2. Построение метода

Исходный замкнутый контур L представим в виде объединения двух разомкнутых гладких контуров Li, L2 (рис. 2), которые задаются функциями y = fk(x), x G [0,1], k =1, 2, связанными следующими условиями:

fi(0) = f2(0)=0, fi(1) = f2 (1) = 1,

lim fi(x) = lim f2(x) = (2.1)

x—0 x—

Предположим, что fk(x) при x ^ 0 имеют асимптотику вида \J~x. Такую асимптотику имеет широкий класс аэродинамических профилей, включая профили Жуковского. С учетом этого представим f1(x), f2 (x) в виде

fk(x) = V^Fk(x), Fk G C2[0,1]. (2.2)

Li : у = fi(x)

О X— x

М : у = f2(x)

Рис. 2. Вид исходного замкнутого контура L

Отметим, что на практике контур L может быть задан таблицей координат. При таком задании необходимо получить аналитическое представление контура с учетом требуемой асимптотики (2.2). Эта задача решена в [5].

Первый шаг в методе криволинейных панелей — аппроксимация исходного контура L кусочно-непрерывным полигоном, составленным из криволинейных панелей. Замена контура L таким полигоном приводит к тому, что интегральные уравнения (1.3), (1.4) начинают зависеть от погрешности аппроксимации как контура, так и углов наклона касательных. Анализ показывает, что представление контура L линейным и квадратичным сплайнами обеспечивает приближение полигона к контуру L, но не способно выработать условие (2.2) о бесконечном значении производной в точке x = 0 при любом конечном числе панелей. Именно это обстоятельство считается основной причиной большой погрешности решения СИУ (1.3), (1.4) стандартным методом панелей вблизи передней кромки, особенно для тонких профилей. Для построения полигона введем равномерную сетку А = {xj = jh, h = 1/N, j = 0,..., N}, на которой определены значения fkj = fk (xj). Исходный контур профиля L заменим полигоном K, составленным из криволинейных элементов Kkj, определяемых следующими уравнениями:

Kkj : y = Ukj(x) = Vx Ukj(x), x e [xj-i,Xj],

Ukj(x) = Fk(xj-i) + Fk(xj) - Fk(xj-i) (x - xj-i), x e [xj-i,xj]. (2.3)

xj xj-i

Уравнения (2.3) при x ^ 0 имеют ту же асимптотику, что и уравнения (2.2). На стыках элементов Ukj (xj-i) = fk (xj-i),Ukj (xj) = fk (xj), но производные от функций Ukj(x) терпят разрыв. Для внутренних точек x e [xj-i,xj] криволинейных элементов Kkj, j > 1, имеют место оценки:

Fk(x) - Ukj(x) = O(h2), Fk(x) - Ukj(x) = O(h),

fk(x) - Ukj(x) = -xO(h2), fk(x) - u'kj(x) = x-i/2 O(h2) + VxO(h). (2.4)

Задание криволинейных панелей в виде (2.3) и оценки (2.4) дают равномерное приближение не только полигона K к контуру L, но и соответствующих производных.

Следующий шаг в предлагаемом методе — выбор специальной аппроксимации для искомой плотности y(s). В предельном случае бесконечно тонкого профиля скорость жидкости на контуре L в окрестности передней кромки имеет асимптотику const/^/x. С учетом этого представим функцию y(s) в виде [6]

Y(s) = -1xjM, gk(x) e Ci(0,1), Jk ^/l + (fk)2, k = 1, 2. (2.5)

Vx Jk (x)

у

Функция 7(5) определена на контуре Ь. На полигоне К заменим ее функцией ф(а), где а — дуговая координата точки т £ К. Следуя (2.5), представим ф(а) на каждом криволинейном элементе Ку в виде

фу (х) = ^ (х-1) +

&(х7) - & О^—О

ху ху_1

(х х у—1), х £ [х у—1, ху ].

Применяя теорему о среднем к функции (х), получим оценку (х) - фу (х) = 0(Л,), х £ (х^_1,х^), к = 1, 2; = 1,

N.

(2.6)

(2.7)

Представление функции ф(а) на криволинейных панелях в виде (2.6) позволяет учесть асимптотику решения в предельном случае бесконечно тонкого профиля, а оценка (2.7) обеспечивает равномерное приближение ф(а) к 7(5) во всех точках исходного контура независимо от его толщины.

Замена в уравнениях (1.3), (1.4) контура Ь полигоном К, а искомой функции 7(5) — на функцию ф(а) приводит к приближенным уравнениям

2д(а) + Б* ^ е^(т)

г><» +

1 Г д(*) ^ т - С(*)

к

0, т,С £ К;

(2.8)

1т ^ е^(т)

^ +

1 Г д(*) ^

т - С (¿)

к

0, т, ( £ К.

(2.9)

Здесь $(т) — угол между касательной к контуру К и осью х в точке с комплексной координатой т(а) £ К.

С учетом (2.6) перейдем в уравнениях (2.8), (2.9) от контура К к криволинейным элементам Ку и заменим ^(¿) на этих элементах функцией фку (х). Тогда

N

т — С (5

к

ТТ) = Е (т), Я«(т) = Ру (т) + Оку_1(т). (2.10)

^ ( ) к=1, 2 7=0

Здесь ду = (ху), а функции Ру, ф^у определяются формулами, интегралы в которых вычисляются точно:

хз+1

Ру (т)

ху+1 х

^х

^^/х т — х — шк7+1(х)

РкМ = 0,

х3+1

(т) =

^х

hv/x т — х — шк7+1(х) к = 1, 2, ^ = 0,...,Ж - 1.

, ^ —1 = 0,

(2.11)

X

3

X

3

Выберем на контуре К множество точек

То = {т, к = 1, 2,] = 1,..., N}, То(к) = хо, + гик,(хо,),

хо, = + V,Л,, 0 < V, < 1. (2.12)

Интегральные уравнения (2.8), (2.9) будем решать методом коллокации, требуя выполнения этих уравнений в точках т, € К. Выбирая в качестве неизвестных величин значения функции Як(ж) в точках х,, к = 1, 2; ] = 0,..., N, приходим к системе линейных алгебраических уравнений, которые соответственно для интегральных уравнений (2.8), (2.9) принимают вид:

1 N

1 ^ ^ ^ С м9(т(г))

! Вгт,к, Як, = 1т | е^(г0т

2 / у / у гт, к] у к] "" ] ^

к=1,2 ,=о

г = 1, 2, т = 0,..., N - 1, (2.13)

1 N (г)

Агт 9тт + Агт+1 9гт+1 + ^ЦТ ^ у ^ ^ Сгт, к, Як, ^е' ( 0т^с

к=1,2 ,=о

г = 1, 2, т = 0,..., N - 1. (2.14)

Коэффициенты Агт, Вгтк,, Сгт)к, определяются по формулам

л __1 Vm л Vm

Агт Т , . Г" ,

2\/Хот Jгm(x0:m) 2\/X0mJгm(x0m)

,к, = -б« {е^-^^т)}, Сгт, к,=е^0т)Дк,(т0(тт^, (2.15)

ег-&(г(г2) = (—1)г+1 1 + гигт (хот) = 11 + гигт (хо

:|

( 1)^+1 ^^от + гигт(жот)/2 + гхоти гт(хот) (2 1д)

1Vжот + '^гт(жот)/2 + гхоти Гт(хот) 1

Системы (2.13), (2.14) содержат 2N — 2 уравнений для 2N неизвестных, поэтому к ним следует добавить еще по два уравнения. При решении задачи стационарного обтекания контура этими уравнениями стали условие непрерывности искомого решения д в точке х = 0 и постулат Жуковского ^1(1) + д2(1) = 0.

Переход от исходных СИУ (1.3), (1.4) на замкнутом контуре Ь к приближенным СИУ (2.8)—(2.11) на полигоне К, составленном из криволинейных панелей Кгт, и дальнейший переход к системе линейных алгебраических уравнений (2.12)-(2.16) обосновываются соответствующей квадратурной формулой, полученной в работе [6].

Вг

3. Численный эксперимент

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

Таблица 1. Относительная погрешность расчета скорости для 20%-го профиля Жуковского

МЛП

X 0 0.025 0.05 0.5 0.95 0.975 1

СИУ-1 0.391 0.022 0.003 0.006 0.284 0.578 125.7

СИУ-2 0.337 0.635 0.266 0.004 0.017 0.019 0.008

СИУ-К 0.500 0.120 0.024 0.002 0.004 0.013 0.004

МК П

X 0 0.025 0.05 0.5 0.95 0.975 1

СИУ-1 0.009 0.022 0.043 0.034 0.912 1.89 6.56

СИУ-2 0.125 0.003 0.007 1Е-4 0.014 0.017 0.002

СИУ-К 0.088 0.009 0.030 0.001 0.003 0.014 8Е-4

Таблица 2. Относительная погрешность расчета скорости для 5%-го профиля Жуковского

МЛП

X 0 0.025 0.05 0.5 0.95 0.975 1

СИУ-1 0.002 0.015 0.005 0.010 0.325 0.605 36.61

СИУ-2 0.397 0.183 0.070 0.002 0.013 0.015 0.022

СИУ-К 0.159 0.009 0.025 0.001 0.003 0.013 0.021

МК П

X 0 0.025 0.05 0.5 0.95 0.975 1

СИУ-1 0.006 0.011 0.031 0.009 0.164 0.330 1.076

СИУ-2 0.004 0.006 0.037 0.006 0.020 0.022 0.001

СИУ-К 0.054 0.011 0.030 0.001 0.007 0.007 4Е-4

Таблица 3. Зависимость максимальной погрешности расчета скорости для 20%-го профиля Жуковского от числа панелей

МЛП

N 20 40 60 80 100 120 140

СИУ-2 1.285 0.635 0.488 0.424 0.389 0.365 0.348

СИУ-К 0.623 0.500 0.424 0.372 0.333 0.304 0.282

МК П

N 20 40 60 80 100 120 140

СИУ-2 0.261 0.125 0.076 0.051 0.035 0.026 0.019

СИУ-К 0.132 0.088 0.068 0.056 0.047 0.040 0.035

Таблица 4. Зависимость максимальной погрешности расчета скорости для 5%-го профиля Жуковского от числа панелей

МЛП

N 20 40 60 80 100 120 140

СИУ-2 0.472 0.397 0.358 0.333 0.313 0.298 0.285

СИУ-К 0.208 0.159 0.145 0.137 0.132 0.127 0.123

МК П

N 20 40 60 80 100 120 140

СИУ-2 0.038 0.037 0.035 0.033 0.032 0.030 0.029

СИУ-К 0.098 0.054 0.037 0.029 0.028 0.027 0.026

уравнений СИУ-К, предложенной Д.Н. Гореловым [3]. Представляет интерес сравнительная оценка точности решения СИУ-1, СИУ-2 и СИУ-К методами криволинейных панелей (МКП) и линейных панелей (МЛП). Следует отметить, что такая оценка уже была ранее для метода линейных панелей [4]. Численный эксперимент проводился для несимметричных профилей Жуковского с относительной толщиной 5 и 20 %, обтекаемых стационарным потоком под геометрическим углом атаки а = 10 Панели на каждой стороне контура выбирались равной длины, а их число в основном расчете N = 40. Результаты расчета представлены на рис. 3 и 4.

В табл. 1 и 2 приведены данные по расчету относительной погрешности е(х) = |^(х)-скорости жидкости в точках верхней стороны контура. Для оценки влияния

МЛП

СИУ-1

МКП

и

2

1

о

О \

1 1 1 о 0.2 0.4 0.6 1 1 д: 1

0 0.2 0.4 0.6

СИУ-2

М

0 0.2 0.4 0.6

П I г

о 0.2 0.4 0.6

И-1-г

о 0.2 0.4 0.6

СИУ-К IV-!

п-1-г

о 0.2 0.4 0.6

Рис. 3. Расчет 20%-го профиля Жуковского методом линейных панелей (слева) и методом криволинейных панелей (справа)

МЛП СИУ-1 мкп

Рис. 4. Расчет 5%-го профиля Жуковского методом линейных панелей (слева) и методом криволинейных панелей (справа)

числа панелей на точность расчета в табл. 3 и 4 дается максимальная относительная погрешность решения СИУ-2 и СИУ-К (в окрестности передней кромки профиля) от числа N. Следует отметить, что все расчеты проводились без какого-либо сгущения расчетных точек (длин панелей) вблизи кромок профилей, как это делается обычно при расчете другими методами.

Эти результаты подтвердили выводы работы [4] о целесообразности сведения краевых задач обтекания замкнутых контуров к решению СИУ-2 или СИУ-К. Главный же результат численного эксперимента — высокая точность решения СИУ-2 и СИУ-К методом криволинейных панелей во всех точках контура профиля, включая его переднюю и заднюю кромки.

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

[1] Белоцерковский С.М., Лифанов И.К. Численные методы в сингулярных интегральных уравнениях. М.: Наука, 1985.

[2] Горелов Д.Н. Методы решения плоских краевых задач теории крыла. Новосибирск: Изд-во СО РАН, 2000.

[3] Горелов Д.Н. Об интегральных уравнениях задачи обтекания профиля // Изв. РАН, МЖГ. 1992. № 2. С. 173-177.

[4] Горелов Д.Н., Смолин Ю.С. Применение системы интегральных уравнений к решению плоских задач теории крыла // Вычисл. технологии. 1999. Т. 4, № 5. С. 24-29.

[5] Горелов Д.Н., Редреев Д.Г. Применение кубических сплайнов для аналитического представления замкнутого контура, заданного таблицей координат // Сиб. журн. индустр. математики. 2005. Т. 8, № 2. С. 26-31.

[6] Горелов Д.Н., Редреев Д.Г. Построение квадратурной формулы для сингулярного интеграла с ядром Коши по контуру крылового профиля // Вычисл. технологии. 2006. Т. 11, № 4. С. 29-36.

Поступила в редакцию 19 февраля 2008 г.