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

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2015

Математика и механика

№ 5(37)

УДК 519.632.4

DOI 10.17223/19988621/37/8

М.А. Пономарева, Е.А. Собко, В.А. Якутенок

РЕШЕНИЕ ОСЕСИММЕТРИЧНЫХ ЗАДАЧ ТЕОРИИ ПОТЕНЦИАЛА НЕПРЯМЫМ МЕТОДОМ ГРАНИЧНЫХ ЭЛЕМЕНТОВ1

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

Ключевые слова: теория потенциала, уравнение Лапласа, осесимметричные задачи, непрямой метод граничных элементов, сингулярные интегралы.

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

Альтернативой непрямому методу граничных элементов выступает метод фундаментальных решений. При его использовании сингулярные интегралы не возникают, однако сложность заключается в выборе точек приложения фиктивных нагрузок. Этот метод для решения осесимметричных задач рассматривается в [5, 6].

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

в области Q, ограниченной осесимметричной поверхностью S с образующей Г

1 Исследование выполнено при финансовой поддержке гранта Президента РФ (МК-3687.2014.1) и РФФИ в рамках научного проекта № 14-08-31579 мол_а.

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

Требуется найти решение уравнения Лапласа

Дм (х) = 0 , х eQ ,

Решение осесимметричных задач теории потенциала

85

(рис. 1) при следующих краевых условиях:

u (х ) = g (х), х e S1; (2)

