CC BY
28
MCS 65X30
СТРУКТУРА ПРИБЛИЖЕННОГО РЕШЕНИЯ В МЕТОДЕ КОНЕЧНЫХ ЭЛЕМЕНТОВ (ПРЯМОУГОЛЬНЫЕ ЭЛЕМЕНТЫ) С ИСПОЛЬЗОВАНИЕМ БИЛИНЕЙНЫХ СПЛАЙНОВ И СПЛАЙН-ИНТЕРЛИНАЦИИ ФУНКЦИИ
О.Н. Литвин, И.В. Нефёдова
Украинская инженерно-педагогическая академия, ул. Университетская, 16, Харьков, 61003, Украина, e-mail: [email protected]
Аннотация. В работе предложен общий метод построения структуры приближенного решения краевой задачи в области Q, согласно которому проводится выбор координат узлов сетки, который сохраняет глобальную непрерывность приближенного решения в области Q, или непрерывность приближенного решения и его частных производных до n — 1 порядка включительно. Метод существенно использует интерполяцию функций, интерлинацию функций двух переменных на системе взаимно перпендикулярных прямых (линиях ректангуляции) и метод построения базисных сплайнов 2-го порядка.
Ключевые слова: метод конечных элементов, структура приближенного решения, ин-терлинация функций, прямоугольные элементы.
Введение. В данное время метод конечных элементов является одним из наиболее широко используемых методов исследования тепловых, электромагнитных и физико-механических нолей. Как отмечалось на международной конференции «Современные проблемы концентрации напряжений», которая проходила в Донецком национальном университете в 1998 году иод руководством акад. О.С. Космодамнанского, современные системы исследования напряженного состояния корпусов самолетов, автомобилей и т.н., требуют решения многих десятков тысяч уравнений. Это связано с тем, что при формировании систем метода конечных элементов (МКЭ) не учитываются известные в теории упругости факты о возможном наличии на поверхности напряженно-деформированного тола точек, в которых возникает концентрация напряжения, т.е. прогиб поверхности в окрестности указанных точек может иметь особенности. На указанной выше конференции, в частности, был сделан вывод о необходимости разработки методов оптимального выбора узлов в МКЭ.
В роботах |1| - |3| предложен метод построения оптимальных схем метода конечных элементов, в которых находятся не только узловые параметры, а также базисные функции и координаты узлов элементов из условия минимума функционала энергии соответствующего решаемой краевой задаче. Одной из нерешенных проблем этого метода на данное время остается проблема оптимального разбиения области интегрирования
на элементы. На практике при оптимальном выборе узлов элементов (прямоугольные элементы) возникает необходимость разбиения отдельно взятого элемента с наибольшим значением функционала на несколько дополнительных элементов. В результате получаем нерегулярное разбиение области на элементы, которое не исследовалось в цитированных выше работах.
В работах Иво Бабушки и его учеников |4|-|5| был разработан метод сгущения сетки в МКЭ в окрестности точек, в которых есть концентрация напряжений. Этот метод хорошо себя зарекомендован, но проблема эффективного выбора узлов в окрестности особых точек остается неисследованной в случае, когда схема МКЭ использует только прямо-угольные элементы. Это утверждение связано с тем, что в цитированных выше работах Иво Бабушки и его сотрудников предлагается переходить к полярной системе координат с центром в особой точке.
В настоящей работе предлагается общий метод построения структуры приближенного решения краевой задачи в области П на базе использования прямоугольных элементов в МКЭ, который сохраняет глобальную непрерывность решения в области. Метод существенным образом использует интер.нинацию функций двух переменных на линиях ректангуляции (системе взаимно перпендикулярных прямых) |6|,
1. Постановка задачи. Предположим, что область П разделена линиями х = хк, (к = 1 ,т), у = у и (I = 1 ,п) на прямоугольные элементы ГЦ.,/ = [я*, Хк+г] х [у/, У1+1] (к = 1, т — 1, / = 1,??. — 1) ив каждом из этих элементов приближенное решение и (х, у) краевой задачи
Ьи(х, у) = /(х, у), (х, у) € П, и(х, у) = 0 , (х, у) € дП ,
т- / \ д ( (1) ди (х,у)\ д ( (2) ди (х,у)\ , . , .
Ьи {х, у) = -— ^ (х, у) дх У) ду У) и & У) .
р(1),р(2) € С1 (П) , q € С (П) , / € Ь2 (П) представляется в виде: и (х, у) = ик>1 (х,у), (х, у) € Пм С П,
ик,г (х,у) = СМ Н™'0 (в) кк2)'0 (*) + Ск+1,1 к[% (в) к^ (*) + (1)
+Ск,1+1 Нк1)+1(в) Нк2)+1 + Ск+1,1+1 Нк+11г+1(в) Нк2+11г+1(^ = 11
= ^ Ск+^,+ Н(к+^1+и (в) ^1+* № = Ык'1 ^ , ^=0 ^=0
где в = —-——, I = —-, функции /г.д1/ ^ (з) , (¿) € С2 [0,1] и имеют евой-
хк+1 — хк уг+1 — уг
ства НкЦ'0 (0) = кк2)'0 (0) = 1, Н^'0 (1) = к^'0 (1) = 0 к«'1 (0) = к^'1 (0) = 0,
° (1) = hfc2)' ° (1) = 1 V (xk,yi) £ Ö.
В ведем обозначение
11 2 , ,\\ 2
Jkl = J J
°°
+qM (s,t) (s,t) - 2fkj (s,t) wM (s,t)] A^A^dsdt,
л(1) Л (2)
где Ayk; = xfc+i - xk, Д; = yi+i - yh
Pk1i (s, t) = p(1) (sAk1) + xk, tA((2) + yi) , pk2)i (s, t) = p(2) (sAk1) + xk, tA((2) + yi) ,
Qk,i (s,t) = q(sA^ + xk, tA(2) + y^ , fk,i (s,t) = f (sAk1) + xk, tA(2) + yi) .
Значения приближенного решения uk;i во всех внутренних узлах сетки находятся путём минимизации функционала
J (U) = У] Jk,i ^ min .
' ^ Ck i
ПМС G M
2. Основные положения предлагаемого метода оптимизации сетки узлов.
Дадим пошаговое описание алгоритма предлагаемого метода.
Шаг 1. Разбиваем область Ö на элементы np,q С Ö прямыми x = xp , y = yq . В данной работе считаем, что Ö является объединением таких прямоугольных элементов.
Шаг 2. При фиксированном разбиении области интегрирования на элементы nki , при известных базисных функциях и узловых параметрах, вычисляем все значения Jk, i и находим среди них максимальное.
Шаг 3. Если выполняется неравенство max | Jpq - Jk>i| < e, k = p, l = q, процесс
(p,q)=(k,i)
разбиения области па элементы прекращаем.
Шаг 4. Если значение Jp,q в элементе npq С Ö является максимальным, т.е. удовлетворяет следующим неравенствам
|Jp,q - Jk,i| >e, k = p,l = q, (3)
где e — заданное число, то считаем элемент np,q базовым и для оптимизации сетки будем разбивать этот элемент па четыре элемента прямыми
_ _ xp + xp+1 _ _ yq + yq+1
x — xp+0,5 — -^- , У — Vq+0,5 — -^- •
Тогда элемент np,q разобьется на четыре элемента (см. рис. 1):
npq = [xp ,xp+°,5] X [yq ,yq+°,5] ; ^ = [xp, xp+°,5] X ^+°,5, yq+1] ; ^ = [xp+°,5,xp+1] X [yq+°,5,yq+1] ; ^ = [xp+°,5, xp+1] X [yq, yq+°,5] ;
р-1.
р-1.
р-
р-1.
'4 + 2 р, Пр-1,9+1 Ш р- ? + 2 р + 1, /> + 0.5,5+1 + 1, ? + 2 /> + 2, 1 +1,5+1 Л! + 2,
р. 1 »г С 5- ПР) .щ ? + 0.5 ' ^;
Щ® 1 пЩ 1 -"-рд
Щ0, р-Ч~ 1 Р, д р+ 0.5.,^ + ? "" 1 р + 1, 5' -1 р + 2,
!? + 2
з1 +1
*
т-1
Рис. 1. Разбиение области О на прямоугольные элементы.
Таким образом, получено новое разбиение области интегрирования на прямоугольные элементы. Приближенное решение, которое соответствует полученному разбиению, имеет новый вид в элементах П^1 ^ , П^ > П^ 5 Пр^ 5 а также в соседних четырех
элементах ^ ПР-1 ^ Пр ПР+1 ч ■> если эти элементы принадлежат области ин-
тегрирования. Во всех других элементах разбиения форма приближенного решения не изменяется.
В элементах Пр^; ПРПР3^ Пприближенные решения будут иметь вид:
Й<1 ) (х, у) = СР,Я к^ 0 (81) к^ 0 (¿1) + Ср+0.5,д й^ (61) к^ (¿1) + (4)
, ^+0.5 кр, д+0.5 (81) (¿1) + (61) (¿1)
(3), 1
,(4)>
ш
(2)
(х, у) — Ср,д+0.5 кР,д+0.5 (б1) кр, д+0.5 (^2) + Ср+0.5,д+0.5 кР+0.5,д+0.5 (б1) кР+0.5,д+0.5 (^2) +
,(4),0
(3), 1
,(4),0
(61) (¿2) + (61) (¿2)
(5)
^Р3) (х> У) — Ср+0.5,д+0.5 ^+05,4+0.5 ЛР+0и5,д+0.5 (¿2) +
,(3),0
(4), 0
+ Ср+1,д+0.5 д+0.5 (^2) (¿2) +
(4), 0
+ Ср+0.5,д+1 Лр+0.5,д+1 (^2) (¿2) + Ы (¿2) , (6)
,(4), 1
(3), 1
,(4),1
^Р4) (х> у) — Ср+0.5, д ^Р+)005,д (^2) Л-Р4^, д (^1) + СР+1,д ^Р+Тд (в2) ^Р+Тд (^1) + (7)
+ Ср+0.5, д+0.5 ЛР+0.5,д+0.5 (¿1) + ы (¿1)
_ Х - ХР _ Х - Хр+0.5 . _ У - Уд , _ У - Уд+0.5
51 — - , «2 — - , 11 — - , Т'2
Хр+0.5 — Хр Хр+1 — Хр+0.5 Уд+0.5 — Уд Уд+1 — Уд+0.5
^ (61) , ^д'" (¿1) , ^ (82) , ^д'" (¿2) € С2 [0, 1]
и имеют место равенства ^Р30 (0) — Л-Р4)'0 (0) — 1, ЛР30 (1) — Лр40 (1) — 0, ^Р3)д 1 (0) — ЛЙ1 (0) — 0 лР3)д 1 (1) — ЛЙ1 (1) — 1 V (хр, Уд) € П.
Для написания структуры приближенного решения в элементах Пр д-1, ПР-1 д-> ПР д+,, ПР+1 д будем использовать интерлинацию функций на четырех сторонах указанных прямоугольников [7]. В элементе ПР-1 д структура приближенного решения будет иметь вид:
4-1,д (Х У) — и (хр-1 , У) ^Р-^!0) (б) + и (хР, У) Л«' 1 (б) + и (ж, Уд) 42-!°д ф + (8)
+ и (х, Уд+1) ЛР--1,д+1 — Ср-1,дЛР--1,д (б) ЛР--1,д — СР-1,д+1 ЛР-) 1, д+1 (б) ^РЛ,д+1 (¿) —
СР ,дЛр,д (б) Лр,д СР ,д+1Лр,д+1 (б) д'-Ы
ВДВ Хр-1 < X < Хр, Уд < У < Уд+Ь
_ х ХР—1 ^ _ У Уд
Хр Хр-1 Уд+1 Уд
и (хр-1,у) — Ср-1 , дЛР- 10 д + Ср-1 , д+1ЛР-) 1,д+1 >
и (х,уд) Ср-1 , дЛр- 1. д (б) + Ср , д(б)
и (х,уд+1) — Cp-l,д+lhР-i,д+l (б) + СР,д+1ЛР>)д'+-1 (б)
СР,д^Р?)д'0 (^1) + СР'д+0.5ЛР4)д'+0.5 (^1) > уд < у < уд+0.5; СР,д+0.5Лр4í)'-°0.5 (^2) + Ср,д+1ЛР4)д+1 (^2) , уд+0.5 < у < уд+1;
■ _ У - Уд . _ У - Уд+0,5 ¿1 — - , ¿2 — - •
и (ХР,У)
ТО
Уд+0,5 - Уд Уд+1 - Уд+0,5
Лемма 1. Если в элементе Пр-1,д структура приближенного решения имеет вид (8),
4-1, д (хр-1,у) — Ср-1,д ЛРЛ,д + Ср-1, д+^рЛ, д+1 > 4-1,) (х,уд) Ср-1,д Лр-1д (б) + Ср, д (б) ,
4-1,) (х,уд+1) — Ср-1, д+1 ЛР-1, д+1 (б) + СР,д+!ЛР,д+1 (б) , „ ( ) — Г Ср,д^Р4)д'0 (^1) + СР,д+0.5Лр4)д'+0.5 (^1) > уд < у < уд+0.5;
^р-1,д (Хр,у)— < (4), 0 (4), 1
к Ср, д+0.5Лр,д+0.5 (^2) + Ср,д+1Лр,
где
■ _ У ~ Уд . _ У ~ Уд+0,5 н — - ) Ъ — -
Уд+0,5 - Уд Уд+1 - Уд+0,5
Аналогичные леммы формулируются и доказываются дня структур приближенных решений в элементах Пр д-1, Пр д+,, ПР+1 д нового разбиения. Для элемента Пр д-1 приближенное решение будет иметь вид:
4,д-1 (Х У) — и (хр, У) ^Рз1)--1(б) + и (хр+1, У) hpp--) 11 д — 1(б) + и (х> Уд-1) ЛР2)д-1 + (9)
+ и (X, Уд) ЛР2 )д 1 (¿) - Ср,д-1ЛР1 )д-01 (6) Лр2)д-01 (¿) - Ср,д ^Р,0 (б) ^д' 1 (¿) -— Ср+1, д-1 ЛрРн) 1, д-1 (б) ЛР+'1,д-1 — СР+1,д^р+Тд (б) ^р+Тд ,
где Хр < х < Жр+Ь уд_1 < у <
Т» _ 'у
р
хр+1 х
У ~ Уд-1 - У9_1
^ (хр> у) — Ср, 9_1кР, Ь + Ср,9 кР2)„ 1
и (х,у9_1) — Ср , )— 1 крэ, )_1
1к(1)Д (б) )_1кр+1)_1(6)
и (х,У))
где
При этом
и (хр+1,у) — Ср+1,)_1к.Р-)'1>)_1 + Ср+1,9кР+1,д >
Ср,9к(Р,(б1) + Ср+0.5,9кP+0.5'<) (б1) > хР < х < хр+0.5 ; Ср+0.5, 9кp-|-0.5') (б2) + Ср+1,9кр+1,о (62) , хр+0.5 < х < хр+1
Т» _ 'у*
р
61
хр+0.5 - хр
62
х - хр+0.5 хр+1 хр+0.5
^^(хр^) — СР,)_1 > ШР,)_1 С^^+ь^) — СР+1,)_1 >
Ч,^! (х,У)|(хр,уд) — Ср,9 > (х,У)|(хр+о.5— СР+0.5,) >
(х, у) |(Хр+1 ) — СР+Т) •
р, ) +1
шр,9+1(х у) — и (хр, у) кР1 )+1(6) + и (хр+1, у) кРР-)'11)+1 (з) + и(х, У)+1) кР2)+1 ф + (10)
+ и (х, у9+2) кр,)) -12 — СР , 9+1 кр ,)+1 (в) кр ,)+1 — СР , 9+2кр ,)+2 (в) кр ,)+2 —
,(1) ,0
(2), 0 р,
р,
(2), 1 р,
Ср+1,9+1 к(Р-|) 1 ,)+1 кР+1, )+1 С3
(2), 0
'p+1')+2кp+i')+2
где хр < х < хр+ь У)+1 < у < У)+2,
<у* _ т>
р
6 —
хр+1 хр
У - Уд+1 У)+2 - У)+1
и (хр,у) — Ср,9+1 кр,)+1 + Ср,9+2кр,))+12
и (х,у9+2) — Ср , 9+2кр;+02
к(1) Д 9+2кр+1 , 9+2
(6)
4
6
£
и (х,Уд+1)
и (хР+1,у) — СР+1, д+1ЛР+'1, д+1 + СР+1,д+2ЛРн)'1,д+2 >
СР,д+1ЛР,)д+1 (б1) + СР+0.5,д+1ЛР+0.5,д+1 (б1) > ХР < Х < Хр+0.5 ;
Ср+0.5,д+1ЛР+0.5,д+1 (б2) + Ср+1, д+1 ^Рн),, д+1 (б2) , Хр+0.5 < Х < Хр+1 ;
где
XX ^^р XX ^^р | 0 5
51 = "-г- ' 52 _
Хр+0.5 Хр Хр+1 Хр+0.5
При этом
д+1 (х,У)|(хр,уд+1) — СР, д+1 > д+1 (х,У)|(хр+о.5,уд+1) — Ср+0.5,д+1 >
4, д+1 (х'у)|(хр+1,уч+1) — Ср+1,д+1 , 4, д+1 (х, у) |(Жр+1,уд+2) — СР+1,д+2 , 4, д+1 (х,У)|(хр,уд+2) — СР,д+2 • Для элемента ПР+1 д приближенное решение будет иметь вид:
4+1, д (х> У) — и (хР+Ъ У) ^Йт) (б) + и (хр+2) У) ^^д (б) + и (Х Уд) ^Зтд + (П)
+ и (х,уд+1) ЛР+'тд+1 (¿) — Ср+1, дЛР+)!, д (б) ЛР+'1, д (¿) — СР+1, д+1 ЛР+1 д+1 (б) ЛР+'тд+1
—С ь(1)' 1 (б) ь(2)'0 т — С Ь(1)' 1 (б) ь(2)'1
Ср+2,дЛр+2,д (6) Лр+2,д (ь) СР+2, д+1 Лр+2, д+1 (6) Лр+2,д+1 (ь) ,
ГДе Хр+1 < X < Хр+2, Уд < у < Уд+Ь
_ х ~ хР+1 ^ _ У ~ Уд
Хр+2 — Хр+1 уд+1 — уд
и (х,уд) — Ср+1,д^^рР —1, д (б) + СР+2,дhРl)2, д (б) ,
и (хр+2, у) — Ср+2,д)2,д + СР+2, д+1 ЛР+2, д+1 ,
и (х,уд+1) — СР+1,д+1ЛР+'1,д+1 (б) + Cp+2,д+lhР—^2, д+1 (б) ,
тт , , . Ср+1,дЛр+Тд (^1) + Cp+l,д+0.5hp—+ '1,д+0.5 (^1) > уд < у < уд+0.5 ;
и (хр+ЪУ) — ^ (4), 0 (4),,
СР+1, д+0.5Лр+1, д+0.5 (^2) + Ср+1, д+1 Лр+1, д+1 (^2) , уд+0.5 < у < уд+1
■ _ У ~ Уд , _ У ~ Уд+0,5
Н — - ) г2 — - •
yq+0,5 - yq yq+1 - Уд+0,5
При этом
WP+1,9 (X,y)|(xp+i,yq) = CP+1,9 , WP+1,9 (x, y) 1 (xp+i,yq+Q.5) = CP+1, 9+0.5 ,
WP+1,9 (X'y)|(xp+i,yq+i) = CP+1,9+1 , WP+1,9 (X'y)|(xp+2,yq+i) = CP+2,9+1 , WP+1,9 (x,y)|(xp+2,yq) = CP+2,9 •
Шаг 5. Находим среди всех элементов нового разбиения в такой элемент Пр/ 9,, в котором выполняются условия
Jp',9' - Jfc,i| > е,V(k, l) G в , k = p',1 = q',
где е — заданное число, и переходим к шагу 2.
Отметим, что последовательность максимальных значений JP9 будет не возрастающей и, кроме того, такой элемент 1 может быть не единственным.
Теорема 1. Если в указанных восьми элементах разбиения структуру приближенного решения u (x, y) записать согласно приведённому выше алгоритму, а во всех других элементах разбиения её оставить неизменной, то независимо от выбора неизвестных параметров C\j (i = l,m—l,j = l,n— 1) , Clhq+0.5, Cp+0.5,9+0.5; Cp+0.5,9; Cp+o.5,g+i; CP+1,9+0.5 и неизвестных базисных функций hp^, hP29'v, hP39'M, hP49'v G C1 [0,1] с соответствующими индексами 0 < ^, v < 1, приближенное решение будет сохранять непрерывность, т.е. u (x,y) G C (П) и u (x,y) G W2 (П).
Замечание. Предложенный метод разбиения на элементы может быть обобщен также на случай силайи-ингерлииации с использованием базисных сплайнов на неравномерной сетке узлов |6|,
Вывод. В работе предложен общий метод построения структуры приближенного решения краевой задачи с использованием силайи-иитерлииации функции двух переменных дня построения адаптивной сетки узлов МКЭ (прямоугольные элементы), сгущающейся в окрестности точек, в которых точное решение имеет особенности.
Литература
1. Литвин О.Н. Оптимальные схемы МКЭ /7 Теоретические и прикладные вопросы дифференциальных уравнений и а;п'ебра. Киев, 1978. С.160-165.
2. Литвин О.Н. Обобщенная нелинейная интерполяция и решение граничных задач /7 До-повщ1 АН УРСР. Сер.А. 1980. №6. С.23-29. (на украинском языке)
3. Литвин О.Н. Оптимальные координатные функции в методе конечных элементов /7 Дифференц. уравнения. 1984. 20, №4. С 677-688.
4. Babuska I. Finite Element Method for domain with corners // Computing. 1970. 6, №3. P.264-273.
5. Zienkiewiez О.С., Gago J.P.R., Ivcllv D.W. The hierarchic concept in finite element analysis /7 Сотр. Struct. 1983. 16. P.53-65.
6. Литвин O.H. Интерлинацня функций и некоторые ее применения / Харьков: Основа, 2002. 544 с. (на украинском языке)
7. Литвин О.Н., Носов К.В. Численная реализация оптимальших) метода конечных элементов для бнгармопической задачи с краевыми условиями второхх) рода /7 Соврем, проблемы конц. напр.: Тр. междунар. науч. конф. Донецк: ДонГУ, «Кассиопея», 1998. С.151-157.
STRUCTURE OF FINITE ELEMENT METHOD (RECTANGULAR ELEMENTS) APPROXIMATIONS ON BASIS OF BILINEAR SPLINES AND SPLINE-INTERLINATION
O.N. Lytvyn , I.V. Nefedova
Ukrainian Engineering Pedagogical Academy, Universitetskaya St., 16, Kharkiv, 61003, Ukraine, e-mail: [email protected]
Abstract. It is proposed the general method of constructing of the approximate solution of boundary value problem in the domain Q. According to this method, the choice of grid nodes is done. It conserves the global continuity of the approximate solution in Q, or the continuity of the approximate solution and its partial derivatives up to order n — 1 inclusively. The method essentially uses the interpolation of functions and also the interlination of functions of two variables in the system of mutually perpendicular lines (lines of rektangulation) together with the construction of basis two-order spiines.
Key words: finite element method, structure approximate solution, interlination of functions, rect angular elements.
