CC BY
40
Вестник КРАУНЦ. Физ.-мат. науки. 2012. № 2 (5). C. 42-50. ISSN 2079-6641
ИНФОРМАЦИОННЫЕ И ВЫЧИСЛИТЕЛЬНЫЕ ТЕХНОЛОГИИ
УДК 533.12
ПРИМЕНЕНИЕ МОДИФИЦИРОВАННОГО МЕТОДА ГРАНИЧНЫХ ЭЛЕМЕНТОВ К РЕШЕНИЮ ЗАДАЧ ТЕОРИИ ПОТЕНЦИАЛА Федотов В.П., Горшков А.В.
Институт Машиноведения УрО РАН, 620049, г. Екатеринбург, ул. Комсомольская, 34.
E-mail: alex.gorshkov@usu.ru
Рассматривается задача теории потенциала, описываемая уравнением Лапласа А и = 0. В работе предложен алгоритм решения задачи потенциала, основанный на методе граничных элементов. Рассмотрены примеры решения задачи Дирихле для круговой области, для куба и решения задачи для куба со смешанными граничными условиями. Проводится сравнение решений, полученными численно-аналитическим методом граничных элементов с аналитическими решениями и решениями, полученными численным интегрированием по граничным элементам.
Ключевые слова: потенциал, граничные элементы, аналитические вычисления
(с) Федотов В.П., Горшков А.В., 2012
INFORMATION AND COMPUTATION TECHNOLOGIES
MSC 65R10
APPLY MODIFIED METHOD OF THE BOUNDARY ELEMENTS TO THEORY POTENTIAL PROBLEMS Fedotov V.P., Gorshkov A.V.
The Russian Academy of Sciences Ural Branch Institute of Ingineering Sciences, Ekaterinburg 620049, Komsomolskay str., 34, Russia E-mail: alex.gorshkov@usu.ru
An algorithm for solving the potential problem, based on the boundary element method, is proposed. Some example of the Dirichlet problem for the circle and the cube are consider. The numerical-analytical solution of the boundary element method were compared with the analytical solution and the numerical solution of the boundary element method.
Key words: potential problem, boundary element, analytic calculation
(c) Fedotov V.P., Gorshkov A.V., 2012
Введение
Многие физические задачи: гидравлики, электростатики, установившейся теплопроводности и фильтрации, равновесной диффузии сводятся к решению задачи потенциала. Как правило, они решаются численно с помощью сеточных методов, которые приводят к достаточно трудоемкому и не всегда корректному счету. Поэтому актуально повышение скорости и точности решения. Одно из направлений повышения скорости вычислений - уменьшение размерности системы разрешающих уравнений. Использование метода граничных элементов позволяет на единицу понизить порядок системы. Другой путь повышения скорости и точности - максимально возможное использование аналитических вычислений. В известных работах [1, 2, 4] интегралы по граничным элементам берутся численно, с использованием квадратурных формул Гаусса. В трехмерном случае это плоские элементы, как правило,
треугольные. Поэтому другой путь повышения скорости и точности - максимально возможное использование аналитических вычислений. Это позволяет получить точные значения интегралов по элементам и их предельные выражения. Такой подход использовался для решения плоских параболических задач [5, 6] и плоских задач теории упругости [7, 8].
Основное уравнение
Рассматривается задача теории потенциала, описываемая уравнением Лапласа:
Ч)и = 0 (1)
в области V с поверхностью Г. На поверхности заданы граничные условия:
при х е Гх, и = и(х); (2)
при х е Г2, q = ц(х), (3)
где Гх и Г2 - части поверхности, Г^Г2 = Г, ГхПГ2 = 0, чертой отмечены заданные функции, q = ди/дп, п - вектор единичной нормали к поверхности, nj его
компоненты. Для решения используется известный подход взвешенных невязок [1]:
IУ2и(х)и*($,х)дУх + 1 (и - й)О*С,х)ёГх- (4)
V Гх
-| ^ - д)и *(£, х)йГх = 0.
Г2
Здесь и*(£,х) - весовая функция, ,х) = ди,х)/дп = ^VjU*. Преобразуем
первое слагаемое в (4), взяв интеграл дважды по частям. Индекс х в дУх и ^Гх показывает, что интеграл берется по переменным х:
I V2ju(x)U*($,х)дУх = I (Vju(x)U*($,х))^дУх - I Vju(x)VjU*($,х)дУх.
V V V
К первому слагаемому правой части применим формулу Гаусса - Остроградского, а интеграл по объему снова проинтегрируем по частям:
У Vju(x)njU*(%,х^х - ! Vju(x)VjU*(%,х^х =
Г V
= 1 q(x)U*($,х)йГх-| (u(x)VjU*($,х)),jdVx +1и(х^^*($,х)д^х.
Г V V
Применяя ко второму слагаемому формулу Гаусса - Остроградского получим выражение:
Iq(x)U*($,х^Гх - I и(х)й*{:,х^Гх + Iи(хУ^р*($,х^х.
ГГ V
После подстановки в (4) получим:
Iи(х^р*(:,х^х - Iq(x)U*(:,х^Гх (5)
V Гх
+ 1 *($, х^Гх + ! й(хт:, х^Гх-] и(х)а*С, х^Гх = 0.
Г2 Г2 Гх
Выберем весовую функцию, удовлетворяющую уравнению
V2U * = -2ап8 £, х),
где 8С,х) - дельта-функция Дирака, для плоских задач а = х, для трехмерных а = 2. Решение этого уравнения в случае трехмерных задач имеет вид:
u*(:,х) = х Л,
4п г(:, х)
где г = а/Сх -хх)2 + (:2 -х2)2 + Сз -хз)2 в плоском случае
и *(:, х) = - 2П1п г(:, x),
г = а/Сх - хх)2 + С2 - х2)2 . В результате подстановки в (5) получим граничное интегральное уравнение [1]:
и(: )= ! О*(:, x)u(x)dГx - ! U *(:, x)q(x)dГx. (6)
ГГ
Здесь интегралы берутся по всей поверхности, включая области, где граничные условия известны и неизвестны.
Согласно методу граничных элементов, разобьем границу на элементов. Удобнее использовать треугольные элементы, а в плоском случае отрезки прямых. Для аппроксимации неизвестных функций выберем наиболее простой вариант - потенциал и(х) и поток q(x) на элементе постоянны. Будем считать, что эти значения соответствуют узлу ^ элемента и равны щ и qk соответственно. В качестве узла ^ выберем середину элемента, в трехмерном случае точку пресечения медиан треугольника, а в плоском - середину отрезка. Положим, что поверхность Гх содержит N элементов, которые обозначим Г(х),Г(2),...Г(^ а поверхность Г2 - М элементов, которые обозначим Г*^+х),Г^+2),...Г^+М), N + М = п. Для узла £р, лежащего на элементе Г(р) запишем:
— ир = I qг [ и (: p, х^ Гх иг I q (:p, х¥Гх | +
г=х V Г(г) Г(г) / (7)
N+M ( С Г \
+ £ Ьг]и*Ср, х^Гх - иг] р, х^Гх\, р = х, 2... N + М
г(г) г(г)
Уравнение (6) свелось к системе линейных уравнений (7). Для решения задачи нужно вычислить интегралы:
Umn — / и {^m, X)dГx, Qmn — / q {^m, x)dГ;
x*
(8)
Интегралы (8) вычислены аналитически. В плоском случае эти интегралы представляют собой линейные комбинации выражений:
Здесь Ь - длина отрезка. Как видно из приведенных выражений, интегралы зависят от логарифмов расстояний от точки £ до концов элемента и от углов между элементом и направлением на точку £.
Пример 1
Стационарная тепловая задача для круговой области при заданных на границе значениях температуры (задача Дирихле для уравнения Лапласа в круге).
Рассмотрим задачу определения температуры в круговой области с радиусом Я при заданных на границе значениях температуры:
Для сравнения найдем приближенное решение с помощью модифицированного метода граничных элементов. Граница области - окружность радиуса Я = х - приближается ломаной, состоящей из п отрезков прямой одинаковой длины. Температура и поток полагаются постоянными на каждом отрезке. Чем больше п, тем точнее будет приближение окружности ломанной. Результаты расчетов. На рис. 1 представлены графики решений при различном числе узлов. Как видно из результатов решений наблюдается очень быстрая численная сходимость метода, что объясняется аналитическими вычислениями в основе общей процедуры счета, причем, уже при разбиении на 32 граничных элемента отличия от аналитического решения практически отсутствуют, что позволяет говорить о снижение порядка разрешающей системы алгебраических уравнений.
Q1 — ln (^12 + £2), q2 — In ^(£1 — L)2 + £|j ,
(9)
кАи — 0,0 < ф < 2п, 0 < r < R,
и — 1 + sinф + 2 sin3^ + cos4^,0 < ф < 2п, r — R. Задача имеет аналитическое решение
Гз
и — 1 + rsinф + ~2 sinЗф + r4cos4ф.
Рис. 1. Графики решений при различном числе узлов: п = 8 (а); п =16 (Ь); п = 32 (с); для сравнения график аналитического решения (й)
Трехмерный случай
Для вычисления интегралов в трехмерном случае вершины элемента нумеруются так, чтобы обход со стороны внешней нормали поверхности был против часовой стрелки. На каждом элементе вводится локальная система координат ОХ1Х2Х3 . Систему выбираем так, что начало о совпадает с третьей вершиной элемента. Ось охз параллельна внешней нормали поверхности и направлена в ту же сторону. Ось ох2 параллельна стороне, противоположной вершине треугольника. Ось 0x1, соответственно перпендикулярна этой стороне (см. рис. 2). Координаты вершин элемента в
Рис. 2. Локальная система координат элемента
локальной системе координат равны (0,0), (h,k1h), (h,k2h), где k1 = tan<p1, k2 = tan(p2, <Ръ Ф2 - углы, которые образуют стороны oa1 и oa2 элемента с осью. Такой выбор системы координат обеспечивает определенную симметрию интегрируемых выражений. Координаты точки наблюдения в локальной системе координат элемента m
обозначим (£тЬ £т2, £тз). Вычисляемые интегралы примут вид:
1
итп —
г») VЙз+(хі- Сіті)2+(х2- Ст^)2
^гх,
бтп —
г(п)
[Стіз+(хі- Стті)2+(х2- Стт2)2]
3/2
^г*
В результате перехода к повторным получим:
к &2*1
итп —
0 к1х1 \/^Іт3 + (х1 - Стті)2 + (х2 - СІ2)2
к /г2х1
-Лх2 (1х1.
йтп —
0 ^1X1
Кз + (х1- Стт1)2+(х2- Стт2)2]3/2
(1х2(1х1.
(10)
(11)
Интегралы (10, 11) после ряда замен приводятся к дробно-рациональным выражениям и интегрируются в элементарных функциях. Обозначим первообразную интеграла (10) через *1(с, £3,г,к), а интеграла (11) через ^2(^1,гвга2, £3,г,к). Функция *1(с, Сз,г, к) представляет комбинацию функций задающими угол между элементом и направлением на точку влияния и логарифм расстояния от узла элемента до точки влияния типа (9). Для функции ^2(^1, ^2, Сз,г,к) - комбинацию функций с аналогичными аргументами. Рассматриваются следующие взаимные положения точки £ и элемента:
1. Точка £ лежит вне плоскости элемента. Тогда подынтегральные функции непрерывны и интегралы (10, 11) непрерывные функции параметров - координат точки £.
2. Точка £ лежит на плоскости элемента, внутри или вне его. Пределы первообразных соответственно равны:
^(с, 0, х, к) — —
2с 1п
л/і2 + с2 — 7
(у1 + к2 — к) V! + к2
ск 1п
2 + 2 (с— (к+ V1 +к2 ) і)(і— V і2+с2 )
с + ТШ2
+11п
с—к7 + Л/72 I с2
VШ2 + Л/ 7 + с .
— 7
>/1 + к2
^2(с, 0, х, к) — 2| яг'^п с ^к + \/1 + к2 ^ + х — л/ с2 + х2
■\/1 + к2 — к) — х + л/ с2 + х2
+
+яг'£п
Графики предельных функций и ^2 представлены на рис. 3а и 3Ь.
3. Точка лежит на границе элемента. Функции 2ши(Съ С2,0) имеют разрыв равный 2п, а функции итп(^1, £2, С3) - устранимые разрывы.
1
1
Пример 2
Для оценки точности построенного метода решена задача Дирихле для уравнения Лапласа в кубе. Заданы следующие граничные условия: поверхность Г - вся поверхность куба, поверхность Г2 = 0.
При
X = -0.5, -0.5 < х2 < 0.5, -0.5 < х3 < 0.5, и = 2;
X = 0.5, -0.5 < х2 < 0.5, -0.5 < х3 < 0.5, и = 1;
-0.5 < х1 < 0.5,х2 ± 0.5,-0.5 < х3 < 0.5, и = 0;
-0.5 < х1 < 0.5, -0.5 < х2 < 0.5,х3 ± 0.5, и = 0,
полученное численное решение сравнивалось с аналитическим, построенным в виде ряда
8 ^ (—1)k+mcosп(2k — 1)x2cos п(2m — 1)x3
г2
и=
^m—l
(2к — l)(2m — 1)
-X
(l2)
Х (3cOSh(ClmkXl) — SІnh(ClmkXl)) ,
где clm = (2к — l)2 + (2m — l)2. На рис. 4a показано сравнение численного реше-
Рис. 4. Сравнение аналитического и приближенного решений
ния с аналитическим на прямой (—0.5,0,0) — (0.5,0,0). Сплошной линией выделено
Применение модифицированного метода граничных элементов ... ISSN 2079-6641
численное решение, пунктиром - аналитическое. На рис. 4Ь показана погрешность численного решения вдоль той же прямой.
Было также построено решение методом граничных элементов с численным интегрированием по элементам. Сравнительные результаты представлены в таблице.
Таблица
Сравнение решений, полученных различными методами
X Аналитическое решение Численно- аналитическое решение Решение по методу граничных элементов
-0.375 1.43016 1.50863 1.34032
-0.250 0.96719 1.0412 0.801924
-0.125 0.65859 0.703483 0.484194
0.0 0.5 0.523011 0.337388
0.125 0.4721 0.487691 0.327176
0.250 0.56012 0.578455 0.440945
0.375 0.74827 0.763672 0.684603
Из приведенной таблицы и графиков видно, что погрешность численно-аналитического решения не превышает 10 %. Погрешность численного решения существенно больше. Время счета численного решения примерно в 3 раза больше, чем численноаналитического.
Пример З
Для этой же области решена задача для уравнения Лапласа с иными краевыми условиями: поверхность Г2 образована двумя квадратами:
х1 = —O.5, —O.5 < x2 < O.5, —O.5 < x3 < O.5, и = 2;
Xl = O.5, —O.5 < X2 < O.5, —O.5 < хз < O.5, и = 1;
поверхность Гі:
—O.5 < х1 < O.5,х2 ± O.5, —O.5 < хз < O.5, и = O;
—O.5 < х1 < O.5, —O.5 < х2 < O.5,хз ± O.5, и = O.
Аналитическое решение при данных краевых условиях имеет вид:
4 ^ (—l)k+mcosп(2к — l)x2cosп(2m — 1)хз sinh(c1mkx1)
и п2 (2к — l)(2m — 1) cosh(clmkXl) (13)
На рис. 5a показано сравнение численного решения с аналитическим на прямой (-0.5, 0, 0) - (0.5, 0, 0). Сплошной линией выделено численное решение, пунктиром
- аналитическое. На рис. 5b показана погрешность численного решения вдоль той же прямой.
“ 4. чх из "■■■. ■% XV 1 Ф.Щ ’ а.ш.
-il4 -Q2 -0.L -0.2 0:2 С 4 -її 2 У - -G.QL k-
Ь
Рис. 5. Сравнение аналитического и приближенного решения
Заключение
В работе рассматривается применение модифицированного метода граничных элементов к некоторым задачам теории потенциала. В изложенном методе максимально используются аналитические вычисления. Их результаты годятся для использования в любых задачах теории потенциала. Основой метода является вычисление интегралов по базовому элементу.
Приведенные примеры решения некоторых задач показывают, что используемый метод позволяет получить достаточно высокую точность при небольшом числе элементов, а скорость вычислений больше, чем при использовании численных методов интгрирования по элементам.
Библиографический список
1. Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов. М.: Мир, 1987. 524 с.
2. Хомяков А.Н. Метод граничных элементов повышенной точности в задачах гидродинамики идеальной жидкости // Вычислительные методы и программирование. 2008. Т. 9. С. 401-404.
3. Михлин С.Г. Лекции по линейным интегральным уравнениям. М.: ГИФМЛ, 1959. 232 с.
4. Зарубин В.С., Кувыркин Г.Н. О возможностях метода граничных элементов при моделировании континуальных систем. - иИЬ: http://technomag.edu.ru/doc/48397.html
5. Федотов В.П., Нефедова О.А. Применение модифицированного метода граничных элементов для решения параболических задач // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. 2011. Вып. 4 (25). С. 93-101.
6. Федотов В.П. Модифицированный метод граничных элементов для решения задач колебания пластин // Вестник КРАУНЦ. Физ.-мат. науки. 2011. № 2 (3). С. 18-32.
7. Федотов В.П., Спевак Л.Ф. Аналитическое интегрирование функций влияния для решения задач упругости и теории потенциала методом граничных элементов // Математическое моделирование. 2007. Т. 19. № 2. С. 87-104.
8. Федотов В.П., Горшков А.В. Численно-аналитический метод решения задач упругости с особенностями // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. 2005. Т. 38. С. 29-34.
Поступила в редакцию / Original article submitted: 08.11.2012
