CC BY
64
УДК 539.3
Д. Ю. Полянский
ТРИАНГУЛЯЦИЯ ПЛОСКИХ ОБЛАСТЕЙ РЕШЕНИЕМ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ В ФОРМЕ ГАЛЕРКИНА ЗАДАЧИ ДИРИХЛЕ
Аннотация. Решена новая задача триангуляции достаточно регулярной ограниченной плоской области. Область представляется объединением непересе-кающихся выпуклых криволинейных подобластей с числом углов от 3 до 6. Каждой подобласти ставится в соответствие эквивалентный аналог - равносторонний треугольник или порождаемый им выпуклый многоугольник. Эквивалентный аналог преобразуется в дискретный, состоящий из равносторонних треугольников. Производится конформное отображение дискретного аналога на исследуемую область решением методом конечных элементов в форме Галеркина краевой задачи Дирихле с использованием Лапласиана. Результатом отображения является дискретная модель исследуемой области с треугольными элементами, близкими к равносторонним.
Ключевые слова: триангуляция, метод конечных элементов, конформное отображение, краевая задача Дирихле.
Abstract. The author has solved a new problem of triangulation of sufficiently regular bounded flat domain. The domain is represented by a combination of nonintersecting convex curvilinear subdomains with the number of angles from 3 to 6. Each subdomain is given a corresponding equivalent analog - an equilateral triangle or a convex polygon generated by it. The equivalent analog is transformed into a discrete one consisting of equilateral triangles. A conformal image of the discrete analog on the domain under analysis is made by the FEA solution in Galyorkin form of the boundary value Dirichlet problem using Laplacian. The result of the imaging is a discrete model of the domain under investigation with the triangular elements being close to equilateral.
Key words: triangulation, method of final elements, conformal image, boundary value Dirichlet problem.
Решение плоских контактных задач термо-упруго-пластичности предъявляет повышенные требования к качеству конечно-элементной сетки (КЭС). Оптимизация формы конечных элементов снижает размерность получаемых в ходе решения матричных уравнений, а значит, значительно увеличивает скорость решения и уменьшает потребности в оперативной и дисковой памяти. Это дает возможность получить решение без применения дополнительных приемов, например параллельных вычислений, используя лишь собственные возможности персональных компьютеров.
Кусочную аппроксимацию исследуемых функций часто строят на треугольных конечных элементах (ТКЭ). Установлено [1-5], что оптимальным видом ТКЭ является равносторонний треугольник. В данной работе предлагается алгоритм триангуляции достаточно регулярной ограниченной области Q на плоскости xy, состоящий из следующих шагов.
Шаг 1. Область Q представляется объединением непересекающихся выпуклых криволинейных подобластей Q с числом углов от 3 до 6.
Шаг 2. Каждой подобласти П, ставится в соответствие эквивалентный аналог ю, на плоскости pq - равносторонний треугольник или порождаемый им выпуклый п -угольник (п = 4,5,6).
Шаг 3. Строится дискретная модель Ю,, состоящая из равносторонних ТКЭ эквивалентного аналога ю,.
Шаг 4. Производится конформное отображение дискретной модели ю, на подобласть П,, сохраняющего углы на ху между пересекающимися кривыми 1}, образующими к -й криволинейный треугольный элемент Т^
(1} е Тк).
Результатом отображения является дискретная модель П, с близкими к равносторонним ТКЭ, за исключением элементов, содержащих вырожденные узловые точки, определяемые особенностями границы ^еП.
В отличие от алгоритма [6], в котором конформное отображение осуществляется с использованием интерполяционных полиномов Лагранжа по методу А. Г. Угодчикова [7], построение отображения ю, (р,q) в П, (х,у) произведем решением краевой задачи Дирихле с использованием Лапласиана:
э2х(р „) + =Дх=0 в П, (1)
др2 дq
х( Р, q) = £ (^) на у.
Искомые функции х(р, q) и у(р, q) преобразования координат pq
в ху являются гармоническими и однозначными. Это означает, что всякой функции х(р, q) можно подобрать ей сопряженную у(р, q), которую на основании условий Коши - Римана
дх = ду дх = ду
др д^ дq др
можно определить как
дх дх
ау =-----ар +----aq. (2)
дq др
Решение задачи Дирихле (1) получим методом конечных элементов в форме Галеркина:
Ц агР( х)а ю = 0, (3)
ю
т
где а, - пробная функция; ю = ^ ю, ; т - число подобластей П,, и тогда
,=1
д 2 х д 2 х ^ др2 +дq2
ю
=я!
-а,
да, дх др 1 ^ др) др др
+
-а,
да, дх дq \ дq) дq дq
= 0.
Воспользовавшись следствием теоремы Остроградского - Гаусса, по-
лучим:
У
дх
------(
др
)а у-Я|
ю
да, дх + да, дх др др дq дq
\
= 0.
(4)
интеграл по границе интеграл по внутренней области
где ю = и(У, ю); — = 0, поскольку координаты узловых точек, принадлежа-
ар
щих границе У , заранее заданы и, следовательно, не зависят от координат р и q . Тогда
л=II
да, дх + да, дх
ю ^ др др дq дq
ю
Л
да, дх + да, дх }=1ю { др др дq дq
ГП
=III
л
интеграл по внутренней области Е
да, дх + да, дх
е=1 и { др др дq ^
=III
> Е
dpdq = I Зех = 0.
е=1
(5)
Переход от интегральной формы (5) к системе линейных алгебраических уравнений (СЛАУ) рассмотрим на примере е-го ТКЭ.
Аппроксимацию функций х и у сужением на ТКЭ Тк представим линейными полиномами, используя аффинное отображение ху в ^п (рис. 1)
х = а + + азП; у = а + + азП,
■} =
координаты узлов Тк ;
Т Т
где 1а} = [5]1х}; 1а} = [5]|у}; Iх} ={xl, x2, хз} и {у} = {у1 , У2, уз}- суть
] =
1 0 0
-1 1 0
-1 0 1
(6)
Тогда
3 3
х(п) = Iакхк; у(п) = Iакук, к=1 к=1
где ак = % + Ь2к£ + Ь3кп .
Аналогичные выражения (рис. 1) получатся прямым отображением
pq в
|р(п) = р1 + (р2 -р )^ + ((3 -р1 )n,
1 q ( п) = ql +(2- ql Н+(з- ql )п;
(7)
ю
и обратным в РЧ
^(Р,Ч) = А(Рр3 -Чі) + ЧР -Р3) + РзЧі -РіЧ3) п^Ч) = А((Чі - Ч2 ) + Ч{Р2 - Р1) + РіЧ2 - Р2Ч1)
(8)
Р2 - Рі Р3 - Р1 Ч2 - Чі Ч3 - Чі
где А =
(j = і, 2, 3) - координаты узловых точек на й =
^2^Х2'У2)
удвоенная площадь треугольника; Рj, Чj
й, .
і=і
Рис. і. Аффинное преобразование координат ху и рч в С учетом выражений (6)-(8) производные
дх(£, п). дх(%, п). Эу(%, п). ду(%, п). да,(^,п).да,(^,п)
Эр дч др дч входящие в уравнение (5), определятся как
др ’ дч
дх др дх = [Е ] 3 ^ Ь2кхк к=і 3 > . < ? дУ др ду = [Е ] 3 ^ Ь2кУк к=і 3 » . < ? д«,' др да,
дч ^ Ь3кхк к=і дч ^ Ь3кУк к=і дч
= [Е ]
[Ь2і
Зі,
4і
1 (3 - 41) (43 -41) е11 е12
А _( - Р2) ((2 - Р)_ _е21 е22 _
Опуская здесь и далее знаки суммирования и считая, что в пределах ТКЭ к = 1,2,3; / = 1,2,3, получим
-(е11Ь2к + е12Ь3к )хк,
дх др дх
= (е21Ь2к + е22Ь3к ) хк;
-Р - ( е11Ь2к + е12Ь3к ) хк,
- ( е21Ь2к + е22Ь3к ) хк;
да
др
да
“ - е11Ь2і + е12Ь3і ,
(9)
- е21Ь2і + е22Ь3і •
д^
Подставим в уравнение (5) найденные производные (9) и получим:
Л - Ц[(е21Ь2і + е22Ь3і ){е21Ь2к + е22Ь3к )хк +
£,Л
(е11Ь2і + е12Ь3і )(е11Ь2к + е12Ь3к )хк ]Лй^Л,
где [Л] -
+ (е11Ь2г
др д4
др д4
дп дп
— А - детерминант матрицы Якоби, а учитывая, что
г 1
в пределах ТКЭ на плоскости ^п (см. рис. 1) I йЪ,йп- —, окончательно
% 2 £,п
получим
ЛХ -"2 [(е21Ь2і + е22Ь3і )е21Ь2к + е22Ь3к ) +
+ (е11Ь2і + е12Ь3і )е11Ь2к + е12Ь3к )]хк •
Представим полученные соотношения в виде
[ке]{хк}-{ек} , (10)
где кі,к - (е21Ь2і + е22Ь3і )е21Ь2к + е22Ь3к ) + (е11Ь2і + е12Ь3і )(е11Ь2к + е12Ь3к ) •
т
Построение глобальной матрицы жесткости [К]- ^ [ке - сборка,
е-1
производится известными в приложениях метода конечных элементов приемами, например поузловая, и не представляет никаких вычислительных трудностей. Тогда из решения СЛАУ вида [К]{X}- 0 находим искомый
вектор {X }Т ={х1,..., Xj,..., хп} , где п - число узловых точек регулярной треугольной КЭС на О( х, у).
• ••, , •••,
Несогласованность узлов по склеиваемым границам фг- подобластей
О, учитывается применением алгоритмов конденсации неизвестных [8] в методе подконструкций [9] с последующей поузловой сшивкой [10].
Координаты у^ узловых точек ТКЭ найдем из условия (2), а с учетом
отображения (7) и производных (9) получим
Проинтегрируем последнее выражение в локальных координатах одного ТКЭ, считая, что координаты У1 и У2 (см^ рис 1) известны (обход определим от граничных узлов с известными координатами внутрь области) и будем иметь
у3 - у2 +((е11Ь2к + е12Ь3к )(2 -41 )-(е21Ь2к + е22Ь3к )((2 - р1 ))хк , к - 1,2,3
Осуществляя обход всех ТКЭ, находим искомый вектор координат
Очевидно, что построение эквивалентного аналога ю (і - 3,4,5,6) на плоскости р4 - базовый равносторонний треугольник (БРТ) или область из него порождаемая - выпуклый п -угольник (п - 4,5,6), опирающийся углами на стороны БРТ (рис 2), не является однозначным^
Таким образом, встает задача построения эквивалентного заданной подобласти 0.і вписанного в БРТ п -угольника юг-, наилучшего в смысле нестрогого равенства соотношения сторон соответственно на ху и р4 •
Если для БРТ задача имеет единственное решение в смысле распределения узловых точек по сторонам, то для п -угольников потребуем соблюдения нестрогого равенства соотношения соответствующих сторон
где 1к - длина к-й стороны подобласти О,; 1к - длины к-й стороны п -угольника ю,; £ - наперед заданная положительная малая величина, определяющая качество соответствия О, с ю,.
Вписанный в БРТ п -угольник ю, получим следующим образом.
Шаг 1. Стороны БРТ разобьем на равные отрезки в заданном отношении, обозначим через N число разбиений по стороне треугольника. Тогда длина стороны БРТ Ь = NА!, где для определенности шаг разбиений по стороне А! = 1.
Шаг 2. Построение п -угольников ю, будем проводить так, чтобы их большая сторона совпадала с базовой стороной БРТ (см. рис. 2).
йУ - ~(е21Ь2к + е22Ь3к )хк [( (2 - р )£ + (3 - р )П] + + (е11Ь2к + е12Ь3к )хк [(42 - 41)х£ + (43 - 41))п]•
п -1 п і . і.
Рис. 2. Построение эквивалентных п -угольников ю, на плоскости pq : а - базовый равносторонний треугольник; б - четырехугольник (Тип 1); в - четырехугольник (Тип 2); г - пятиугольник (Тип 1); д - пятиугольник (Тип 2); е - шестиугольник
Шаг 3^ Из очевидных геометрических и пропорциональных соотношений находим длины сторон п -угольников Юі на р4, как показано в табл^ 1
Таблица 1
Типы п -угольников Длины сторон п -угольника юг-
Треугольник - - 3
Четырехугольник - тип 1 (параллелограмм) _ 3М- - - 1, , і 1, 2; Ї3 іі , Ї4 і2 1 11 + 212 + 213 +14 3 4 2
Четырехугольник - тип 2 (трапеция) _ 3М- - - 1 - г , і -1,2; 13 - N 12,14 - 12 11 +12 + 213 + 14 3 2 4 2
Пятиугольник - тип 1 % - Щ , і -1,3; І1 + 2І2 + І3 + 2І4 +15 І2 - N - /[, 14 - N -12 -13, "15 - N -ї4
Пятиугольник - тип 2 ? - Щ , і -1,2; І1 + 2І2 + 2/3 +14 + 2/5 /3 - N -12,15 - N -11 -12,14 - N -13 -15
Шестиугольник - 3М- 1 і , і -1,2,3; 11 + 212 + 213 +14 + 215 + 21б 14 - N -12 -13,16 - N -11 - 12,15 - N -14 -16
^ ^ li j
Шаг 4. Определим относительную погрешность еs = V V-------------------— .
j=1 k=j+1 lk lk
Если найденная погрешность не удовлетворяет условию es < е, то следует изменить количество разбиений N = N + AN, где AN - шаг по числу разбиений. В противном случае считаем n -угольник юг- с найденным числом разбиений по сторонам построенным.
Построение дискретного аналога й)г- проводится следующим образом:
через полученные точки (pj,qj) еФг-, где j = 1,2,...,3N; Ф, - граница ёг-,
последовательно проводятся прямые линии параллельно сторонам БРТ. Пересечением прямых внутри n -угольника юг- получаем узловые точки ТКЭ дискретного аналога йг- на pq. Необходимые сгущения КЭС достигаются на
0. методом последовательных подразбиений [11] или методом «зон и сторон» [12].
Список литературы
1. Зенкевич, О. С. Метод конечных элементов в технике / О. С. Зенкевич. - М. : Мир, 1975. - 296 с.
2. Оден, Дж. Конечные элементы в нелинейной механике сплошных сред / Дж. Оден. - М. : Мир, 1976. - 486 с.
3. Ворошко, П. П. Полуавтоматическая триангуляция плоской области по заданным углам. Алгоритмы и программы по расчету на прочность и исследование напряженно-деформированного состояния элементов конструкций / П. П. Ворошко, О. Н. Петренко. - Киев : Наукова думка, 1979. - 112 с.
4. Cavendish, J. С. Automatic triangulation of arbitrary planar domains for the finite element method / J. С. Cavendish // Int. J. Num. Meth. Eng. - 1974. - V. 8. -P. 679-696.
5. Sadek, E. A. A Scheme for the automatic generation of triangular finite elements / E. A. Sadek // Int. J. Num. Meth. Eng. - 1980. - V. 16. - P. 1813-1822.
6. Кравченко, А. А. Автоматическая дискретизация двумерной области на конечные элементы с использованием конформно отображающей функции / А. А. Кравченко, Д. Ю. Полянский, М. А. Тарасов // Прикладные проблемы прочности и пластичности. Методы решения задач упругости и пластичности : Всесо-юз. межвуз. сб. - Горький : Изд-во Горьк. ун-та, 1985. - С. 115-119.
7. Угодчиков, А. Г. Решение краевых задач теории упругости на цифровых и аналоговых машинах / А. Г. Угодчиков, М. И. Длугач, А. Е. Степанов. - М. : Высшая школа, 1970. - 246 с.
8. Wilson, E. L. The static condensation algorithm / E. L. Wilson // Int. J. Num. Meth. Eng. - 1974. - V. 8, № 1. - P. 198-203.
9. Кравченко, А. А. К вопросу реализации метода подконструкций / А. А. Кравченко, Д. Ю. Полянский // Прикладные проблемы прочности и пластичности : Всесоюз. межвуз. сб. - Горький : Изд-во Горьк. ун-та, 1983. - С. 145-152.
10. Taylor, R. L. On completeness of shape functions for finite element analysis / R. L. Taylor // Int. J. Num. Meth. Eng. - 1972. - V. 4, № 1. - P. 17-22.
11. Темам, Р. Уравнения Навье-Стокса. Теория и численный анализ / Р. Темам. -М. : Мир, 1981. - 408 с.
12. Sienkiewicz, O. C. An automatic mesh generation scheme for plane and curved surfaces by isoparametric co-ordinates / O. C. Sienkiewicz, D. V. Phillips // Int. J. Num. Meth. Eng. - 1971. - V. 3, № 4. - P. 519-528.
Полянский Дмитрий Юрьевич
кандидат технических наук, доцент, ректор Владимирского политехнического колледжа
E-mail: [email protected]
Polyansky Dmitry Yuryevich Candidate of engineering sciences, associate professor, rector of Vladimir polytechnic college
УДК 539.3 Полянский, Д. Ю.
Триангуляция плоских областей решением методом конечных элементов в форме Галеркина задачи Дирихле / Д. Ю. Полянский // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2012. - № 3 (23). - С. 38-46.
