Численное исследование метода конечных суперэлементов на примере решения задачи о скважине для уравнения Лапласа
( Numerical investigation of the Finite Superelement Method for the solution of the “well” problem for Laplace equation
Preprint, Inst. Appl. Math., the Russian Academy of Science)

Введение

Метод конечных суперэлементов (МКСЭ) предложен в работах Л.Г. Страховской и Р.П. Федоренко ([1]) и развит затем в последующих работах ([2]-[5]) как метод дискретизации дифференциальных уравнений, учитывающий мелкомасштабные (относительно шага сетки) неоднородности задачи. Существует широкий класс задач, решение которых содержит резкие неоднородности, проявляющиеся на мелких пространственных масштабах не только по отношению к размеру об­ласти, но даже по отношению к приемлемому шагу сетки H. Такие неоднородности задачи могут быть выделены явно, как, например, инородные включения в композитных материалах, стержни в ядерном реакторе и т.п. Численное решение этих задач се­точными методами требует специальных сеток для разрешения особен­ностей. Для этого необходимо использовать либо адаптивные к решению сетки, сгущающиеся в окрестности особенностей, либо достаточно мел­кие сетки с шагом h и огромным количеством точек. Первый вариант требует применения специальных алгоритмов, второй – соответствую­щей памяти ЭВМ. В то же время проявления особенностей зачастую являются локальными, сосредоточенными в мелкомасштабных подобла­стях. Подтверждением тому является характерный вид функций Грина и типичных решений многих задач математической физики ([6]), а также физические эффекты такие, как принцип Сен-Венана в теории упругости ([7]) и другие. Для решения подобных задач и предложен метод конечных суперэлементов.

В соответствии с МКСЭ ([1]-[5]) расчетная область разбивается на некоторое количество подобластей - суперэлементов таким образом, чтобы особенности решения задачи были целиком сосредоточены внутри суперэлементов, а на их границах решение являлось достаточно гладкой функцией. Каждый суперэлемент (СЭ) оснащается системой “базисных” функций, которые являются точными решениями рассматриваемого уравнения в суперэлементе с некоторыми граничными данными. Решение задачи ищется как линейная оболочка построенной таким образом системы “базисных” функций. Это позволяет, с одной стороны, использовать грубую суперэлементную сетку, а с другой, за счет выбора специальных “базисных” функций, учитывать все особенности решения задачи. При этом дифференциальные свойства искомой функции и ее ограничений на ребра суперэлементной сетки существенно различны. Функция, рассматриваемая лишь на ребрах, может быть достаточно гладкой, значения ее в узлах суперэлементной сетки дают весьма полное представление о поведении функции на линиях этой сетки, но не во всей области. Эффективность же рассматриваемого метода напрямую зависит от качества аппроксимации искомой функции на границах суперэлементов ([2]), и при проведении численного исследования МКСЭ встает вопрос об определении зависимостей ошибок МКСЭ от метода интерполяции следов решения на этих границах. Представляет также интерес прямое сравнение эффективности МКСЭ и обычного метода конечных элементов (МКЭ).

В работе представлены результаты численного исследования МКСЭ. Проведен сравнительный анализ различных вариантов МКСЭ на примере решения модельной задачи. Рассмотренные варианты алгоритма основаны на различных возможных способах задания интерполяции следов решения на границе СЭ и задании регулярного разбиения области на СЭ. Осуществлено сравнение МКСЭ и МКЭ. Получены точные количественные данные об ошибках для модельной задачи.

Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (проект РФФИ № 03 - 01 - 00461).

§ 1. Модельная задача о скважине

Сформулируем следующую модельную задачу для уравнения Лапласа в квадратной области с одной скважиной. Требуется определить функцию , являющуюся решением задачи Дирихле для уравнения Лапласа в двумерной области G. Будем считать, что стороны области G имеют длину  и параллельны осям координат, а скважина представляет собой удаленный из области круг малого радиуса , расположенный в точке с координатами ,  (рис. 1). Границу квадрата обозначим , границу скважины – . На границе скважины поставим краевое условие . В качестве точного решения рассмотрим функцию

