Программные продукты и системы /Software & Systems
№ 4 (108), 2014
УДК 519.67 Дата подачи статьи: 05.06.2014
DOI: 10.15827/0236-235X.108.178-182
ПРОГРАММНЫЙ КОМПЛЕКС ДЛЯ РЕШЕНИЯ ЗАДАЧ ТЕОРИИ ПОТЕНЦИАЛА МЕТОДОМ ГРАНИЧНЫХ ЭЛЕМЕНТОВ
(Работа выполнена при поддержке РЦПразвития вычислительных, телекоммуникационных и информационных
ресурсов УрО РАН - РЦП-2014, проект № РЦП-14-П1)
В.П. Федотов, д.т.н., старший научный сотрудник, главный научный сотрудник, fedotov_vp@mail.ru;
Л.Ф. Спевак, к.т.н., доцент, зав. лабораторией, lfs@mach.uran.ru;
О.А. Нефедова, к.т.н., научный сотрудник, nefedova@mach.uran.ru (Институт машиноведения Уральского отделения РАН, ул. Комсомольская, 34, г. Екатеринбург, 620049, Россия)
Работа посвящена развитию подхода к реализации метода граничных элементов, направленного на распараллеливание вычислений, для решения двухмерных задач об установившихся потенциальных течениях. Главной особенностью применяемых алгоритмов решения является точное вычисление всех интегралов по граничным элементам с помощью выведенных авторами аналитических формул. Это позволяет заметно повысить точность решения и сократить время расчета. Разработанный программный комплекс предназначен для решения двухмерных задач теории потенциала в области произвольной геометрии при заданных граничных условиях для искомой функции или потока. В комплекс заложена возможность решения однородных задач (при отсутствии внутренних источников), а также задач с заданными функциями источников. Функция источника задается отдельным программным модулем. Особо рассмотрен случай, когда функция источника является гармонической. Это позволяет свести все вычисления при решении неоднородной задачи на границу области. Для задания геометрии исследуемой области и ввода параметров задачи и граничных условий в программный комплекс включен графический редактор. Область задается своей внешней границей, состоящей из прямолинейных и круговых участков. Для каждого участка задаются граничные условия и количество граничных элементов на нем. С помощью графического редактора внутри расчетной области может быть задана зона, в которой требуется вычислить значения искомой функции. Эта зона также задается границей из прямолинейных и круговых участков. Частота расчетной сетки в зоне задается произвольно, пересчет для новой сетки не требует нового решения исходной задачи. Программный комплекс был реализован на суперкомпьютере «Уран» ИММ УрО РАН. В качестве примеров рассмотрены задачи распространения тепла в пластинах квадратной и эллиптической форм при различных функциях источника.
Ключевые слова: параллельные вычисления, стационарные задачи теории потенциала, функция источника, метод граничных элементов, аналитическое интегрирование.
Все большая доступность супервычислителей в настоящее время повышает интерес к методам решения краевых задач, допускающим распараллеливание счета на уровне алгоритма. В связи с этим разработка подобных алгоритмов и их программная реализация на многопроцессорных комплексах находят все большее распространение и применение на практике. В данной работе представлены алгоритм и комплекс программ для решения двухмерных задач теории потенциала.
В основу алгоритма положена методика применения метода граничных элементов [1], основанная на точном вычислении интегралов по граничным элементам с помощью полученных авторами аналитических формул и параллельных вычислений на каждом этапе решения. Разработанный подход и созданные на его основе программы показали свою эффективность для решения линейных задач теории упругости [1], а также нестационарных задач диффузии и теплопроводности [2].
Краевая задача потенциального течения
Основным дифференциальным уравнением задачи теории потенциала является уравнение Пуассона, в двухмерном случае имеющее вид
Д9(х) = b(x), x(xb x2) е Q, (1)
где 0(x) - значение искомой функции в точке x; A=S2/5xi2+S2/5x22 - оператор Лапласа; Q - рассматриваемая область; b(x) - заданная функция источника, распределенная на области Q.
Стационарное потенциальное уравнение (1) описывает многие физические процессы, такие как кручение призматического стержня, установившееся течение жидкости, а также стационарное тепловое поле, электростатическое поле и т.д. Задача заключается в определении функции 0(x) (например, электрического или гидравлического потенциала либо температуры), удовлетворяющей уравнению (1) и граничным условиям:
9(x0) = 9*(x0), x0(x01, x02) е Sj; q(x0) = q\x0), x0(x01, x02) е S2. (2)
Здесь q(x)=d0(x)/dn(x) - поток (соответственно плотность электрического тока, скорость течения жидкости или поток тепла); n(x) - внешняя нормаль к границе S = S1u S2 области Q; звездочкой отмечены известные величины.
Алгоритм решения
В соответствии с методом граничных элементов [3] осуществляется переход от основных урав-
178
Программные продукты и системы /Software & Systems
№ 4 (108), 2014
нений краевой задачи (1) и (2) к граничному интегральному уравнению для граничной точки x0:
19( x0) = J q (x) u (x0, x) dS (x ) —
2 / . (3)
-J9(x) f (x,x)dS(x) — Jb (x)u (x,x)dQ(x),
S Q
где u (x0, x), f (x0, x) - функции влияния. Для двухмерной задачи теории потенциала
u'(м x) = Ь“ (1);
f * ( x0.x ) =
2 п
du* (x,x) dn ( x)
1
2nr
rm = Xm — x0m . (4)
Здесь и далее по — овторяющемуся индпксу производится суммирование от 1 до 2.
Для решения граничного интегрального уравнения необходима дискретизация границы S, то есть построение модели области Q с границей, разбитой на граничные элементы. Пусть N + M -общее количество граничных элементов в предположении, что часть границы Sj разбита на N граничных элементов, а часть границы S2 - на M элементов. В данной работе используются прямолинейные граничные элементы и принимается, что функции 0(x) и q(x) имеют постоянные значения на каждом элементе
9(x;) = 9; = const; q(x) = q, = const, i = J, 2, ..., N+M. (5)
Здесь x, - узел, находящийся в середине элемента [a;-J, a,].
Уравнение (3) в этом случае приводится к виду
в . . N+M . .
-9(xj)= Ё qt J u*(x,x)dS(x)—
? (x.
2 t=в [al—i,a, ]
N+M
(6)
— Ё9,. J f (xj, x) dS (x) —
t=i [^t ]
—J b (x) u (x, x) dQ(x), j=1, 2, ..., N + M.
Q
Вычисление интеграла по области (последнее слагаемое в правой части уравнения (6)) зависит от вида функции источника b(x). Рассмотрим два случая.
J. b(x) - гармоническая функция в области Q, то есть имеет место равенство Ab(x) = 0. Тогда согласно [3] интеграл по области в правой части (6) преобразуется в интеграл по границе области следующего вида:
J b (x) u* (x, x) d Q( x) =
Q
= IT l|b ( x )f21nB + 1] rmnm — r 2 fln1 + 1]тТГУу? ( x )■
8nJ| 4 /f r ) m m ( r )dn(x
2. b(x) - произвольная числовая функция. В этом случае область Q разбивается на конечные элементы AQJ, Л^2, ., AQK, в каждом из которых функция b(x) полагается постоянной. Тогда интеграл по области в правой части уравнения (6) вы-
числяется следующим образом:
Jb(x)u*(xJ,x)dQ(x) = ^b(x) J M*(x,x)dQ(x). (7b)
Q k=1 AQ
Здесь xk - узел элемента AQk.
С учетом (7a) и (7b) уравнение (6) принимает соответственно вид
в . . N+M . .
-9(xj)= Ё qt J u*(xj,x)dS(x)—
2
i=1 [ai—i,at ]
N +M
— Ё9, J f * (xj,x)dS (x) —
i=1 [a,.! ,a, ]
1 N+M ( ( 1
—— Ё f у (уУ 2У11 +1
8л Ё п П П ( \ r
1 )db (x) У у у
-+11—И Уу у x)
r ) dn (x )J V ’
rmnm — (8a)
или
1 , , N+M , ,
T9(xj)= Ё qt J u*(xj,x)dS(x)—
2 i=1 [al—1,a, ]
N+M , ,
— Ё 9 J u* (xj, x) dS (x ) —
i=1 [a,—1,a, ]
(8b)
—Ё b (x) J u* (x, x)dQ(x).
k=1 AQk
Интегралы по граничным элементам от функций влияния вычисляются аналитически по формулам, полученным в [1] с помощью перехода от произвольного граничного элемента к базовому отрезку [0, L], лежащему на оси абсцисс:
L ___
J u* (x, x) dS (x) = J u (xj, x) dx1 ;
[a—1, ]
J f * (xj, x) dS (x) = J f * (xj, x) d^. (9)
K-t, ] 0
Соотношения (8), записанные для всех граничных узлов, образуют систему линейных алгебраических уравнений относительно узловых значений функций 0(x) и q(x), не заданных граничными условиями:
P*0 + R*q + b* = 0. (J0)
Здесь P и R - квадратные матрицы размером (N + M)X(N + M), элементы которых p и r( являются интегралами от функций влияния по граничным элементам [a;-J, a,]; 0 - (N + М)-мерный вектор с элементами 0, - значениями функции 0(x) в узлах, причем первые N значений заданы граничными условиями, а последующие M неизвестны; q - (N + М)-мерный вектор с элементами q, -значениями потока q(x) в узлах границы, при этом первые N значений неизвестны, а последующие M заданы граничными условиями; b* - вектор размерности (N+M), элементы которого b( =
= J b ( x) u ( x, x) dQ( x)
определяются численным
интегрированием в соответствии с (8a) или (8b).
r n
Q
J79
Программные продукты и системы /Software & Systems
№ 4 (108), 2014
После определения всех граничных значений значения функции 0 во внутренних точках области рассчитываются по формуле
N+M
0(^)= Ё q J «* *)dS(*)-
'=! \a,,,a,]
11 (11)
N+M v '
- Ё J /*(^,*)dS(x)-\b(*)«*(£,,*)dQ(*).
1=1 \aM ,at ] Q
Программная реализация
Представленный алгоритм был реализован в виде программного модуля. Языки программирования выбирались с требованием максимальной переносимости комплекта программ под различные операционные системы. Для создания графического интерфейса был использован язык программирования Java. Java-модуль ru.uran.imach.po-tential состоит из одного пакета и включает в себя шесть базовых и несколько вспомогательных классов. Назначение базовых классов:
- класс FrontFrame отвечает за главное окно программы, описывает вид и стиль программы, панель инструментов и вспомогательные диалоговые окна, отвечает за многооконный интерфейс в целом;
- класс Expoint2D описывает узловую точку;
- класс BorderModel отвечает за построение внешней и внутренней границ модели, содержит методы добавления, перемещения, удаления узловых точек и метод перемещения модели, включает в себя методы сохранения объекта и загрузки объекта из файла;
- класс MainPanel описывает графическое представление объекта, содержит методы отри-совки отдельных частей модели и работы с масштабом;
- класс PieceSegment описывает конкретную границу модели и произвольный граничный элемент;
- класс MGEModel представляет математическую модель объекта, таблицы с результатами расчетов, содержит возможность сохранения результатов в разных форматах на жестком диске.
Для создания программы использовалась среда разработки Eclipse. Модуль, написанный на языке Java, был скомпилирован и экспортирован в jar-архив. Программа компилировалась и выполнялась в среде Oracle Java SE 7 [4]. На рисунке 1 представлен экран ввода данных графического интерфейса.
Расчетный модуль программы написан на языке программирования C++ с использованием библиотеки параллельных вычислений MPI [5, 6] и включает в себя непосредственно алгоритм расчета. С++-модуль mge.potential состоит из четырех базовых классов:
- класс PieceSegment описывает граничные элементы объекта;
- класс Triangle содержит методы разбиения внутренней области объекта на конечные элементы;
- класс MGEModel описывает математическую модель объекта, содержит методы загрузки объекта из файла и методы разбиения границы объекта;
Рис. 1. Экран ввода данных Fig. 1. Data entry screen
180
Программные продукты и системы /Software & Systems
№ 4 (108), 2014
- класс PotentialSolver решает задачу теории потенциала.
Программный код модуля на C++ после компиляции и проверки на машине разработчика копировался на кластер umt («Уран») ИММ УрО РАН, где компилировался вновь, уже в окружении библиотеки MPI-кластера. При запуске модуля задавалось количество вычислительных ядер на кластере, кратное числу точек разбиения границы области. Выходные данные были сформированы в виде таблицы Excel. При написании модуля использовались также библиотеки BOOST [7] и GSL [8]. Программный модуль компилировался компилятором GCC версии 4.6.2. На рисунке 2 показана общая структура программного модуля.
Java-модуль
“ru.uran.imach. potential”
f Программный \ модуль
«C++ расчетный модуль “mge.potential”
Рис. 2. Общая структура программного модуля Fig. 2. Software module general structure
Примеры
Для иллюстрации работы разработанной программы приведем результаты решения двух тестовых задач.
Пример 1. Рассмотрим задачу о распространении тепла в пластине, имеющей форму эллипса с полуосями a = 5 м и d = 3 м. Внутри пластины действует источник тепла, заданный в виде гармонической функции b(x)=6x-l+6x2+2. На границе задан тепловой поток
q (х)|
4+4=1
27 х3 + 75х23 + 50х2: ^8Ц2 + 625х2
В результате решения задачи с помощью созданного пакета программ было определено распределение температуры на границе пластины и вычислены значения температуры во внутренних точках пластины. На рисунке 3 показаны результаты сравнения значений температуры вдоль оси абсцисс, рассчитанных для 20 и 100 граничных элементов, с точным решением 0(x)=xj3+x23+x22.
Пример 2. Рассмотрим задачу о распространении тепла в квадратной пластине со стороной
h=1 м. Пусть на вертикальных участках границы пластины задана температура 9* (х) ^ =100 °С,
0* (х) ^ =200 °С, а на горизонтальных - нулевой поток: q* (х) ^ = 0, q* (х)| ^ = 0, соответствующий отсутствию теплообмена с окружающей средой. Внутри пластины действует источник тепла b(x)=-103x12. Задача была решена для 4 и 20 граничных элементов. На рисунке 4 приводится сравнение изменения полученных значений температуры вдоль отрезка AB, концы которого имеют координаты (0.5,0) и (0.5,1), с точным решением 0(x)=100+183.3x1-83.3x14.
Рис. 3. Распространение тепла в эллипсе: изменение температуры вдоль оси абсцисс
Fig. 3. Heat propagation in elipsys: temperature change along an X line
Рис. 4. Распространение тепла в квадратной пластине: изменение температуры вдоль отрезка AB
Fig. 4. Heat propagation in square platten: temperature change along AB jog
Проведенные расчеты показали, что решения с помощью разработанного программного комплекса рассмотренных примерных задач с увеличением количества граничных элементов приближаются к соответствующим точным решениям. Реализованный в программном комплексе алгоритм позволяет сократить время расчета задачи по сравнению с классическим методом граничных
181
Программные продукты и системы /Software & Systems
№ 4 (108), 2014
элементов за счет использования аналитического интегрирования и технологий параллельного программирования.
Литература
1. Федотов В.П., Спевак Л.Ф. Модифицированный метод граничных элементов в задачах механики, теплопроводности и диффузии. Екатеринбург: Изд-во УрО РАН, 2009. 164 с.
2. Федотов В.П., Нефедова О.А. Применение модифицированного метода граничных элементов для решения параболических задач // Вестник СамГТУ: Сер. Физ.-мат. науки. Самара, 2011. Т. 25. № 4. С. 93-101.
3. Бреббия К., Теллес Ж., Вроубел Л. Метод граничных элементов. М.: Мир, 1987. 526 с.
4. Oracle Technology Network for Java Developers. URL: http://www.orade.com/technetwork/java/index.html (дата обращения: 11.12.2013).
5. Pacheco P. An Introduction to Parallel Programming. Morgan Kaufmann, 2011.
6. The Message Passing Interface (MPI) standard. URL: http://www.mcs.anl.gov/research/projects/mpi/ (дата обращения: 22.04.2013).
7. Boost C++ Libraries. URL: http://www.boost.org (дата обращения: 09.12.2013).
8. GSL-GNU Scientific Library. URL: http://www.gnu.org/ software/gsl (дата обращения: 26.11.2013).
DOI: 10.15827/0236-235X.108.178-182 Received 5.06.2014
A SOFTWARE PACKAGE DESIGNED TO SOLVE THE POTENTIAL THEORY PROBLEMS BY THE BOUNDARY ELEMENT METHOD Fedotov V.P., Dr.Sc. (Engineering), Senior Researcher, Chief Researcher, fedotov_vp@mail.ru;
Spevak L.F., Ph.D. (Engineering), Associate Professor, Head of Laboratory, lfs@imach.uran.ru;
Nefedova O.A., Ph.D. (Engineering), Research Associate, nefedova@imach.uran.ru (Institute of Engineering Science of the Ural Branch of the RAS,
Komsomolskaya St. 34, Ekaterinburg, 620049, Russian Federation)
Abstract. The paper develops the approach to the implementation of the boundary element method aimed at the paralleling of computations for solving two-dimensional problems on steady-state potential flows. The key feature of the applied solution algorithms is the exact computation of all integrals over boundary elements by analytical formulae derived by the authors. This noticeably increases solution accuracy and reduces calculation time. The designed software package is intended for solving two-dimensional problems of the potential theory in arbitrary geometry with specified boundary conditions for a required function or flow. A possibility of solving homogeneous problems (without internal sources) is implied, as well as problems with specified source functions. A source function is specified by a separate program module. Special attention is given to the case when the source function is harmonic. This enables all the calculations in solving an inhomogeneous problem to be brought to the area boundary. A graphic editor is included in the software package to specify the geometry of the area under study and to introduce problem parameters and boundary conditions. An area is specified by its external boundary consisting of rectilinear and circular portions. Boundary conditions and the number of boundary elements are specified for each portion. Besides, by means of the graphic editor, inside the calculated area a zone can be specified where the values of the required mesh are to be calculated. This zone is also specified by a boundary consisting of rectangular and circular portions. The mesh fineness in the zone is specified arbitrarily, and the recalculation for a new mesh does not require a new solution to the initial problem.
The software package was implemented on the Uran supercomputer installed at IMM, UB RAS. Problems on heat distribution in square and elliptic plates for various source functions are considered as examples.
Keywords: parallel calculations, stationary problems of the potential theory, source function, boundary element method, analytical integration.
References
1. Fedotov V.P., Spevak L.F. Modifitsirovanny metod granichnykh elementov v zadachakh mekhaniki, teploprovodnosti i diffuzii [A Modified Boundary Element Method in Mechanics, Transcalency and Diffusion Problems]. Ekaterinburg, Ural Branch of RAS Publ., 2009, 164 p.
2. Fedotov V.P., Nefedova O.A. Vestn. Samarskogo gos. tekh. univ. [The Bulletin of the Samara State Technical Univ.]. 2011, vol. 25, no. 4, pp. 93-101 (in Russ.).
3. Brebbia C.A., Telles J.C.F., Wrobel L.C. Boundary Element Techniques. Springer, Berlin, Heidelberg Publ., 1984, 526 p.
4. Oracle Technology Network for Java Developers. Available at: http://www.oracle.com/technetwork/java/in-dex.html (accessed December 11, 2013).
5. Pacheco P. An Introduction to Parallel Programming. Morgan Kaufmann Publ., 2011.
6. The Message Passing Interface (MPI) standard. Available at: http://www.mcs.anl.gov/research/projects/mpi/ (accessed April 22, 2013).
7. Boost C++ Libraries. Available at: http://www.boost.org (accessed December 9, 2013).
8. GSL-GNU Scientific Library. Available at: http://www.gnu.org/software/gsl (accessed November 26, 2013).
182