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

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

Том 11, № 4, 2006

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

Д. Н. ГОРЕЛОВ, Д. Г. РЕДРЕЕВ Омский филиал института математики им. С.Л. Соболева

СО РАН, Россия e-mail: gorelov@iitam.omsk.net.ru, d.redreev@mail.ru

A quadrature formula for a singular integral with the Cauchy kernel is constructed. Quadrature considers large curvature of contour in the vicinity of leading edge as well as large gradients of density when tightening contour to the slit. Uniform convergence of formula on a discrete set of contour points is proved.

Введение

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

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

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

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

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

окрестности передней кромки, с учетом сингулярной особенности плотности, которая может иметь место при стягивании контура к разрезу.

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

В комплексной плоскости £ = х+гу рассмотрим сингулярный интеграл Ф(2) с ядром Коши по замкнутому контуру Ь:

Здесь в — дуговая координата, отсчитываемая от некоторой точки 20 Е Ь вдоль контура Ь по часовой стрелке; 7(в) — плотность интеграла. В теории крыла Ф(^) определяет комплексную скорость в точках контура Ь, обтекаемого плоским потенциальным потоком идеальной жидкости. В этом случае функция 7(в) представляет собой касательную составляющую скорости жидкости на Ь или интенсивность вихревого слоя, моделирующего контур Ь.

Предположим, что контур Ь является крыловым профилем с одной угловой точкой 24, вне которой Ь — гладкая кривая с непрерывно меняющейся кривизной. Назовем точку 21 задней кромкой, а точку 20, максимально удаленную от передней кромкой. Расстояние между этими точками Ь = ^ — 201 называют хордой профиля. Без нарушения общности будем полагать координаты х, у безразмерными, отнесенными к хорде профиля, а точки 20 = 0, 21 = 1 (см. рисунок).

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

Но прежде чем дать математическую формулировку поставленной задачи, преобразуем выражение (1.1) к более удобному виду. Представим замкнутый контур Ь в виде объединения двух разомкнутых контуров Ь1, Ь2, имеющих общие концы 20 = 0, 21 = 1. Пусть Ь1, Ь2 задаются однозначными непрерывными функциями

(1.1)

ь

Ьк : у = Д(х), х Е [0,1], к = 1, 2,

и : у = /г(х)

Задание контура крылового профиля.

связанными в точках z0, z1 условиями

fi(0) = /2(0) = 0, /1 (1) = /2(1) = 0,

lim /'(x) = limo /2(x) = -ro. (1.2)

Пусть функции /k(x) при x ^ 1 имеют асимптотику вида -^/x. Такую асимптотику имеет широкий класс аэродинамических профилей. С учетом этого будем полагать, что

/(x) = VxFfc(x), Ffc G C2[0,1]. (1.3)

В соответствии с (1.3) при x ^ 1 запишем

Fk(x) = Fk(0) + Fk(0) x + O(x2),

/fc(x) = Fk(0) Vi + F'(0) x3/2 + O(x5/2). (1.4)

В дальнейшем в (1.1) интегрирование по дуговой координате s будет заменено на интегрирование по переменной x G [0,1]. Якобиан преобразования координат Jk = + (/)2 при x ^ 1 имеет своей асимптотикой выражение

Jk(x) = 2^ ^4x + F|(0) + ... (1.5)

В предельном случае бесконечно тонкого профиля,

когда контуры Lk вырождаются в разрез L, относительная толщина профиля 5 = 0. В этом случае

lim Fk(0) = 0, Jk(x) = 1 + ..., x < 1, k =1, 2. (1.6)

Рассмотрим теперь функцию y(s). Из теории крыла известно, что на гладком контуре при 5 = 0 скорость жидкости конечна и непрерывна вместе со своей первой производной (вне задней кромки). С уменьшением толщины профиля увеличивается кривизна передней кромки, что приводит к безграничному росту скорости при 5 ^ 0. В результате в предельном случае (5 = 0) скорость жидкости на контуре L в окрестности передней кромки имеет асимптотику const /^/x. Функция y(s) определяет скорость жидкости на контуре L и имеет ту же асимптотику. С учетом сделанных замечаний и асимптотических выражений (1.4)—( 1.6) представим функцию 7(s) в виде

Y(s) = 4= Jk7x, gk(x) G C 1(0,1), k =1, 2. (1.7)

Vx Jk (x)

Здесь индекс к соответствует точкам контура Ьк. При х ^ 1

5 = 0 : 7(.) = + .... * = 0 : 7(.) = ^ + ... (1.8)