,                                      (1.1)

где – полярные координаты на плоскости , – радиус-вектор центра скважины.

Эта функция удовлетворяет уравнению Лапласа всюду в области G и принимает значение  на границе скважины. Функция f, определяющая граничное условие на внешней границе квадрата, задается как ограничение функции  на границу области:

   .

Имеем следующую задачу, записанную в системе координат :

;  , ,

,                                                 (1.2) .

   Рис. 1. Область G со скважиной.

§ 2. Определение вариантов МКСЭ

Построение вариантов алгоритма МКСЭ для модельной задачи основано на подходе, рассмотренном ранее в работах ([2]-[5]). Дадим описание метода их получения для данной задачи.

Расчетную область G разобьем на суперэлементы прямыми, параллельными осям координат

,  , если .                          (2.1)

Будем считать, что расстояние между этими прямыми равно ,  Каждый суперэлемент  является квадратом со стороной H, скважина находится в центре одного из суперэлементов. Общее количество суперэлементов ,  (рис. 2). Рассмотрены разбиения области до значения .

   Рис. 2. Разбиение области на суперэлементы.

  Рис. 3. Суперэлемент и суперэлементные узлы на границе одного СЭ.

Рассмотрим теперь некоторый суперэлемент . На стороне суперэлемента  введем некоторое количество счетных точек, которые вместе с его вершинами образуют набор узлов суперэлемента ,  (рис. 3).

На границе СЭ зададим набор граничных базисных функций ,  ,  , где  – граница СЭ . В узлах СЭ они принимают значения:

,   ,  ,                                (2.2)

,   ,                                               (2.3)

и с этих узлов на грани суперэлемента проинтерполированы одним из следующих способов ([8]-[9]):

·        глобальная полиномиальная интерполяция,

·        кусочно-линейная интерполяция,

·        сплайн-интерполяция,

·        тригонометрическая интерполяция,

·        алгебраическая интерполяция с узлами в точках экстремума многочлена Чебышева.

Более точное описание перечисленных здесь способов интерполяции приведено в §3. При их использовании рассмотрено до 64 узлов в одном СЭ.

В соответствии с МКСЭ решение ищется в виде линейной комбинации “базисных” функций , , которые являются точными решениями рассматриваемого уравнения внутри СЭ, а на его границе совпадают с функциями , то есть для СЭ  и произвольного узла  выполнено

   ,                                 (2.4)

   .

Такая процедура построения “базиса” проведена для каждого СЭ.

Определим сеточную функцию  в суперэлементных узлах . Тогда внутрь каждого суперэлемента искомая функция  проинтерполирована с помощью построенного “базиса” этого элемента следующим образом:

,                              (2.5)

где у “базисных” функций  ради простоты опущен индекс, определяющий номер суперэлемента.

Узловые значения , , соответствующие суперэлементным узлам, определены из уравнений, формально совпадающих с условиями метода Бубнова-Галеркина при выборе в качестве базисных и пробных функций функций  ([2]-[5]):

   ,                                             (2.6)

где – скалярное произведение в L₂, или

   ,       (2.7)

                                        ,                

,                                                       (2.8)

где , .

В дальнейшем в работе решены следующие задачи:

·        Исследование зависимости точности численного решения от способа интерполяции решения на границе суперэлемента при неизменном регулярном разбиении области на суперэлементы.

·        Определение зависимости точности от количества суперэлементов при регулярном разбиении и заданном способе интерполяции на границе.

·        Сравнение эффективности МКСЭ с обычным МКЭ.

Точность расчетов определена в нормах C(G), L₁(G), L₂(G), (G), H(), C() как норма разности точного и приближенного решений в соответствующих пространствах.

Здесь

·        C(G) – пространство непрерывных в G функций;

·        L₁(G) – пространство интегрируемых в G функций;

·        L₂(G) – пространство интегрируемых с квадратом в G функций;

