CC BY
36
УДК 621.315.5:681.3.06
УЧЕТ УСЛОВИЯ ПОСТОЯНСТВА НЕИЗВЕСТНОЙ ФУНКЦИИ В КОНЕЧНО-ЭЛЕМЕНТНОМ КОМПЛЕКСЕ ПРОГРАММ ЕЕМРБЕ80ЬУЕЯ М.И. Батаронова, С.А. Кострюков, В.В. Пешков, Г.Е. Шунин
В рамках вариационной формулировки получены дискретные уравнения метода конечных элементов, учитывающие условие постоянства неизвестной функции и интегральное условие на заданной поверхности
Ключевые слова: дифференциальные уравнения в частных производных, граничные условия, метод конечных элементов
В некоторых физических задачах, описываемых дифференциальными уравнениями в частных производных, может ставиться дополнительное условие равенства нулю тангенциальной производной неизвестной функции ф на внутренней или внешней границе Г1
дф/дт = 0, которое равносильно условию
ф = const (1)
на этой границе.
Данное условие отличается от краевого условия Дирихле тем, что требует лишь постоянство граничных значений функции, но не определяет эти значения. Другими словами, на заданной границе функция сохраняет неизменное значение, которое заранее неизвестно и должно определиться в результате решения краевой задачи.
В ряде случаев условие (1) допускает вариационную формулировку [1]. К примеру, возьмем функционал, соответствующий уравнению Лапласа, в виде
F(Ф) = 2 j ^ф)2 dQ + Q ф|Г] , (2)
2 Q 1
где Q - некоторое число. Если его рассматривать на множестве всех гладких функций, принимающих постоянное значение на Г1, то функция ф, доставляющая минимум функционалу (2), будет одновременно удовлетворять в области Q уравнению Лапласа
V> = 0,
(3)
на части границе, отличной от Г1, - естественному граничному условию
дф
дп
= 0,
и интегральному соотношению ■ дф
(4)
(5)
Последнее равенство возникает только в связи с учетом условия (1). Таким образом, реализация ус-
Батаронова Маргарита Игоревна - ВГТУ, аспирант, тел. (4732) 46-42-22
Кострюков Сергей Александрович - ВГТУ, канд. техн. наук, доцент, тел. (4732) 46-42-22
Пешков Вадим Вячеславович - ВГТУ, канд. техн. наук, доцент, тел. (4732) 46-42-22
Шунин Генадий Евгеньевич - ВГТУ, канд. физ.-мат. наук, профессор, тел. (4732) 46-42-22
ловия постоянства функции автоматически влечет за собой выполнение интегрального условия (5).
Отметим, что приведенная формулировка не исключает условия Дирихле, играющего, как и ф = const, роль главного краевого условия.
Из методов численного решения краевых задач, допускающих вариационную формулировку, наиболее эффективным является метод конечных элементов (МКЭ). Рассмотрим особенности применения МКЭ к решению краевой задачи для уравнения Лапласа с учетом дополнительных условий типа постоянства функции и соотношения (5).
Потребуем, чтобы триангуляция была выполнена так, чтобы элементы, примыкающие к Г1, имели в пересечении с Г1 точку или сторону в двумерном случае, точку, либо ребро, либо грань - в трехмерном.
Как известно, стандартная конечно-элементная процедура состоит в вычислении матричных элементов и вектора правых частей для каждого конечного элемента (т.н. локальной системы линейных алгебраических уравнений) и добавлении их в глобальную систему, относящуюся ко всей конечноэлементной сетке.
Сформируем СЛАУ для трех конечных элементов на рисунке, примыкающих к границе Г1, где задано условие ф = const, а также интегральное условие (5). Функционал для этих элементов представляется в виде
F = + F® + F3 + Q ф|г ,
где Fe) - функционал e-го элемента:
F( ■) = 1 Цффj
где
S(e) = j VN(e)VN(e)dQ - элементы матрицы жест-
кости е-го элемента; ^(е) - базисная функция (функции формы) элемента е, соответствующая узлу /.
Условие ф|г = const приводит к тому, что
ф1 = фз = ф5 = ф*.
Теперь функционал запишется в виде
F=1 ж(!У2+1 5'22)ф2+2 5,зз)ф*2+ф 2+ж(ЗУ2+
+5'2з)ф* ф 2+2 s222y+2 ^зз2)ф*2+1 + ^^фч
+5'242)ф2ф4 + 5'з(42)ф* ф4 + 25'33)ф*2 + 2SS4 + 1 ^5(53)ф*2 +
+ ^з(43)ф * ф4 + ^353)ф *2 + ^453)ф * ф4 + Qф * .
Имеем теперь три узловые переменные (ф2, ф4, ф ) вместо первоначальных пяти (ф1, ... , ф5). Применяя условие минимума функции F^ , ф2, ф4), получим
JF.=(s1(11)+s33)+s33)+^3з3)+s533)+23?+2S353))+
- (+s23)+s22) )ф2+(+s33)+s™ )+q=o,
JF=S+s23+s23)) ф*+ ((+s22) ) ф2+s24) ф4=o,
=(s32)+s33)+s43)) ф*+s22y+(s42)+s43)) ф4=o. (6)
дР_
дФ4
Матрица этой системы остается симметричной и положительно определенной.
Без учета условий (1) и (5) (при этом в функционале (2) нет второго слагаемого) глобальная система этих элементов имела бы вид
^14 + £«ф2 + ЯЗУ = 0,
^12)Ф1 + (('' + Ж22 )ф4 + (и + Я<4з4> )ф3 + ^22)Ф4 = 0 ^1<3)ф1 + (Я23) + ) ф2 + (^33) + ^3<2) + ) ф3 +
+(3? + я33) )ф4+жТфз = о,
42)Ф2+( Я?+Я33) )фз+(^42)+Я43) )ф4+^>5 = о,
Яз(53)Фз + Я43)Ф4 + 4з)Ф5 = 0. (7)
Таким образом, новую систему (6), учитывающую данные дополнительные условия, можно получить простой модификацией системы (7), где этого учета нет. Для этого нужно выполнить следующие действия:
а) все узловые переменные (ф,}, соответствующие узлам на Г1, заменить одной переменной ф*,
б) вместо уравнений, номера которых совпадают с номерами указанных узлов, записать одно новое уравнение, полученное сложением левых частей этих уравнений и последующим приравниванием результата числу (-0).
С точки зрения программной реализации МКЭ важно иметь правило ансамблирования элементов, т. е. формирования данной глобальной системы из локальных систем, относящихся к отдельным конечным элементам. Поэтому приведем локальные системы для всех трех элементов, учитывающие дополнительные условия (кстати, их легко получить из стандартных конечно-элементных уравнений, пользуясь приведенным выше алгоритмом).
Локальная система 1-го элемента:
(+ Я33) + 2 ) Ф * + (Я? + Я23)) Ф 2 = 0,
(Я? + )ф * + Я4Ч = 0.
Локальная система 2-го элемента:
яЗз2)ф * +я22)Ф2 + = 0,
я232)ф * +я322)ф2+ж^у = 0, яЗ42)ф * + ^+^У = 0.
Локальная система 3-го элемента:
(Я3з3)+ж®+4яЗзЗ) )ф *+(Я® + ж? )ф4 = 0, (ж33) + ж43) )ф * + Я43)Ф4 = 0.
Прежде чем осуществить ансамблирование, необходимо каждому уравнению локальной системы поставить в соответствие определенный узел, точнее, степень свободы. Это соответствие определяется по неизвестной, которая умножается на диагональный в данной строке элемент матрицы системы. Например, в первом уравнении 1-й системы такой неизвестной является ассоциированная узловая переменная ф*, значит, этому уравнению следует поставить в соответствие именно эту степень свободы. Во втором уравнении диагональный элемент умножается на ф2, следовательно, этому уравнению соответствует второй узел и т. д.
После проведенной указанной нумерации уравнений построение глобальной системы сводится к стандартной процедуре: уравнение, которое соответствует ,-му узлу, добавляется в ,-е уравнение глобальной системы. Что касается ассоциированной переменной ф*, связанной с дополнительными условиями на Г1, то ей целесообразно дать какой-нибудь уникальный номер (например, номер одного из узлов, лежащих на Г1) и пользоваться им при формировании глобальной системы.
Поскольку все дискретные уравнения (6) получены из соотношений, накладывающие ограничения только на разность функции, а не на саму функции, то возникает неоднозначность в определении МКЭ-решения, для исключения которой необходимо зафиксировать одну из узловых величин (например, принять ф2 = 0). Заметим, что необходимость в такой фиксации отпадает, если в задаче есть краевые условия Дирихле.
Исследуем первое уравнение системы (6). Очевидно, это уравнение учитывает интегральное условие (5): там явно фигурирует заданная величина 0. Перенесем эту величину в правую часть:
(ж®+ж33) + яЗ4)+ж®+ж®+2S1<2)+2ж33)) ф *+ +(Я? + ж43)+ж42) ) ф4+(34)+Я3+Я3) ф4 = _б-
Левая часть этого уравнения преобразуется к виду | N^г+ | ^Г+ | N3^Г = 01 (8)
дп
(1-2) (1-3) (2-3)
Действуя аналогично, определим вклады в левую часть интегрального соотношения (5) от второго и третьего элемента соответственно:
В процессе ансамблирования эти вклады будут суммироваться в полный интеграл (5). В левой части уравнения глобальной системы, соответствующего степени свободы ф*, будет сумма 01 + 02 + 03. Правая часть этого уравнения неизменна - это величина (-0). Видно, что помимо основных интегралов (по сторонам 1-3 и 3-5, целиком лежащим на Г1, от дф/дп), вклад в полный интеграл дают посторонние слагаемые - линейные интегралы по сторонам треугольников, имеющим только один узел на линии Г1.
Если считать функцию дф/дп строго непрерывной в каждой точке сетки, то нетрудно видеть, что при суммировании вкладов соседних элементов, примыкающих к линии Г1, посторонние интегралы по смежным сторонам будут взаимно уничтожаться. В итоге остаться могут лишь несколько слагаемых от элементов со сторонами, являющимися частью внешней границы, не связанной с Г1 (из-за отсутствия смежного элемента). В рассматриваемом случае это интегралы вдоль 4-5 и 1-2. Влияние этой ошибки невелико, так как на внешней границе приближенно выполняется условие дф/дп = 0. Поэтому фактически окончательно формировать уравнение глобальной системы, соответствующее ассоциированной степени свободы ф*, будут линейные интегралы по сторонам, целиком приналежащим Г1.
Очень важный момент связан с тем, что предположение о непрерывности нормальной производной верно только для эрмитовых конечных элементов. Для лагранжевых элементов градиент Уф и тем
самым дф/дп терпит скачок при переходе от элемента к элементу, поэтому о полном уничтожении интегралов по смежным сторонам двух соседних элементов говорить нельзя - присутствует некоторая погрешность.
С целью выяснения характера этой погрешности рассмотрим интегралы по линиям 2-3 и 3-4, которые являются внутренними в рассматриваемом фрагменте конечно-элементной сетки, и определим их суммарный вклад с учетом лагранжевой линейной аппроксимации внутри элементов. Он может быть преобразован в виду
+2(Ф* Ф2)(с!-^а42) -^а®) (11)
(здесь а^е) - внутренний угол треугольника е при вершине, совпадающей с узлом т).
Как уже отмечалось, если бы функция ф была достаточной степени гладкости (в частности класса
С1), то эта сумма непременно равнялась нулю. Однако на линейных лагранжевых элементах градиент Ф - функция кусочно-постоянная, и эта сумма, вообще говоря, отлична от нуля. Интересно отметить, что данные интегралы дают величины, не зависящие от линейных размеров конечных элементов. Влияние геометрии элемента идет только через его внутренние углы. Причем эти углы входят в выражения в виде котангенсов, а это означает, что слишком маленькие углы могут спровоцировать неоправданный рост указанных величин.
Несложно указать один простой частный случай, когда эта сумма точно обращается в ноль. Это когда все элементы - равносторонние треугольники, тогда все углы равны 60° и их котангенсы все одинаковы. Еще в одном случае - когда узлы 1, 3, 5, принадлежащие Г1, с одной стороны, и не лежащие на Г1 узлы 2, 4, с другой, образуют параллельные линии - тогда она достаточно мала.
Приведем примеры вычислений для вполне определенных значений координат узлов. Пусть для рассмотренного ранее фрагмента сетки (рис. 1) эти координаты следующие: (0, 1), (л/3/2, 0.5), (0, 0), (\/3/2, -0.5), (0, -1). Это соответствует случаю строго равномерной сетки: все три конечных элемента правильной формы. Выпишем полный вклад в заряд, разделяя квадратными скобками на три составляющие - от основных интегралов по части линии Г1, включающей 1-3 и 3-5, от интегралов по смежным сторонам 2-3, 3-4 и, наконец, по внешним границам 1-2 и 4-5:
[2.3094ф* - 1.1547ф2 - 1.1547ф4] + [0] +
+ [-0.5774ф* + 0.2887ф2 + 0.2887ф4].
То же самое для слегка измененной, но достаточно регулярной сетки, с координатами узлов (0, 1.2), (л/3/2, 0.4), (0, 0), (-4/3/2, -0.6), (0, -1):
[2.5403ф* - 1.38564ф2 - 1.1547ф4] +
+ [-0.1155ф2 + 0.1155ф4] +
+ [-0.6928ф* + 0.4619ф2 + 0.2309ф4].
А теперь для нерегулярной сетки с координатами узлов (0, 1.2), (1, 0.4), (-0.3, 0), (1.3, -0.6), (0.2, -1):
[2.0240ф* - 1.0625ф2- 0.9615ф4] +
+[-0.4567ф* - 0.0782ф2 + 0.5349ф4] +
+ [0.1715ф* + 0.2292ф2 - 0.0577ф4].
В целом можно сделать вывод, что результаты аппроксимации интегрального условия (5) будут зависеть от качества конечно-элементной сетки, а также её регулярности. Что касается качества сетки, то здесь ничего нового нет, поскольку это характерная особенность МКЭ, состоящая в росте ошибки решения при наличии плохих (вытянутых) элементов. Регулярность же сетки обеспечивает более плавное изменение ф, а значит, более точное вычисление градиента. С другой стороны, с измельчением сетки в МКЭ точность вычисления решения задачи должна повышаться, а значит, и точность учета до-
полнительного условия. В этой связи заметим, что уменьшение размеров рассмотренных треугольников означает, что | ф * -ф20, а это ведет к убыванию лишних слагаемых, представленных в (11).
Если вместо линейных рассмотреть квадратичные элементы и провести для них аналогичные выкладки, обнаружим во многом сходные результаты. В частности, формирование дискретного уравнения, отвечающее за условие (5), обеспечивают не только треугольные элементы со сторонами, целиком принадлежащими Гь но и которые имеют на этой границе свой единственный узел. Причем соотношения, учитывающие вклад в виде линейных интегралов по сторонам этих элементов, в основном повторяют (8-10). Если элементы далеки от правильной формы, то взаимного уничтожения интегралов по смежным сторонам элементов здесь тоже не гарантируется. Однако градиент на квадратичных элементах меняется линейно, его скачки при переходе от элемента к элементу меньше, величина нормальной производной дф/дп под знаком интеграла учитывается точнее, поэтому следует ожидать уменьшения поправок, вносимых указанными интегралами, по сравнению с линейными элементами.
Приведенный алгоритм учета дополнительного условия постоянства функции также легко переносится на трехмерные конечные элементы. Действительно, учет данного условия не меняет исходных элементов матрицы жесткости, а приводит к появлению дополнительных строки и столбца, которые, в свою очередь, формируются исключительно за счет перераспределения имеющихся матричных элементов. Таким образом, внешнее формальное сходство локальных систем треугольника и тетраэдра позволяет считать, что способ учета интегрального условия для двумерного симплекса может быть реализован и для трехмерного элемента.
Как и раньше, перестройке подлежат локальные матрицы только тех тетраэдров, которые примыкают к поверхности Гь на которой задано дополнительное условие.
Воронежский государственный технический университет
Представим ключевые моменты алгоритма учета рассматриваемого дополнительного условия в трехмерном случае:
1) каждому отдельному интегральному условию вида (5) ставится в соответствие ассоциированная переменная ф* - в реальности узловая переменная, связанная с одним из узлов на Г!; в глобальной системе уравнений резервируется соответствующее ей уравнение, под исключенные степени свободы место не выделяется;
2) цикл по всем конечным элементам - формирование локальных систем происходит обычным образом - как в стандартном МКЭ;
3) если хотя бы один узел обрабатываемого тетраэдра принадлежит поверхности Гь данная локальная система подвергается перестройке по одной из схем - в зависимости от того, сколько узлов тетраэдра лежит на Гь
4) локальные системы добавляются в глобальную систему дискретных конечно-элементных уравнений по обычным правилам - в соответствии с локальной и глобальной нумерацией узлов, причем уравнения, которые соответствуют ассоциированной степени свободы ф*, добавляются в соответствующее уравнение глобальной системы;
5) на завершающем этапе в правой части уравнения глобальной системы, которое соответствует Ф*, записывается соответствующее значение величины заряда со знаком минус (-0).
Выводы в отношении погрешности аппроксимации соотношения (5) применительно к фрагменту тетраэдральной сетки остаются верными: для тетраэдров правильной формы, а также расположенных регулярно внутри сплошного слоя указанная погрешность близка к нулю.
Литература
1. Кострюков С.А., Пешков В.В., Шунин Г.Е. Конечно-элементная формулировка дополнительных условий в задачах электро- и магнитостатики сверхпроводников // Изв. РАН. Сер. физ. 2004. Т.68. №° 7. С. 1053.
ACCOUNTING OF CONSTANCY CONDITION FOR UNKNOWN FUNCTIONS IN THE FE PROGRAM PACKAGE FEMPDESOLVER M.I. Bataronova, S.A. Kostryukov, V.V. Peshkov, G.E. Shunin
Using variational formulation, the finite element discrete equations accounting the condition of unknown function constancy and the integral condition on a given surface are obtained
Key words: partial differential equations, boundary conditions, finite element method