4х + ^2(0) Vх

Отсюда следует, что функция 7(5) в передней кромке профиля является ограниченной при 5 = 0 и не ограниченной при 5 = 0.

Вернемся к рассмотрению исходного сингулярного интеграла (1.1). При сделанных предположениях его можно представить суммой двух интегралов

ФМ = ±фИ. (;^ Ь. (1.9)

Фк(2) является сингулярным интегралом, если 2 Е Ь. Перейдем в (1.9) от дуговой координаты в к переменной £ Е [0,1]. Тогда

^ = Л (£) С а(в) = £ + г/к (£)

и с учетом (1.7)

1

Фк(2) = 2^2—/ ■ (1-ю)

0

Введем на [0,1] равномерную сетку А = {х, = 3к, к = 1/^, 3 = 0,..., N} и представим

интеграл (1.10) в виде

Ф(2) = £Ф"(2),Ф"и = ^ / / • (111)

, х-

На этом преобразование исходного сингулярного интеграла (1.1) заканчивается. Перейдем непосредственно к постановке задачи о построении квадратурной формулы для сингулярного интеграла (1.1). Пусть контур Ь заменяется на полигон К, составленный из криволинейных элементов:

К,- : у = и, (х) = л/х Ц, (х), х Е [хз-1 ,х,],

(х) Е С1 (х,-1,х,), (х,-1) = ^(х,-1), Ц,(х,) = ^(х,), (1.12)

а заданные функции (х) аппроксимируются функциями

Як,(х) Е С1 (х,-1,х,), (х,-1) = ^(х,-1), ?к,(х,) = (х,). (1.13)

Введем в рассмотрение интегралы, зависящие от комплексной координаты т Е К:

2 N

Ф(т) = $>*(т), Фк(т) = $>,(т), к=1 ,=1

хз

ф, (т) = — [ -__^ (£К К (114)

^(т) 2пг У — £ — ги,(£) ,т Е К (1.14)

Выберем на контурах Ь, К два множества точек:

^0 = {г0к); к =1, 2,3 = 1,..., N}, То = {т0?); к = 1, 2, 3 = 1,•••,N},

2ок) = х0, + г/к(х0,), т0к) = х0, + гик,(х0,), х0, = х,- + к (1.15)

к 2.

Поставим следующую задачу. Пусть задан сингулярный интеграл Ф(2), определяемый формулами (1.1), (1.9)—(1.11), в которых функции /(х),$к(х) удовлетворяют соотношениям (1.3), (1.7). Требуется аппроксимировать Ф(^) интегралами (1.14) путем выбора функций Ц, (х), (х) таким образом, чтобы выполнялось условие

¿ш |Ф(т0к)) — Ф(20к))1 = 0, Е ^0, т0,к) Е Т0. (1.16)

2. Аппроксимация контура и плотности

Исходный контур профиля Ь заменяем полигоном К, составленным из криволинейных элементов Ку, определяемых уравнениями (1.12). Будем задавать функции Цу(х) линейным интерполяционным сплайном вида

Цу(х) = ^(ху-1) + ^(х) - ^(Х-1) (х - ху-1), х € [ху-1,ху]. (2.1)

X у х у_1

Выясним основные особенности аппроксимации контура Ь криволинейными элементами Ку, задаваемыми уравнениями (1.14), (2.1). На стыках элементов иу(ху_ 1) = /(ху_ 1), и^у(ху) = /(ху), но производные от функций (х) терпят разрыв, величина которого имеет порядок О(Л). Для внутренних точек х € [ху_ 1,ху] криволинейных элементов Ку, ^ > 1, имеют место оценки

^(х) - Цу(х) = 0(Л2), ^(х) - цу(х) = О(Л),

Д(х) - иу(х) = ^х 0(Л2), /(х) - ику(х) = х-1/2 0(Л2) + ^х О(Л). (2.2)

Такие же оценки вне передней кромки имеют место, по существу, и при замене контура Ь на полигон, составленный из прямолинейных элементов. Однако в окрестности передней кромки аппроксимация контура прямолинейными элементами дает большие погрешности, особенно для производных. Предлагаемая аппроксимация оказывается существенно лучше. Действительно, из (1.12), (2.1) следует, что уравнения криволинейных элементов, примыкающих к передней кромке, определяются функциями ик1(х) = ^(0)^/х + ^(0) х3/2 + 0(Л5/2). Это выражение совпадает с асимпотическим представлением (1.6) для уравнения исходного контура при х ^ 1. Поэтому для х € [0, Л]