·        (G) – пространство Соболева функций, имеющих в G суммируемые с квадратом первые производные;

H() определим таким образом ([2]). Рассмотрим гильбертово пространство , определенное как , где  – внешняя граница  СЭ,  – пространство интегрируемых с квадратом функций на , . Элемент  пространства  представляет собой упорядоченный набор функций . Введем скалярное произведение и норму прямого произведения: , . В пространстве  подпространство Н определяется следующим образом:

.

Норма по узлам суперэлементов в C() задается как  , где m – общее количество суперэлементных узлов.

§ 3. Реализация различных типов интерполяции

Перейдем к решению первой задачи. Она состоит в определении зависимости ошибки метода от способа интерполяции на границе СЭ для модельной задачи (1.2). Приведем необходимые уточнения для расчетной схемы, описанной ранее в §2, и покажем результаты расчетов.

1. Построение граничных базисных функций

Приведем описание способов интерполяции, которые использованы при построении граничных базисных функций . Значения этих функций в узлах СЭ заданы в соответствии с (2.2). Рассмотрим сначала методы интерполяции на одном ребре границы СЭ.

Пусть N – число суперэлементных узлов  на одном ребре СЭ (N = 2, ...),   – координата узла  в некоторой системе отсчета , связанной с ребром.

Глобальная полиномиальная интерполяция. С узлов  на рассматриваемое ребро СЭ значения  проинтерполированы алгебраическим многочленом степени , определяемым многочленом Лагранжа ([9]):

Кусочно-линейная интерполяция. На каждом из отрезков , , , искомая базисная функция проинтерполирована линейно с учетом заданного условия  ([9]).

Сплайн-интерполяция. Интерполяционный кубический сплайн, определяющий базисную функцию , на каждом отрезке , ,  , задан многочленом ([9])

где , , и наклоны  выбраны из условий

, ,

,

Тригонометрическая интерполяция состоит в замене искомой функции  тригонометрическим многочленом ([9])

,

, коэффициенты которого найдены из заданного условия , , 

Рис. 4.   в СЭ                    Рис. 5.   в СЭ                          Рис. 6.   в СЭ

с глоб.полиномиальной      с глоб.полиномиальной            с кусочно-линейной

интерполяцией 2-го             интерполяцией 2-го                  интерполяцией 

порядка.                                порядка.                                     при n = 12.

Для рассмотренных выше способов интерполирования использована сетка с равноотстоящими узлами . То есть  .

Рассмотрим теперь метод алгебраической интерполяции с узлами в точках экстремума многочлена Чебышева. В этом случае использована следующая неравномерная сетка ([8]):

,  , ,           (3.1.1)

Через эти узлы интерполяции построен алгебраический многочлен по обычной схеме, описанной выше. Такой алгебраический интерполяционный многочлен обладает некоторыми “хорошими” свойствами интерполяции, проявляющимися на чебышевской сетке ([8]). Выбор же в качестве узлов нулей многочлена Чебышева приводит, в отличие от данного способа, к функции, не являющейся непрерывной в углах СЭ. Отметим также, что в узлах  многочлен Чебышева достигает своих экстремальных на отрезке  значений, что и оправдывает название этого способа интерполяции.

Для каждого СЭ построение базисных функций  проведено далее по следующему плану, исходя из алгоритма МКСЭ (§2). Функция  построена в ячейке со скважиной очевидным образом в соответствии с определением (2.3). Число прочих базисных функций (БФ) определяется заданным при разбиении области количеством узлов СЭ n. Каждая БФ  сопоставлена с узлом суперэлемента . На ребрах, содержащих данный узел,  – интерполянт по значениям (2.2) в узлах этого ребра, построенный одним из описанных выше способов. На остальных же ребрах БФ  обращается в ноль  (см. рис. 4-6).

2. Результаты расчетов