q(х) =----(—^ = h (х), х e S2, (3)

dn

где S = S1 u S2, q - поток, n - внешняя нормаль к S, g (х) и h (х) - известные распределения потенциала и потока на границах S1 и S2 соответственно. Причем функции g (х) и h (х) не зависят от угловой координаты 0.

Рис. 1. Область решения и система координат (r,0,z)

Гранично-интегральная формулировка задачи

Потенциал, создаваемый в точке х источником единичной интенсивности, сосредоточенным в точке 4, может быть найден из уравнения

АО (х, 4) = -5 (х - 4),

где 5 (х - 4) - дельта-функция Дирака. Решением этого уравнения является ньютоновский потенциал

О (х, 4)= 4-, (4)

4пр

где р - расстояние между точками х и 4. Тогда, если считать, что потенциальное поле создается фиктивными источниками, непрерывно расположенными на границе S, то суммарный потенциал, создаваемый в точке х всеми фиктивными источниками, можно получить интегрированием выражения (4), предварительно умножив его на плотность распределения фиктивных источников ф (4):

u(х) = |ф(4)О(х,4)dS^ ,

S

(5)

86

М.А. Пономарева, Е.А. Собко, В.А. Якутенок

где обозначение dS^ означает, что при интегрировании именно точка 4 «пробегает» поверхность S. В дальнейшем индекс t в подобных случаях будет опущен, однако подразумеваться будет то же самое.

С учетом того, что dS^ = rtd0dГ , выражение (5) преобразуется к виду

и (х) = 11 Ф(4)G(х,4)'d0dr .

Г 0

где Г - образующая поверхности S, расположенная для простоты в плоскости r0z, при этом 0 x = 0 . Так как задача осесимметрична, то ' и ф (4) не зависят от угла

0 и могут быть вынесены за знак соответствующего интеграла:

f 2п 2

и (х) = |1 | G(х, 4)d0 'ф(4)dr .

Г V 0

Введя обозначение

2п

G (х, 4)= | G (х, 4)d0 .

выражение (5) окончательно можно записать в форме

и (х) = | G(х, 4)г^ф(4)dr.

Г

Как видно из рис. 1:

'-fax - )2

+ rX + rt - 2rn cos (

P = yi\Zx~ZS, 'x^'\~ 2'x't

где rx , zx, ', z^ - радиальные и аксиальные координаты точек х и 4 соответственно. Тогда

d0

G,-4)==

222

+ rx + ' - 2rxrt cos 0

()

Пусть

' = (Zx - Zt )2 + rx2 + rt 2rxrt

2b

b = 2rxrt,

(7)

m =----.

a + b

Тогда формула (6) принимает вид

2п d0 = K (m)

G(х,4) = 4- |т 4n Wс

- b cos 0 %4a + b

()

n/2

где

K (m) = I

d0

полный эллиптический интеграл первого рода.

0 л/1 - m sin2 0

Для рассматриваемого осесимметричного случая поток в точке х еГ в направлении внешней нормали n (х), обусловленный наличием единичного сосре-

Решение осесимметричных задач теории потенциала

87

доточенного источника в точке 4, определяется выражением

F (x, 4)

dG

dn

dG dG

---nr +----n_

dr dz

где nr и nz - компоненты внешней нормали n (x). Совершенно аналогично предыдущему:

q (x ) = J F (x, 4 )г=ф (4 )d Г,

Г

где

F (x, 4)= J F (x, 4)d0 .

Нетрудно получить, что

F (x,4)=

1

W a+b

1

2r

2 2 r2 - к +

(- z? )2

a - b

E (m)- K (m)

nr +—x---E(m)nz , (9)

a - b

n/2 _________

где E(m)= I v1 -msin2 0d0 - полный эллиптический интеграл второго рода.

о

Таким образом, гранично-интегральная формулировка задачи с учетом условий (2) и (3) представляется в виде системы уравнений

g(x) = JG(x,4)г?ф(4)dГ , x e S1; (10)

Г

h(x) = JF(x,4)^ф(4)dr, x e S2, (11)

Г

где G (x, 4) и F (x, 4) определяются выражениями (8) и (9), ф (4) - искомая плотность распределения фиктивных источников. После нахождения функции ф (4) становится возможным вычисление потенциала и потока в произвольном направлении в любой точке области Q.

Дискретизация границы и метод решения

Для решения интегральных уравнений (10), (11) образующая Г разбивается на N граничных элементов. Используются «постоянные» граничные элементы, т. е. элементы считаются отрезками прямых, а величина ф внутри каждого элемента считается постоянной (рис. 2). Тогда интегралы (10) и (11) преобразуются в конечные суммы:

N

g (x ) = Х ф, JG (-^4 К dr; (12)

г=1 г,.

N

h (x )=Х ф, JF (x4 kdr, (13)

г=1 г,

где Г, - ,-й граничный элемент, ф, - значение функции ф (4) на ,-м граничном

элементе.

88

М.А. Пономарева, Е.А. Собко, В.А. Якутенок

Теперь, последовательно располагая точки наблюдения х в узлах (т. е. серединах) граничных элементов и записывая для этих точек одно из уравнений (12) или (13) (в зависимости от принадлежности точки к Sj или S2), можно получить систему из N линейных алгебраических уравнений вида

N

gj = Х 6 ф, ; (14)

i=1 N

h =Ё F Ф, (15)

i =1

относительно N неизвестных ф,, где gj и hj - известные значения потенциала или потока в узле j-го граничного элемента,

Gt = i6(j■ 5Isdг. F = J F(хj.5)dr.

Г г

Xj - узел j-го граничного элемента (рис. 2), а G (х, 5) и F (х, 5) определяются выражениями (8) и (9). Значения ф, находятся из системы (14), (15). После этого потенциал и поток в любой точке х eQ могут быть найдены по формулам

NN

u (х ) = Х ф, 16 (х■5 )rsdr ^ q (х ) = Х ф, IF (х5 )rsdr.

i=1 г i=1 Г

Вычисление интегралов (7-, Fj

Из выражений (8) и (9) видно, что в интегралах Gj и Fj подынтегральные

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

Решение осесимметричных задач теории потенциала

89

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

Для вычисления интеграла Git используется координата п (рис. 3). Начало оси П располагается в узле элемента. Тогда от интегрирования по Г можно перейти к интегрированию по п в пределах от -Li до Li, где Li - половина длины i-го элемента.

Из рисунка видно, что

Г - rк = п cos а =-пег; (16)

= п sin а = nez, (17)

где r, zi - координаты точки xi (узла i-го элемента), п - координата точки £, er, ez - компоненты орта е. Выражения a, b и m, определяемые формулами (7), при x = xt с учетом (16) и (17) приобретают вид

a = 2т-2 + 2пгег + п2,

b = 2r2 + 2пгег, (18)

2b 4r2 + 4пгег

m =-----= —т2-------' т 2

a + b 4r2 + 4пгег + п2

и являются функциями переменной п.

При подстановке (18) в (8) с учетом (16) интеграл Gii записывается как

G =

L K

1!-

тт •>

f 4r2 + 4пг-ег ^ у 4r2 + 4пгег + п2

•(ler + r )п

V4r2 +4пг

ег + п

Выражение, стоящее под знаком интеграла, является произведением двух функций:

1 Ц -

Gu = - I K (п)A (п)^п =

IT J

(19)

90

М.А. Пономарева, Е.А. Собко, В.А. Якутенок

где

K (n) = K А (п) =

^ 4^2 + 4щег ^

v 4гг2 + 4щег + п2

п- + r

V4r2+4nr

2

er + п

Эллиптический интеграл (20) можно аппроксимировать многочленом [7]:

K (п)“ Ьп - 1п8г •

(20)

(21)

(22)

а функция А (п) после разложения в ряд Тейлора с учетом слагаемых порядка не выше п имеет вид

1е

А (п)~ — +—п •

2 4r

Тогда из (19) находится приближенное значение Gii :

G,

П)[^п-in8-Vп)=Lr ~2r2 n L12^ 8r) {2 4r 1)1 n

e2 L2 L

1 + -^L- - 1n -

12r

8r

(23)

(24)

i /

Интеграл Fu имеет значение 0.5 в соответствии с теоремой о разрыве производной потенциала простого слоя [8].

Решение тестовых примеров

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

формулам E

и q

.(Ч)

(А) (А)

численно полученные значения потенциала или потока в i-й точке, и) ’ и q) ’ -их точные значения, n - число рассматриваемых точек. При вычислении интегралов Gj , Fj в несингулярном случае применяются квадратурные формулы Гаусса

с четырьмя и восемью точками. Полученные результаты представлены на графиках черными и серыми линиями соответственно. Это сделано для определения оптимального количества точек интегрирования.

Пример 1. Рассматривается круговой цилиндр единичной длины и единичного радиуса. На основаниях определены потенциалы 0 и 1. На боковой поверхности задан нулевой поток (рис. 4). Точное решение для потенциала: и = z. Численный расчет проведен на трех прямых: r = 0

Решение осесимметричных задач теории потенциала

91

(ось цилиндра), r = 0.5 (потенциал вычисляется в 101 точке: zi = 0.01- [i -1], 1 < i < 101) и r = 1 (боковая поверхность, вычисление потенциала в узлах элементов). На рис. 5 приведены графики зависимости погрешности для потенциала от числа граничных элементов. Причем N означает число элементов на каждом из трех участков границы, т.е. общее число элементов - 3N. Видно, что на границе решение сходится медленнее, чем внутри области решения. Это может быть объяснено тем, что при вычислении потенциала в узлах необходимо вычислять сингулярные интегралы, которые, по-видимому, определяются с меньшей точностью, чем несингулярные. Вообще говоря, эти интегралы имеют слабые особенности порядка логарифмической и поэтому могут быть посчитаны непосредственно по квадратурной формуле Гаусса. Рис. 6 показывает незначительные отличия в результатах при вычислении этими двумя способами, что подтверждает правильность формулы (24). Также проведен расчет потока на прямой z = 1 (верхнее

Рис. 5. Зависимости погрешности E для потенциала от числа элементов N для примера 1: точки - r = 0 , сплошные линии - r = 0.5 , пунктирные линии - r = 1. Черный цвет - формулы Гаусса с четырьмя точками, серый цвет - формулы Гаусса с восемью точками

Рис. 6. Зависимости погрешности E для потенциала от числа элементов N для примера 1: точки - при использовании формулы (24), сплошные линии - при использовании формулы Гаусса. Черный цвет -формулы Гаусса с четырьмя точками, серый цвет -формулы Гаусса с восемью точками

92

М.А. Пономарева, Е.А. Собко, В.А. Якутенок

основание цилиндра, вычисление потока в узлах элементов). Значение потока для внутренних узлов (т. е. всех узлов, кроме ближайшего к оси и ближайшего к боковой поверхности) сходится к аналитическому решению, однако значение в при-осевом узле отдаляется от аналитического (рис. 7). В крайнем правом узле погрешность еще в 10-15 раз больше, чем в приосевом (погрешность в узле рассчитана по формуле E = |д(Ч^ - д(А)|). Такое поведение решения для потока обусловлено наличием особенности в угловой точке.

Рис. 7. Зависимости погрешности E для потока от числа элементов N для примера 1: сплошные линии - во внутренних узлах, пунктирные линии - в приосевом узле. Черный цвет - формулы Гаусса с четырьмя точками, серый цвет - формулы Гаусса с восемью точками

Пример 2. Рассматривается область между двумя концентрическими сферами. На внутренней, единичного радиуса, задан потенциал 0; на внешней, с радиусом 2, потенциал равен 1 (рис. 8). Точное решение для потенциала

и = 21 1 -

л1г 2 + z2

Рис. 8. Постановка задачи для примера 2

Решение осесимметричных задач теории потенциала

93

Численный расчет проведен на двух прямых: r = 0 (вычисление потенциала в 101 точке: zi = 1 + 0.01- [i -1] , 1 < i < 101 - рис. 9) и z = 0 (вычисление потенциала в 101 точке: ri = 1 + 0.01 - [i -1], 1 < i < 101 - рис. 10). На рисунках приведены графики зависимости погрешности для потенциала от числа граничных элементов. Обозначено: N - число элементов на каждой из двух полуокружностей, следовательно, общее число элементов - 2N. Видно, что решения быстро сходятся к аналитическому, как при использовании формулы Гаусса для вычисления сингулярных интегралов, так и при использовании формулы (24). Полученный результат также подтверждает достоверность формулы (24).

Рис. 9. Зависимости погрешности E для потенциала от числа элементов N для примера 2 (на прямой r = 0 ): сплошные линии - при использовании формулы (24), точки - при использовании формулы Гаусса. Черный цвет - формулы Гаусса с четырьмя точками, серый цвет - формулы Гаусса с восемью точками

Рис. 10. Зависимости погрешности E для потенциала от числа элементов N для примера 2 (на прямой z = 0 ): сплошные линии - при использовании формулы (24), точки - при использовании формулы Гаусса. Черный цвет - формулы Гаусса с четырьмя точками, серый цвет - формулы Гаусса с восемью точками

94

М.А. Пономарева, Е.А. Собко, В.А. Якутенок

Пример 3. Рассматривается область внутри сферы единичного радиуса. На нижней полусфере задан потенциал 0; на верхней - потенциал равен 1 (рис. 11). Точное решение для потенциала

1

и = — 2

1 -

Ж+1 (0)-Pk-1 (0)]-

k=1

где Pk (x) - полином Лежандра порядка k. Численный расчет проведен на прямой r = 0 (вычисление потенциала в 101 точке: zi =-1 + 0.02•[/-1], 1 <i < 101).

Рис. 11. Постановка задачи для примера 3

На рис. 12 приведены графики зависимости погрешности для потенциала от числа граничных элементов (сумма ряда для получения точного решения рассчитана с точностью 10-6). Обозначено: N - число элементов на каждой из двух

Рис. 12. Зависимости погрешности E для потенциала от числа элементов N для примера 3: сплошные линии - при использовании формулы (24), точки - при использовании формулы Гаусса. Черный цвет - формулы Гаусса с четырьмя точками, серый цвет - формулы Гаусса с восемью точками

Решение осесимметричных задач теории потенциала

95

частей границы, следовательно, общее число элементов - 2N. Видно, что решения быстро сходятся к аналитическому. Это также подтверждает достоверность формулы (24).

Заключение

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

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

ЛИТЕРАТУРА

1. Бенерджи П., Баттерфилд Р. Методы граничных элементов в прикладных науках: пер. с англ. М.: Мир, 1984.

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

1987.

3. Rui Z., Jin H., Tao L. Mechanical quadrature methods and their splitting extrapolations for solving boundary integral equations of axisymmetric Laplace mixed boundary value problems // Engineering Analysis with Boundary Elements. 2006. No. 30. P. 391-398.

4. Singh J., Gliere A., Achard J. A multipole expantion-based boundary element method for axi-symmetric potential problem // Engineering Analysis with Boundary Elements. 2009. No. 33. P. 654-660.

5. Smyrilis Y.-S., Karageorghis A. A matrix decomposition MFS algorithm for axisymmetric potential problems // Engineering Analysis with Boundary Elements. 2004. No. 28. P. 463-474.

6. Reutskiy S. The method of approximate fundamental solutions for axisymmetric problems with Laplace operator // Engineering Analysis with Boundary Elements. 2007. No. 31. P. 410-415.

7. Абрамовиц М., Стиган И. Справочник по специальным функциям с формулами, графиками и математическими таблицами: пер. с англ. М.: Наука, 1979.

8. Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1977.

Статья поступила 16.09.2015 г.

Ponomareva M. A., Sobko E. A., Yakutenok V. A. SOLVING AXISYMMETRIC POTENTIAL PROBLEMS USING THE INDIRECT BOUNDARY ELEMENT METHOD

DOI 10.17223/19988621/37/8

The paper contains all necessary relations to implement the indirect boundary element method in order to solve axisymmetric potential problems. It includes fundamental solutions of Laplace’s equation for the potential and flux in the axisymmetric case. These solutions contain complete elliptic integrals of the first and second kinds. Based on this, boundary integral equations were written corresponding to the boundary value problem. The equations were quantized by means of constant elements. The approximate formula for the integral of the fundamental solution for the potential along the element with singularity was obtained using truncated Taylor series and complete elliptic integral of the first kind approximation by a polynomial. In a similar case for the

96

М.А. Пономарева, Е.А. Собко, В.А. Якутенок

flux, a value of 0.5 was used according to the theorem about the discontinuity in the derivative of the simple layer potential. Approximate convergence of the method was explored based on three test examples. Attained results demonstrate a good convergence of the method, except for the flux computation in a close proximity to the axis of symmetry and corner points, which have nonremovable singularities. It was also shown that using the Gauss quadrature formula with four points is sufficient to estimate the nonsingular integral along the elements with a sufficient level of accuracy.

Keywords: potential theory, Laplace’s equation, axisymmetric problems, indirect boundary element method, singular integrals.

PONOMAREVA Maria Andreevna (Candidate of Physics and Mathematics, Assoc. Prof.,

Tomsk State University, Tomsk, Russian Federation)

E-mail: [email protected]

SOBKO Evgeniy Alekseevich (Tomsk State University, Tomsk, Russian Federation)

E-mail: [email protected]

YAKUTENOK Vladimir Albertovich (Doctor of Physics and Mathematics, Prof.,

Tomsk State University, Tomsk, Russian Federation)

E-mail: [email protected]

REFERENCES

1. Banerjee P.K., Butterfield R. Boundary Element Methods in Engineering Science. London, McGraw-Hill Company, 1981.

2. Brebbia C.A., Telles J. C.F., Wrobel L.C. Boundary Element Techniques. Theory and Applications in Engineering. Berlin, Springer-Verlag, 1984.

3. Rui Z., Jin H., Tao L. Mechanical quadrature methods and their splitting extrapolations for solving boundary inte-gral equations of axisymmetric Laplace mixed boundary value problems. Engineering Analysis with Boundary El-ements, 2006, no. 30, pp. 391-398.

4. Singh J., Gliere A., Achard J. A multipole expantion-based boundary element method for axi-symmetric potential problem. Engineering Analysis with Boundary Elements, 2009, no. 33, pp. 654-660.

5. Smyrilis Y.-S., Karageorghis A. A matrix decomposition MFS algorithm for axisymmetric potential problems. Engineering Analysis with Boundary Elements, 2004, no. 28, pp. 463-474.

6. Reutskiy S. The method of approximate fundamental solutions for axisymmetric problems with Laplace operator. Engineering Analysis with Boundary Elements, 2007, no. 31, pp. 410-415.

7. Abramowitz M., Stegun I. A. Handbook of mathematical functions with formulas, graphs and mathematical tables. New York, Dover, 1974.

8. Tikhonov A.N., Samarskiy A.A. Uravneniya matematicheskoy fiziki. Moskow, Nauka Publ., 1977. (in Russian)