Д(х) - ик 1(х) = 0(Л5/2), Д(х) - ик 1(х) = 0(Л3/2). (2.3)

Аппроксимируем теперь плотность 7(5) интеграла (1.1) на полигоне К функцией ф(а), где а — дуговая координата точки т € К. Следуя (1.9), представим ф(а) на каждом криволинейном элементе Ку в виде

9(а) = ^ у)' ^(х) = + К)2- ^(х) € С 1[ху-1,ху]. (2.4)

С учетом (1.13) будем задавать функции фу(х) интерполяционным сплайном вида (2.1):

Фку(х) = ^(ху-1) + (хх) - хк(ху-1) (х - ху-1), х € [ху-1,ху]. (2.5)

х у х у_1

Применяя теорему о среднем к функции (х), получим оценку

^ (х) - фу (х) = О(Л), х € (ху-1,ху), к =1, 2; ^ = 1,...,Ж (2.6)

Приведенные оценки (2.2), (2.6) показывают равномерное приближение предложенных аппроксимаций контура (1.12), (2.1) и плотности интеграла (2.4), (2.5) во всех точках исходного контура независимо от его толщины. Этот результат обусловлен выделением множителя л/х в уравнениях контура (1.3), (1.12), учитывающего асимптотику контура Ь вблизи передней кромки, и якобиана преобразования координат в выражениях (1.7), (2.4) для плотности. Следует отметить, что выделение якобиана в качестве мультипликативной функции позволяет не только учесть особенность поведения плотности в предельном случае стягивания контура к разрезу, но и свести исходный интеграл (1.1) к сумме интегралов, которые вычисляются точно.

3. Квадратурная формула

Перейдем к непосредственному построению квадратурной формулы. Для этого нужно вычислить интеграл (1.14), полагая т = т^, г = 1,2; р = 1,...,N, и аппроксимируя контур и плотность интеграла соотношениями (1.12), (2.1), (2.5). Представим функции

Фк, (£), и, (£) в виде

х * — £ £ — х_ 1 /—

Фк,(£) = ^к(х,-1) 3_ . + ^к(х,) . _ 3~ , ик,(£) = V £ (ак, + Ьк,£), х, — х, 1 х, — х, 1

х, х, 1 1 О:, = -,-(х,-1)--,-(х,), Ьк, = -[^к(х,) — ^к(х,-1)].

Тогда

х, — х, 1 х, — х, 1 х, — х, 1

Фк, (т^?) = £к (х,-^-!^) + £к (х, ^(т^); (3.1)

р(-)1(т0Г)) = (—1)к+1 ^ [Мк,(т0(?) — ^(т0^))

2 П г I- х, х,_1 х, х,_1

= (—1)к+1 ^ [ — Мк,(т0^)) + (т0;))1; (3.2)

2П ^ х, х,_1 х, х,_1

[ 1 _

Мк,(т0Р)= / — (г

V7? т01) — £ — г^(ак, + Ьк,£)'

' 0р

Х-1 х3

^(т0?) = / - (г) £ £1£( + Ь £). (3.3)

Здесь тОр = х0р + г—х0р (агр + Ьгрх0р). При г = к,р = 3 (точка т0р) Е Кк,) интегралы Мк,, становятся сингулярными и их следует понимать в смысле главного значения по Коши. В общем случае для их вычисления удобно перейти от переменной £ к £ = В результате

Мк, (г'г))

^(#)= / (г) И Д + Ь £2,. (3.4)

т0Р — £2 — г£(ак, + Ьк, £2)

Интегралы (3.4) вычисляются точно, допуская различные варианты в зависимости от значений коэффициентов а^,, Ьк,. В частности, при Ьк, = 0 знаменатель в дробных выражениях под интегралом является кубическим полиномом по £ и его можно представить в виде

т^ — £2 — ¿¿К- + Ьк, £2) = —гЬк, (£ — ¿1;))(£ — ^ — 4?),

где , т = 1, 2, 3, — корни полинома.