Рассмотренные граничные базисные функции  задают в соответствии с алгоритмом метода конечных суперэлементов краевые условия в задаче (2.4) для определения “базиса” . При программной реализации метода “базисные” функции рассчитаны приближенно. При их расчете в данном случае использован обычный МКЭ с кусочно-линейными базисными функциями на треугольной сетке. Приведем на рис. 7, 8 некоторые результаты расчета этих задач в суперэлементе со скважиной для значения .

Решение задачи (1.2) найдено в виде линейной комбинации построенной системы “базисных” функций  (2.5). Сравним с точным решением задачи некоторые приближенные решения, полученные при реализации различных вариантов задания интерполяции следов решения на границе СЭ (рис. 9-11). На рисунках приведены линии уровня  при разбиении области на K = 4 СЭ. На представленных рисунках видно значительное ухудшение качества решения при использовании линейной интерполяции, а также хорошие результаты интерполяции через большее число узлов в СЭ.

Рис. 7а.  Функция              Рис. 7б.   с глоб.            Рис. 7в.   со сплайн-

                                                полиномиальной                 интерполяцией на

                                                 интерполяцией 2-го           границе СЭ при n = 16  

                                                 порядка на границе СЭ      

Рис. 8а.   с кусочно-           Рис. 8б.   с тригоно-     Рис. 8в.   с алгебр.

линейной интерполяцией      метрической интерпо-      интерполяцией порядка

на границе СЭ при n = 28       ляцией на границе СЭ     5 в точках экстремума

                                                 с учетом 2-х гармоник       мн-на Чебышева на гр.СЭ

Рис. 9. Точное

                                                                                                решение

    Рис. 10. Решение задачи,                          Рис. 11. Решение задачи,

      полученное с использованием                  полученное с использованием

      линейной интерполяции на                       сплайн-интерполяции на границе

      границе СЭ (n = 4)                                      СЭ при n = 8

3. Данные об ошибках

Используя полученные при реализации различных типов интерполяции на границе СЭ приближенные решения модельной задачи, а также известное точное решение (1.1), далее подсчитаны точные значения ошибок в нормах пространств, ранее перечисленных в §2:

   ,                                        (3.3.1)

  ,                                       (3.3.2)

      ,                              (3.3.3)

  ,              (3.3.4)

  ,                  (3.3.5)

  .                             (3.3.6)

Приведены зависимости относительных ошибок в нормах (3.3.1)-(3.3.6) от числа узлов в суперэлементе для каждого из рассматриваемых способов интерполяции (рис. 12-14). Разбиение области фиксировано, общее количество суперэлементов .

На рисунке 12 показаны зависимости при использовании алгоритма глобальной полиномиальной интерполяции на равномерной сетке и сетке (3.1.1). На рисунке 13 приведены зависимости относительных ошибок от числа узлов в СЭ при кусочно-линейной, сплайн- и тригонометрической интерполяции на границе.

Остановимся подробнее на ошибках численного решения, полученных при фиксировании кусочно-линейного способа интерполяции решения на границе СЭ. Увеличение числа узлов n в СЭ в этом случае уменьшает ошибку и приводит ее (для больших значений n) к величине “неустранимой” ошибки численного решения (см. рис. 13). Наличие “неустранимой” ошибки связано здесь с влиянием ошибок приближенного расчета задач (2.4) для определения “базисных” функций  и не является свойством непосредственно МКСЭ. Отметим зависимость достигнутой ошибки от качества расчета граничных базисных функций . Граничные базисные функции  вычислены в конечно-элементных узлах на границе СЭ. Избавление от несовпадения узлов P_i и конечно-элементных узлов на ребрах СЭ при расчете  уменьшает “неустранимую” ошибку численного решения (см. табл. 1).

Табл. 1. Значения “неустранимых” ошибок.

На рис. 14 представлены зависимости относительной ошибки от числа узлов в СЭ в логарифмическом масштабе при использовании алгоритма полиномиальной, сплайн  и кусочно-линейной интерполяции на ребрах СЭ.

Из анализа данных рисунка 14 имеем для относительной ошибки численного решения:

                                                 (3.3.7)

