Application of hermitian biquadratic finite element

Full Text

Билинейные конечные элементы на прямоугольниках давно и успешно используются для решения двумерных стационарных и нестационарных задач [1-3]. C помощью аффинных, изопараметрических и других преобразований область их применения расширена до широкого круга двумерных областей, в том числе с криволинейной границей [1; 3-6]. Однако точность аппроксимации этими конечными элементами невысока: второй порядок в L2 -норме и только первый в H1 -норме. Поэтому интенсивно развились конеч ные элементы с базисными функциями-многочленами более высокой степени, обеспечивающими и более высокий порядок аппроксимации. Причем развитие шло в направлении как лагран-жевых, так и эрмитовых элементов. Сопоставление двух типов элементов дает основание утверждать о большей эффективности эрмитовых элементов по сравнению с лагранжевыми элементами ввиду меньшей размерности порождаемых систем дискретных алгебраических уравнений при равных свойствах 80 Математика, механика, информатика аппроксимации [7]. Более того, для некоторых эрмитовых элементов достигнута не только межэлементная непрерывность, но и межэлементная C'-гладкость, включающая непрерывность первых (частных) производных [1; 8; 9]. Поэтому повышение интереса к эрмитовым конечным элементам остается актуальным. Вместе с тем описанная в литературе линейка эрмитовых элементов начинается с бикубических элементов и продолжается только по нечетным степеням. В этой статье мы опишем двумерный эрмитов элемент второй степени на прямоугольнике, начинающий линейку эрмитовых конечных элементов. Описание элемента. Сначала, используя терминологию работ [1; 2], построим референтный элемент (e, Pê, Z|) как тройку, состоящую из ячейки ê, пространства функций Pê и множества степеней свободы £e. В качестве референтной ячейки возьмем единичный квадрат e = [0, 1] х [0, 1] с четырьмя вершинами aj = (1,1), a2 = (1,0), a3 = (0,0), a4 = (0,1) (рис. і). (однозначной разрешимости) пары (P~, Se) достаточно построить базис Лагранжа {фij (x, у) є P- , j = 1,2, і = 1,...,4} на e, удовлетворяющий условию [1] <V i, j(ф k ,l) = si,k 8 j,l, (3) где 8ik - символ Кронекера. Прямая проверка показывает, что базис Лагранжа имеет следующий вид: фи = x)(1 - x + y), (ф1 2 = xy(x -1), CP21 = x(1 - у)(x + У), cp>2,2 = xy(1 - j)), (4) ф3,1 = (1 - x)(1 - i))(1 + x - j)), ф,2 = x(1 - x)(1 -Ä ф4,1 = (1 - x) j)(2 - x - У), ф4,2 = (1 - x)j)(j) -1). Для проверки интерполяционных свойств этого элемента используем обычные обозначения для пространств Соболева. Пусть Z2(Q) - гильбертово пространство функций, измеримых по Лебегу в области Q, со скалярным произведением (u,v)Q = f uvdQ, u,v є L2(Q) J Q и конечной нормой ll0,Q = (Xu)Q , u Є L2(Q). 0 1 jr Рис. 1. Референтная ячейка Мы определим пространство функций Pê как линейную оболочку восьми полиномиальных одночленов: P) = span {і, x, у, x2, xy, y2, x2 y, xy2}, (1) а множество степеней свободы Se складывается из значения функции и одной из частных производных в каждой вершине квадрата: ={')i,1(jp) = p(aX i = 1,...,4, Vi,2(p) = dp(â)/dx, і = 1,3, (2) Vi,2(p) = dp(a )/ду, і = 2,4, p є P) }. Отметим, что каждой вершине соответствуют две степени свободы, но направления производных различны в вершинах с четными и нечетными номерами. Покажем, что этот набор действительно является корректным конечным элементом в смысле монографии [1]. Лемма 1 . Тройка (e, P~, Se) представляет собой конечный элемент. Доказательство. Размерность пространства P) совпадает с количеством элементов множества Se. Поэтому для доказательства унисольвентности Для целого неотрицательного k обозначим через Hk (Q) гильбертово пространство множества функций u є L2(Q), слабые производные которых тоже принадлежат L2(Q) до порядка k включительно. Норма в этом пространстве определяется формулой lu llk,Q S 0<s+r<k ds+ru dx[ dx2 1/2 (5) Введем также полезную полунорму Ik ,Q S s+r=k дs+ru dx1s dx2 1/2 u є Hk (Q). Пусть û - произвольная функция из H (e). По теореме вложения пространств Соболева H 3(e) непрерывно вложено в Cj()) [10], поэтому u є c1()). В итоге мы можем построить интерполянт u)i є P() : 4 2 uI (^ x2 ) = SS Vi, j (г))іфi, j (^ x2 ). i=1 j=1 Теорема 1. Пусть г) є H3(e). Тогда для любого целого да < 3 справедлива оценка lu - u,| ,< c1 \uV„ (6) I 1 \m,e И I3,e v 7 с константой cj, не зависящей от u. Доказательство. Максимальный порядок частных производных в определении множества Se равен 81 Вестник СибГАУ. № 2(54). 2014 единице. А как уже упоминалось, пространство H3(e) вложено в Cl(e). Кроме того, из (1) следует, что Pê з P2(ê), где P2(e) - пространство многочленов суммарной степени не выше двух. Таким образом, выполнены все условия теоремы 3.1.5 в монографии [1], из которой и следует оценка (6). К этому конечному элементу возможно применение аффинных и изопараметрических преобразований для аппроксимации границы области [1]. Поскольку вдоль границы ячейки базисные функции являются квадратичными, то они предоставляют возможность более точной аппроксимации границы, чем билинейные или линейные элементы. Вместе с тем из-за неоднородности степеней свободы для этого элемента полезно использовать еще одно простое преобразование. Для иллюстрации его необходимости рассмотрим разбиение области Q = (0,1) х(0,1) на элементарные квадратные ячейки (рис. 2), проведя два семейства параллельных прямых xi = ih, і = 1, ..., n -1, и yj = jh, j = 1, ..., n -1, с шагом h = 1/ n. x = xi + hx, y = yj + hy. (7) Рассмотрим две соседние элементарные ячейки в разбиении исходной геометрической области (рис. 3). При использовании преобразования вида (7) (рис. 3, a) получается рассогласование степеней свободы в общих узлах соседних элементов. В принципе, можно ввести формулы пересчета производных из одного элемента в другой. Но это усложнит реализацию метода. Поэтому мы введем еще одно простое преобразование: x = xi + hy, y = yj + hx. (8) * (7) (7) (7) (8) б 0 1 x Рис. 2. Разбиение прямоугольника на элементарные ячейки Обычное преобразование референтного элемента на элементарную ячейку выглядит следующим образом: Применяя его в одном из соседних элементов, мы получим совпадение степеней свободы в узлах сетки (рис. 3, б). Рис. 3. Соседние ячейки с одинаковыми и разными преобразованиями: a - одинаковые преобразования; б - разные преобразования Итак, предложенный эрмитов биквадратный конечный элемент имеет 8 степеней свободы в каждой элементарной ячейке, а каждому узлу ( xt, y. ) соответствует комбинация всего двух базисных функций фі . и у{ J, которые строятся следующим образом. Базисная функция ф, j принимает значение 1 в узле ( xi, yj ) и 0 в других узлах сетки так же, как и ее производные дфі] ( xk, yl )/дх = 0 в узлах с четной комбинацией k +1 и дфі]- (xk, yl )/ôy = 0 в узлах с нечетной комбинацией k +1. Базисная функция равна нулю во всех узлах сетки. Но при четной сумме і + j ее производная по x d^ij ( xk, yl )jdx принимает значение 1, если k = і и l = j, и 0 во всех остальных узлах с четной суммой k +1, а ее производная по y дуi. (xk, yl )/dy обращается в нуль во всех узлах с нечетной суммой k +1. А при нечетной сумме і + j ее производная по y ду ij ( xk, yl )fdy принимает значение 1, если k = і и l = j, и 0 во всех остальных узлах с нечетной суммой k +1, а ее производная по x ду,.(xk,yl)/dx обращается в нуль во всех узлах с четной суммой k +1. В итоге базисные функции имеют следующий вид. Для узла с четной комбинацией индексов і + j 82 Математика, механика, информатика (xi+1 - x ) ( у j+1 - y ) ( hy (x - xi ) + hx ( yj+1 - y ))/ hl Щ при ( x, y) є [xi, xi+1 ]х [ y J, yj+1 ], (x-xi-1 )(.yj+1 -y)(xi -x) + hx (yj+1 -y))/) при (x,y) xi]х[yJ,yj+1 ], (xi+1 -x)(y-yj-1 )((x-(i) + hx (y-yj-1 ))/) при (x,y) фі,xi+1 ]х[УJ-l,yJ], (x - xi-1 ) (y - yj-1 ) (hy (xi - x( + hx (y - yj-1 ))hx2hy2 при (x, y) є [xi-1, xi ]х [yj-1, yJ], иначе; у* = ( x - x)( x+1 - x)( yJ+1 - y V hxh (x - xi)(x - xi-1 )(yJ +1 - y)/hxh (x - xi)( xi+1- x)(y - y.-O/hxh, (x - xi )(x - x-1 )(y - yj-1 )/hxh 0 при (xy) є [xi,x+1 ]х[yJ, при ( x, y) є [xi-1, xi ] X [ yJ, yJ+1 ], при (x, y) є [xi, x+1 ]х [yJ-1, yJ], при (x y) є x-1, xi ] х [yj-1, yJ ], а для узла с нечетной комбинацией индексов і + j (x+1 - x)(yJ+1 - y)(hy(x - x) + hx(y - yJ-1)Vhxh (x - x-1)(yJ+1 - y)(hy(x - xt ) + hx(y - yJ-1 ))/hlh ( xi+1- x)( y - yj-1 )(hy( xi- x )+hx( y j+1- y »/ hxh ( x - xi-1)( y - yj-1 )(hy( x - xi) + hx( y.+1 - y »/ hx2 h 0 ,~sh Фі,] = при (x, y) є h, xi+1 ]х [yJ, yJ+1], при ( x, y) є x-1, xt ]х [ yJ, yJ+1 ], при (x, y) Ф, xi+1 ] х [yJ-1, yJ ], при (x, y) є xxi-1, xt ]х [yJ-1, yJ], иначе; у. = ( y - yJ)( x+1 - x )( yJ+1 - y V hxh ( y - yJ)( x - xi-1)( yJ+1 - y V hxh. ( y - yj )( xi+1- x )( y - yj-О/hxh. ( y - yJ)( x - x-1)( y - yJ-1V hxh 0 при (x, y) Ф, xi+1 ]х [ yJ, yJ+1 ], при (x, y) є x-1, x ]х [yJ, yJ+1 ], при (x, y) є x, x+1 ]х [ yJ-1, yJ], при (x, y ) є xxi-1, xi ]х [yJ-1, yJ ], (9) (10) (11) (12) 0 u ( x, y ) = на Г Численный пример. Проиллюстрируем свойства предлагаемого конечного элемента на следующем примере. Пусть Q = (0,1) х (0,1) - квадрат (см. рис. 2) с границей Г. Рассмотрим краевую задачу d ( du \ d f du Л hxJ-іуl"dyJ=f в n' 0 при x = 0, 0 при y = 0, - y sin y при x = 1, - x sin x при y = 1, с правой частью f ( x, y) = 2 ( x + x2 + y + 3xy + y2 ) cos (1 - x - y ) + + (-y + 2x2y + x (-1 + 2y + 2y2 ) sin (1 - x - y) и коэффициентом ц(x, y) = x + y +1. Точным решением этой задачи является функция u(x, y) = xy sin (1 - x - y). Разделим область Q. на элементарные квадраты, проведя два семейства параллельных прямых xi = ih, i = 1, ..., n -1, и y- = jh, j = 1, ..., n -1, с шагом h = 1/ n. Для выяснения порядка точности при уменьшении размера сетки построим систему линейных алгебраических уравнений методом конечных элементов с использованием базисных функций (9)-(12) для n = 10,20,40. Поскольку точное решение априори известно, то разность u - uh между точным и приближенным решением можно выразить в явном виде. Рассмотрим следующие нормы - дискретные аналоги норм в L2 и H1 : Il0,h 1,h Z (u(xi,y-)-uh(-xi,y-))h2 1<i< n-1, 1< j<n-1 Z (u(xi,y.) - (xi, y.)) h2 + 1<i<n-1,1< j<n-1 + Z (d^dx(xi,(-)-duhjdx(x,,y.)j 2h2 + 1<i <n-1, 1< j<n-1 i + j - четные + Z (^y (xi, (j) -duVdy (xi, y.o)2 2h2 1<i <n-1, 1< j<n-1 i + j - нечетные Напомним, что используемые в них значения производных совпадают со степенями свободы и не требуют дополнительных вычислений или аппроксимаций. 83 Вестник СибГАУ. № 2(54). 2014 Точность приближенного решения h 5h = llu - uh II h II Il0,h ah = u - uh| h II Il1,h 52h 1 5h °2h ! ° h log2(52h 1 5h ) !°g2(°2h < °h ) 0,1 1,15 х 10-6 0,00084 14,7 3,5 3,88 1,8 0,05 7,8 х10-8 0,00024 15,4 3,7 3,94 1,9 0,025 5,07 х 10-9 6,48 х 10-5 Отметим, что с теоретической точки зрения для достаточно гладкого решения задачи гарантируются следующие порядки точности. Теорема 2. Пусть u є H3(Q). Тогда справедливы оценки u - ul < c2h2 Hull Q (13) IIi,q 2 11 |I3,Q и u - ul < c3h3 \\u\\ Q (14) ll0,Q 3 11 ll3,Q с константами c2 и c3, независящими от u и h. Доказательство. Оценка (13) получается стандартным образом [1; 2; 6] из теоремы 1 путем масштабирования и применения к совокупности элементарных квадратов. А оценка (14) вытекает из нее на основании приема Нитше [1; 2]. То есть в нашем примере мы должны получить второй порядок сходимости в норме H1 и третий порядок в норме L2. А практически для дискретных аналогов получаем следующее (см. таблицу). Что касается поведения погрешности CTh = j |u - uh І то ее порядок действительно близок к двум. А вот погрешность 8h = u - uh\ ведет себя 0,h гораздо лучше теоретически предсказанного третьего порядка, демонстрируя близость к четвертому порядку. Это объясняется следующим образом. Оценки (13) и (14) справедливы, вообще говоря, на неравномерных сетках. На неравномерной сетке погрешность в дискретной среднеквадратичной норме действительно будет лишь третьего порядка малости. Но что касается равномерной сетки, то получающаяся конечноразностная схема имеет симметричный шаблон и потому не может быть нечетного порядка точности ввиду сокращения нечетных степеней в разложении Тейлора для погрешности аппроксимации. Поэтому после сокращения слагаемых третьего порядка аппроксимации остаются лишь слагаемые четвертого порядка малости, которые и определяют четвертый порядок сходимости для дискретного набора значений. Но при вычислении (недискретной) нормы ||u - uh II Q она оказывается лишь третьего порядка, как и предсказывается теоремой 2. Итак, в статье представлен новый эрмитов биквадратный элемент на прямоугольнике. До сих пор эрмитовым конечным элементом наименьшей степени был бикубический элемент. Поскольку биквадратный элемент проще бикубического, то для решений класса H 3(Q) он оказывается более экономичным. Относительно неожиданным свойством оказался его повышенный порядок точности в дискретной среднеквадратичной норме на равномерной сетке. Вместо третьего получается четвертый порядок точности. Это объясняется симметрией шаблона получающейся конечно-разностной схемы. В итоге, в этой дискретной норме на равномерной сетке порядок точности биквадратного и бикубического элемента совпадают, что делает первый элемент более предпочтительным ввиду меньшего числа степеней свободы и более простой структуры дискретных уравнений.

About the authors

Vladimir Viktorovich Shaydurov

Institute of Computational Modeling of Siberian Branch of Russian Academy of Sciences; Siberian State Aerospace University named after academician M.F. Reshetnev

Email: shaidurov04@mail.ru
Corresponding member of Russian Academy of Sciences, Doctor of Sciences (physics and mathematics), Professor, Director of Institute of Computational Modeling of Siberian Branch of Russian Academy of Sciences; Director of Institute of Space Research and High Technologies of Siberian State Aerospace University named after academician M.F. Reshetnev.

Sergey Vladimirovich Shut

Siberian Federal University

Email: seshoot@mail.ru
postgraduate student

References

Supplementary files