(г)) = /* __

У т0г) — ¿2 — «¿(ак, + Ьк, ¿2)

X

3

Выражение (3.1) позволяет представить функцию Ф(т0(г)) в виде

2 N

фО = ££(ху) [ру(к)(т^) + д;к)(т0(;))],

д0к) = = 0. (3.5)

Функция Ф(т0(г)), определяемая формулами (3.2)-(3.5), дает искомую квадратурную рмулу для сингулярного интеграла (1.1). Для оценки точности этой формулы Ф(^0^)), т0р ) представим в виде

ф(4)) = дф(г0;)) + ФгР(г0;)), Ф^р) = ДФ^Р) + Ф^р). (3.6)

Здесь ^0р) € Ьгр, т0(р) € Кгр; сингулярные интегралы Фгр(^0р)), Фгр(т0р)) определяются

формулами (1.11), (1.14), а функции ДФ(г0р), ДФ(т0р)) определяют Ф(г0р), Ф(т0р)) соответственно по контурам ДЬ = Ь - Ьгр и ДК = К - Кгр. Подынтегральные функции в

дф(4?), дф(т0р) ) ограничены, так как точки 4р\ т0р) находятся вне контура интегрирования. В параметрической форме задания оба интеграла сводятся к интегрированию по переменной £ € [0,1]. С учетом (2.2), (2.5) для рассматриваемых подынтегральных функций имеем оценки

|*0? - £| > Л/2, |т0р) - £| > Л/2, $ - £ - /(£) = Т0(Р) - £ - ги(£) + -хР0(Л2),

- -= 4 Т^- [1 + О(Л)], (3.7)

^^ - £ - г/(£) ^т0р) - £ - ш(£) 1 1 '

которые позволяют сделать вывод, что

|ДФ(4?) - ДФ(т0р))|< СЛ 1пЛ, (3.8)

где постоянная С не зависит от Л.

Осталось оценить |Фгр(^0р)) -Фгр(т0р))|. Для вычисления сингулярных интегралов удобно возвратиться к комплексной переменной, полагая

Л , (г)ч 1 Г 1 £г(£) .

Фгр<2»р' = 2пг У ^•' = - г/у•с =£+г/г^

,т, / (г)ч 1 [ ф(т)^т ч 1 Ф р(£) . .

^ ) = 2П^У ТМГ7' ф(т' = - т^игрс£> •т =£ + ги'рК).

Кг 0р

V гр

Представим интегралы ФГр(^0р)), Фгр(т0р)) в виде

Фгр(г0р))= 2тп\ I Л + 2пг / '

0р р 0р ^(4?) г ^ , 1 />к) - ^м

^0р Ч / ^0р Ч

ш , (г)ч Ф(Т0(Р)) Г ^т 1 Г [ф(т) - Ф(т0р))]^т

Фгр(т°Р ) = 2пг I (г) _ т + Т(^ - Т ■ ( )

т т т т

0р ; 2п т(О _ т ^ 2пг У т(Г) т0р - т К т0р

Ьгр

Первые интегралы в (3.9) равны

л/ (г) (г) ГА (г) (г)

"Ц _ г0р — гр-1 — 1п Т°Р — ^Р-1

/ г(г) _ /- г(г) _ г(г) ' / т(г) _ т г(г) _ т(г)

/ г0р Ч гр г0р К ' 0р т гр ' 0р

^гр Кгр

но

4? - г0р) — 0(к2), *(;$) - Ф(т0рг)) — О(к),

г(г) _ г(г) _(г) _ г(г)

^0р_5р=1 — т0р_5р=1 [1 + 0(к)]

(г) (г) (г) (г) [1 + °(к)].

гр ^0р гр ' 0р

Отсюда следует, что разность первых слагаемых в (3.9) является величиной порядка 0(к).

Что касается вторых слагаемых, то их разность имеет порядок 0(ка), 0 < а < 1, если

предположить гёльдеровость функций ^(г), ф(т). Суммируя все оценки, приходим к выводу, что

|Ф(г0р)) - Ф(г0р))| < Ака, гОр) € 20, ^ € Т0, (3.10)

где А — положительная константа, не зависящая от к.

Оценка (3.10) доказывает равномерное приближение построенной квадратурной формулы для сингулярного интеграла (1.1) по замкнутому контуру Ь на множестве точек

(г)

г0р € Ь. Следует подчеркнуть, что это множество точек соответствует равномерной сетке

А — (ж — ^к, к — 1/^, ] — 0,..., N} на интервале [0,1].

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

[1] Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов. М.: Мир, 1987.

[2] Метод граничных интегральных уравнений. Вычислительные аспекты и приложения в механике / Под ред. Т. Круза, Ф. Риццо. М.: Мир, 1979.

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

[4] Габдулхаев Б.Г. Оптимальные аппроксимации решений линейных задач. Казань: Изд-во Казан. гос. ун-та, 1980.

[5] Иванов В.В. Теория приближенных методов и ее применение к численному решению сингулярных уравнений. Киев: Наук. думка, 1968.

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

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

Поступила в редакцию 22 марта 2006 г.