где n – число узлов в суперэлементе, и рассчитанное k приведено в табл. 2.

Рис. 12.  Зависимость относительной ошибки от числа узлов в СЭ в нормах

(3.3.1)-(3.3.6) при использовании следующих типов интерполяции:

Фиксировано разбиение области K = 4.

Рис. 13.  Зависимость относительной ошибки от числа узлов в СЭ в нормах

(3.3.1)-(3.3.6) при использовании следующих типов интерполяции:

Фиксировано разбиение области K = 4.

Рис. 14.  Зависимость относительной ошибки от числа узлов в СЭ в логарифмическом масштабе для следующих типов интерполяции на границе СЭ:

Фиксировано разбиение области K = 4.

Табл. 2. Данные о значении коэффициента k соотношения (3.3.7).

Сравнение различных способов интерполяции при использовании небольшого числа узлов в суперэлементе показывает заметное преимущество сплайн-интерполяции в минимизации ошибки численного решения (см. рис. 12, 13). Небольшой по сравнению с интерполянтами других типов коэффициент линейной регрессии k, несвойственный сплайн-интерполяции (см. табл. 2), вызван изначальной близостью ошибки, полученной на небольшом количестве узлов, к величине “неустранимой” ошибки.

В табл. 3 приведены минимальные значения относительных ошибок для модельной задачи, полученные при использовании различных способов интерполяции на границе СЭ и фиксированном разбиении области при .

Табл. 3. Максимальная достигнутая точность на фиксированном разбиении К=4.

Итак, для исследуемой задачи получены численные данные об ошибках при реализации вариантов МКСЭ, основанных на различных возможных способах интерполирования решения на границе СЭ. Проанализировав и структурировав полученные из графиков данные (рис. 12-14), можно отметить, что, как и следовало ожидать, при достижении хорошей аппроксимации решения на границах СЭ независимо от способа интерполяции величина ошибки приближается к своему минимальному значению. За счет выбора варианта задания интерполяции следов решения на ребрах суперэлементной сетки получено значительное уменьшение ошибки во всех нормах, а величина достигнутой ошибки при единственном рассмотренном здесь разбиении области на 4 СЭ позволяет судить об эффективности исследуемого метода (табл. 3). В то же время необходимо отметить, что достигнутая здесь точность не претендует на максимальную при использовании МКСЭ для изначально поставленной расчетной задачи (см. далее §4).

§ 4. Зависимость точности от количества суперэлементов

Следующий шаг – определение зависимости ошибки метода от числа СЭ в области G для модельной задачи (1.2). Зафиксируем конкретный тип интерполяции на границе СЭ. Будем изменять количество СЭ во всей области. Число СЭ , , выбрано так, чтобы скважина при разбиении попала в центр ячейки (§2). Оценка эффективности метода, как и в предыдущем разделе,  проведена в нормах (3.3.1 – 3.3.6). Построены графики зависимостей ошибки численного решения от числа СЭ в области. Зависимости абсолютной ошибки численного решения от числа СЭ в области при использовании алгоритмов глобальной полиномиальной интерполяции первого и второго порядков на границе СЭ приведены ниже (рис. 15 – 16).

Результаты согласуются с теоретическими исследованиями ([4]). При увеличении числа СЭ в области и  приближение границы СЭ к скважине сопровождается возрастанием градиента решения на границе СЭ, что может привести к неудовлетворительной аппроксимации следа решения на этой границе с помощью заданной схемы интерполяции. По приведенным ниже зависимостям видно возрастание величины ошибки при увеличении числа СЭ в области и .

В табл. 4 представлена минимальная достигнутая ошибка численного решения модельной задачи по всем типам разбиений и способам интерполяции на границе СЭ.

Возникает следующая задача – сравнение полученных с помощью МКСЭ ошибок с точностью “обычных” методов. Следующий раздел посвящен результатам решения модельной задачи обычным методом конечных элементов и дальнейшему сравнению полученных результатов.

Табл. 4. Минимальная достигнутая ошибка численного решения модельной задачи по всем типам разбиений и способам интерполяции на границе СЭ

Рис. 16. Зависимость абсолютной ошибки численного решения в нормах (3.3.1 – 3.3.6) от числа суперэлементов в области при заданной линейной интерполяции на границе СЭ.

Рис. 15. Зависимость абсолютной ошибки численного решения в нормах (3.3.1 – 3.3.4) от числа суперэлементов в области при заданной глобальной полиномиальной  интерполяции 2 порядка на границе СЭ.

§ 5. Сравнение эффективности с обычным МКЭ

Задача (1.2) решена методом конечных элементов на нескольких мелких КЭ сетках, с характерными размерами (шаг h):  , , , . Использована равномерная треугольная сетка с кусочно-линейными базисными функциями. Определены ошибки в нормах (3.3.1 – 3.3.4). Исходя из этих результатов, далее полученные величины проанализированы в комбинации с ошибками МКСЭ для модельной задачи. Результаты показаны на графиках (рис. 17-19), отображающих зависимость относительных ошибок МКСЭ для модельной задачи от количества узлов на границе одного СЭ при фиксированном разбиении области для К = 36. Горизонтальными линиями на рисунке отмечены величины относительных ошибок МКЭ. Приведены численные значения ошибок (табл. 5).

Табл. 5. Значения ошибок МКЭ и МКСЭ.

Рис. 17, 18.  Сравнение относительных ошибок численного решения задачи (1.2) при использовании МКСЭ на фиксированном разбиении для К = 36 и относительных ошибок, полученных при расчете этой задачи методом КЭ с различным шагом h. 

Рис. 19, 20.  Сравнение относительных ошибок численного решения задачи (1.2) при использовании МКСЭ на фиксированном разбиении для К = 36 и относительных ошибок, полученных при расчете этой задачи методом КЭ с различным шагом h. 

Заключение

В работе представлены результаты численного исследования МКСЭ. Опробованы различные варианты построения расчетной схемы МКСЭ и проведен их сравнительный анализ на примере решения задачи о скважине для уравнения Лапласа. Построены варианты алгоритма МКСЭ для данной задачи, основанные на возможных способах интерполирования решения на границе суперэлемента и задания регулярного разбиения области на суперэлементы. Для модельной задачи получены точные количественные данные об ошибках. Осуществлено сравнение МКСЭ с обычным методом конечных элементов при различных характерных размерах конечно-элементной сетки. Исследуемые варианты показали свою эффективность при правильном выборе способа построения расчетной схемы МКСЭ.

Список литературы

1. Федоренко Р.П. Введение в вычислительную физику // М.: МФТИ, 1994. – 528 с.

2. Галанин М.П., Савенков Е.Б. О связи метода конечных суперэлементов Федоренко и проекционно-сеточных методов // Препринт ИПМ им. М.В.Келдыша РАН, 2001, № 67.

3. Галанин М.П., Савенков Е.Б., Темис Ю.М. Метод конечных суперэлементов для задач теории упругости // Препринт ИПМ им. М.В. Келдыша РАН, 2004, № 38.

4. Галанин М.П., Савенков Е.Б. Совместное использование метода конечных элементов и метода конечных суперэлементов // Препринт ИПМ им. М.В. Келдыша РАН, 2004, № 13.

5. Галанин М.П., Савенков Е.Б. Метод конечных суперэлементов для задачи о скоростном скин-слое // Препринт ИПМ им. М.В. Келдыша РАН, 2004, № 3.

6. Тихонов А.Н., Самарский А.А. Уравнения математической физики // М.: Наука, 1972. – 736 с.

7. Физический энциклопедический словарь // М.: Советская энциклопедия, 1984. – 944 с. С. 675.

8. Рябенький В.С. Введение в вычислительную математику // М.: Физматлит, 1994. – 336 с.

9. Самарский А.А., Гулин А.В. Численные методы // М.: Наука, 1989. – 432 с.