Об одном варианте МКСЭ для уравнений Навье-Стокса

( One FSEM variant for Navier-Stokes equations
Preprint, Inst. Appl. Math., the Russian Academy of Science)

СОДЕРЖАНИЕ

   Введение……………………………………………….………………........3

1. Постановка задачи……………………...…..………………………………4

2. Слабая постановка задачи.....………………………………………………6                

3. Аппроксимация.…..………………………………………………………...9

4. Метод конечных суперэлементов………..…………………………...…..12

5. Вычислительный эксперимент……….……….………..…………………16

6. Заключение…..………………………….………………………………….17

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

         Рисунки ………………………………………………………....................20

ВВЕДЕНИЕ

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

В стандартном  МКСЭ [1] считается, что внутри ячейки решение является «точным» решением рассматриваемого уравнения в классическом смысле, а на границе ячеек удовлетворяет условиям слабой непрерывности потока.

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

В работах [2,3]  схемы МКСЭ строятся применительно к двумерному уравнению конвекции-диффузии, которое следует рассматривать как фрагмент нелинейной системы уравнений Навье-Стокса. Такой фрагмент может появиться при итерационном решении стационарной задачи, основанном на общей идее расщепления системы уравнений и линеаризации, или при решении нестационарной задачи по какой-либо неявной схеме. Опыт решения уравнения конвекции-диффузии был использован при решении системы уравнений Навье-Стокса.

МКСЭ был предложен в 1976 году [4] для расчета задач, моделирующих нейтронно-физические процессы в ядерных реакторах [5], а затем использовался для решения задач теории упругости [6]. МКСЭ - это проекционно-сеточный метод, использующий идеи МКЭ, но в деталях существенно отличающийся от стандартных конструкций МКЭ [7]. Основная цель МКСЭ - это расчет задач на сетках с большим шагом, большим относительно степени гладкости искомого решения, но в то же время – эффективный учет мелкомасштабных неоднородностей внутри ячейки, играющих важную роль в рассчитываемых процессах.

В данной работе представлены расчеты трех модельных задач: задача о тепловой конвекции (GAMM-test 89), течение в каверне с движущейся верхней стенкой и задача с известным точным решением.     Решение нестационарной задачи осуществляется по неявной схеме, на каждом шаге по времени для разрешения нелинейности делается небольшое число (3-5) итераций типа Ньютона. Используется треугольная не разнесенная сетка, одна и та же для скоростей и давления, узлы сетки расположены в вершинах и на сторонах треугольников. Основное внимание уделяется построению пространственной аппроксимации стационарной линеаризованной системы Навье-Стокса методом конечных суперэлементов. В каждом треугольнике строится векторный суперэлемент: по числу узлов на границе треугольника строятся три векторные трехкомпонентные базисные функции с полюсом в данном узле.

Подробные эксперименты по построению неявной схемы, по решению нелинейной задачи и конструированию векторного суперэлемента будут представлены в следующих публикациях.

1. ПОСТАНОВКА ЗАДАЧИ

Система уравнений Навье-Стокса, описывающая течение вязкой  несжимаемой жидкости в ограниченной области  с границей , имеет вид:

,        ;     ,                  (1)

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

функция источника  и     -  заданные функции, а

.

Нестационарную задачу (1) будем решать по полностью неявной схеме, тогда на каждом шаге по времени   необходимо решить стационарную нелинейную задачу:

                                       (2)

Для ее решения строится итерационный процесс, в котором на каждой итерации решается линеаризованная краевая задача [8]. Используются итерации типа

.

 Полученную линеаризованную дифференциальную систему уравнений Навье – Стокса на    итерации  запишем в виде: 

                         (3)      

  .

Здесь       -   известные функции   с предыдущей итерации.

Введём три вектора    , опуская индекс :

    ,      ,                       (4)

При условии  уравнения (3) могут быть записаны в «почти дивергентной форме»:

 ,                                                    (5)                      

 где матрица        и правая часть:       .

Для (5) будем использовать также представление

,                                                    (5а)                      

Аппроксимация векторного уравнения (5) на  треугольной неструктурированной сетке, полученной некоторой правильной триангуляцией области , строится методом конечных суперэлементов.

2.СЛАБАЯ ПОСТАНОВКА ЗАДАЧИ

Рассмотрим формальное интегральное тождество, на которое опирается определение обобщенного (слабого) решения уравнения (5): для любой подобласти

    (6)

здесь - произвольная гладкая вектор-функция с носителем  (пока без требования обращения в ноль на ),    - внешняя нормаль к . 

2.1. В МКЭ алгоритм построения приближенного решения уравнения (5) состоит в следующем:

а) В области  вводится треугольная сетка, в локальной нумерации счетные узлы расположены в вершинах треугольной ячейки (схема 1-го порядка), на сторонах треугольника и, для схем выше второго порядка точности, еще  внутри . В узлах нужно вычислить сеточные значения .

б) В выбирается интерполяционный базис , состоящий, например, из полиномов Лагранжа, или Чебышева, или экспоненциальных функций. Внутрь ячейки приближенное решение продолжается интерполированием, пока неизвестных,  сеточных значений

 .                                                   (7)

Заметим, что в случае неоднородной области при недостаточно малом размере ячейки, даже, если известны точные сеточные значения  решения, приближенное решение (7)  является грубым, не учитывающим сложное поведение решения внутри .

в) Определяется билинейная форма

и скалярное произведение

,

слабым решением уравнения (5) называется [9] непрерывная кусочно-дифференцируемая функция  , удовлетворяющая уравнению

,       ,     ,                                (8)

 т.е. функция должна принадлежать тестовому пространству.

Решение (7) подставляется в (8). Матрица жесткости системы алгебраических уравнений относительно неизвестных собирается из функционалов .

2.2. МКСЭ отличается от МКЭ следующим:

а) счетные узлы расположены в вершинах и на сторонах треугольника, внутри счетных узлов нет.

б)  Врассчитываются два интерполяционных векторных базиса  ,     из специальных краевых задач:

,

             ,     ,        (9)              

3N - число узлов на  ,   M - число узлов в .

Внутри ячейки приближенное решение представляется в виде

,                                           (10)

 - искомые сеточные значения,  - известные сеточные значения правой части.

Задачи (9) в решаются «точно» на подробной внутренней сетке.

в) Слабым решением является непрерывная, кусочно-гладкая функция  , которая в точках гладкости (внутри ) является «точным» классическим решением уравнения (5), а на линиях разрыва первых производных (на  ) удовлетворяет условию слабой непрерывности потока [1],[2],[10]:

                                               (11)

Решение (10) подставляется в (11). Матрица жесткости системы алгебраических уравнений относительно неизвестных собирается из функционалов типа  (11)  от базисных функций . Построенная конечно-элементная схема отличается от (8) более точной аппроксимацией решения внутри и, следовательно, на. .

   Уравнение (11) можно трактовать как условие того, что в качестве носителя  в (5) выбирается крестообразная область меры ноль с центром в узле счетной сетки, а в качестве - характеристическая функция области  [11].

2.3.               В данной работе предлагается следующее определение слабого решения задачи (5). Определим форму 

        (12)

и  тестовое пространство функций  с компонентами   и , поэтому  в контурном интеграле  формулы  (6) после умножения векторов останется только,  где   - внешняя нормаль к .

Слабым решением задачи (5) назовем непрерывную  функцию ==, которая на границе ячеек состоит из  интерполяционных полиномов заданной степени (порядок схемы). Внутрь ячейки она продолжается как приближенное решение уравнения (5) в классе полиномов более высокой заданной степени  и удовлетворяет уравнению

,                                                      (13)

,    

:  ,  ;     .

3.     АППРОКСИМАЦИЯ

3.1. Основная сетка. Предполагаем, что с помощью правильной триангуляции в области    можно построить сетку, в общем случае неструктурированную, с треугольными ячейками  , причем  и любые два треугольника либо не пересекаются, либо имеют общее ребро, либо общую вершину.  Узлы сетки расположены в вершинах треугольников (схема первого порядка), а также, в зависимости от порядка схемы N, на каждой стороне треугольника вводится N-1 узел. Всего на границе треугольника вводится 3N счетных узлов (назовем их объединение ), в которых должны быть определены значения сеточной функции , , здесь - локальный номер узла в ячейке. Объединение  по всем ячейкам назовем сеткой .  Полином , интерполирующий след функции  на  по 3N  узлам, будем записывать в виде:

                    ,                                                                (14)

   -  система функций, образующая на полный  интерполяционный       базис порядка N относительно системы узлов ,  и  . Узел , в котором функция  называется ее полюсом.

3.2. Вспомогательная сетка. Для кусочно-гладкого восполнения известных сеточных функций   внутри ячейки будем использовать также сетку , которая совпадает с сеткой  для , а для  сетка  включает еще внутренние  узлы в треугольнике. Число узлов M сетки  равно , из них  узла лежит внутри треугольника. Полином , интерполирующий, например, функцию  по M узлам в ячейке, будем записывать в виде:

              ,                                                                 (15)

  -  система функций, образующая в полный  полиномиальный базис степени N относительно системы узлов , где .  Для   системы базисных функций   и   совпадают. Для   

3.3.  Конструкция скалярных  функций  , , .

Прежде всего, делается преобразование системы координат    в систему координат  так, чтобы треугольник с вершинами:  перешёл в треугольник  с вершинами  :

            (16)

  -  Якобиан  преобразования.

1

                                                          

Рис. 1

      Обратное преобразование имеет вид:

Координаты нормали , например, ко второй стороне треугольника имеют выражения:

.

Исходным для построения функций ,   является степенной предбазис, состоящий из мономов – степеней переменных , которые обозначаем    .  Структуру предбазиса представим в виде «бесконечного» треугольника Паскаля, занумеровав функции  по строкам слева направо и сверху вниз.

Мономы  разделяем на граничные  и внутренние. Функции    интерполяционного базиса порядка N на границе треугольника  строятся из мономов, расположенных по границе  первых N+1 строк треугольника Паскаля, в виде

,  .                                     (17)

Коэффициенты  находятся обычным способом построения интерполяционного базиса из условия:  . Порядок расположения узлов вдоль границы треугольника и номера граничных мономовдля каждого порядка схемы N строго определены. Узлы занумерованы в направлении против часовой стрелки, начиная с узла , список номеров    задаётся таблицей. Например, чтобы построить функции  для схемы третьего порядка, надо выделить в бесконечном треугольнике Паскаля первые четыре строки и взять в (17) все граничные мономы, начиная с 1. Коэффициенты ,  не зависят от номера треугольника , вычисляются для каждого N один раз и хранятся на диске.

Таким же образом строятся функции  полного полиномиального базиса степени N в треугольнике , в которые входят уже все мономы из первых N+1 строк треугольника Паскаля, включая внутренние:

,  .                                     (18)

Нумерация узлов и мономов идет сначала по границе как в (17), затем по внутренним мономам подобно исходной нумерации . Коэффициенты ,  вычисляются для каждого N также один раз и хранятся на диске. Будет использоваться также система линейно-независимых функций, обращающихся в ноль на границе: 

                                                                              (19)

4. МЕТОД КОНЕЧНЫХ СУПЕРЭЛЕМЕНТОВ

Метод конечных суперэлементов является проекционно-сеточным методом. Построение схемы состоит из двух этапов: построение базиса в ячейке (суперэлемент) и вычисление коэффициентов схемы. Основной конструкцией является построение специального базиса в счетной ячейке, размеры которой велики по сравнению со степенью гладкости искомых функций. Базисные функции должны  отражать эту негладкость решения внутри ячейки. В стандартном МКСЭ для этой цели в ячейке вводится подробная сетка, на которой базис рассчитывается каким-либо конечно-разностным методом, либо МКЭ. В данной работе для построения базиса используется проекционный метод. Для задачи (3) впервые предлагается конструировать векторный суперэлемент [8].

4.1 Векторный суперэлемент. В каждой треугольной ячейке , на границе которой расположено 3N счетных узлов (сетка ), построим систему  из 3×3N  базисных  вектор-функций    (они зависят и от )    как решения элементарных краевых задач:

                (20)

Кроме того, определим еще базисные функции   как решения неоднородного уравнения в с нулевыми значениями на :

,             (21)      

где - полный скалярный интерполяционный базис на   и  - полный скалярный интерполяционный базис в  , являются заданными функциями.

Решение краевых задач в треугольной ячейке осуществлялось: методом коллокаций, Галеркина, методом минимальных интегральных невязок, МКЭ.    Исследовалась точность решения в зависимости от порядка схемы N и способа решения.

Векторным суперэлементом  (СЭ) порядка N назовем треугольную ячейку , оснащенную системами базисных вектор-функций :

                                       (22)

4.2. Приближенное решение в ячейке.

Приближенное решение в ячейке ищется в виде

                                      (23)

здесь  - неизвестные пока сеточные значения искомого решения ,  - локальный номер узла на , - номер компоненты;  - известные значения сеточного представления правой части, - локальный номер узла в . Индекс у  вектор-функции  обозначает номер базисной функции с граничными значениями в (20), а у вектор-функции  индекс обозначает номер базисной функции, являющейся решением уравнения с правой частью в (21).

4.3.  Конструкция базисных функций для однородной задачи. 

Элементарная ячейка сетки   рассматривается как отдельная изолированная область, в которой надо решить серию краевых задач (20), отличающихся данными на границе . Базисные трёхкомпонентные вектор-функции  ищутся как слабое решение (20)  в виде полиномиальных функций:

                          (24)                                         

Пока полагаем , но рассматривается вариант с .

Система уравнений для нахождения коэффициентов  ,   в (24) строится в соответствии с уравнением (13) для однородной задачи с тестовыми функциями .  

                                                          (25)

-  для   это система из  уравнений с   неизвестными  с одной и той же матрицей   с элементами 

Для дальнейшего нам будет удобно представить матрицу системы (25),  состоящей из блоков, сгруппировав блок в соответствии с видом уравнения (5а):

              (26)

Таким образом, для построения функций  надо один раз обратить матрицу  и 3×3N раз умножить   на вектор правой части   .

4.4.  Конструкция базисных функций неоднородной задачи. 

Базисные трёхкомпонентные вектор-функции  ищутся как слабое решение неоднородного уравнения    (21)  в с нулевыми значениями на  в виде:

                                              (27)                                        

Система уравнений для нахождения коэффициентов  ,   строится в соответствии с (13) с тестовыми функциями  .

,                                                 (28)

где   .

Это система уравнений с той же матрицей , что и (25).  Таким образом, для построения функций  надо 2M раз умножить матрицу  на вектор правой части  .

4.5 Вычисление коэффициентов разностной схемы.  После того как найдены коэффициенты .  выражения (23) подставляются в (13), в качестве    берутся  . Уравнение (13) записывается в виде:

,                                                             (29)

Для формирования разностных уравнений предварительно в каждой  ячейке вычисляются функционалы:

.

Собирая матрицу жесткости из функционалов при коэффициентах   по шаблону типа «звезда»  или  «ромб», мы получим систему алгебраических уравнений в   на  итерации (3),   здесь  – номер узла глобальной сетки   в .

                                                           (30)

4.6 Вычислительный алгоритм. К уравнениям (30) надо добавить еще одно уравнение, условие нормировки для : . В модельных задачах для решения переопределенной системы (30) используется стандартная программа  DLSBRR  из  MSIMSL. Сделав несколько итераций () решения нелинейной задачи (2), мы получаем решение нестационарной задачи на  шаге по времени.

5. ВЫЧИСЛИТЕЛЬНЫЙ ЭКСПЕРИМЕНТ

5.1. Задача о тепловой конвекции

Результаты расчетов модельной задачи о тепловой конвекции (GAMM-test 89) по схеме МКСЭ 3-го порядка  приведены на рисунках 2,3:  - прямоугольник X´Y со сторонами  4´1, , , где  - число Грасгофа.  При   устанавливается один вихрь, при – два вихря. Число Грасгофа характеризует степень нелинейности задачи [12]. Расчеты проводились на ортогональной сетке с числом треугольников равным (область делится на прямоугольников, которые делятся диагональю из левого нижнего угла в верхний правый на два треугольника). В первом расчете для установления решения требовалось 4 шага по времени, во втором – 60.

1.           ,           .

2.        .

На рисунках представлена также функция тока :  Stream.fun =

.

5.2. Точное решение модельной задачи:

Исследовалась зависимость точности решения от порядка схемы [13]. На рис. 4 приведен расчет модельной стационарной нелинейной задачи с известным точным решением на той же сетке:

.

.

В таблице приводится зависимость числа итераций  от порядка схемы  и числа Рейнольдса  ,  .  Итерации велись до:  ,        .

5.3. Течение в каверне с движущейся верхней стенкой. На рис. 5,6 приведены два расчета течения в каверне с движущейся верхней стенкой для  и . Аналогичные результаты в [14] получены на более подробных сетках модифицированным вариантом RFB (residual-free bubbles), который близок рассматриваемому в данной работе варианту МКСЭ.

На рис.7 приведены базисные функции в  треугольнике. - первая базисная функция в скалярном полиномиальном базисе в форме Лагранжа для схемы 3-го порядка, используемая в МКЭ и три компоненты первой векторной базисной функции , построенные в предлагаемом варианте МКСЭ.

6. ЗАКЛЮЧЕНИЕ

1. Построен алгоритм  решения двумерной нестационарной системы уравнений Навье-Стокса по неявной схеме на треугольной неструктурированной сетке, общей   для скоростей и давления.

2. Для решения на каждом шаге по времени стационарной нелинейной системы уравнений построено три варианта итерационного процесса, в котором на каждой итерации решается  линеаризованная краевая задача.

3. Для построения схемы МКСЭ базис в треугольной ячейке строится либо как решение скалярного уравнения конвекции-диффузии (элемент системы Навье–Стокса), либо  трехкомпонентный (две скорости и давление) векторный базис строится как решение системы Навье–Стокса, что особенно эффективно при использовании  более точной аппроксимации конвективных членов уравнений. Принципиально трудная задача построения векторного базиса для интерполяции решения уравнений Навье-Стокса другими авторами нам не известна.

4. Создан автоматизированный алгоритм построения схем высокого порядка до 10-го включительно, при этом отсутствуют внутренние счетные узлы  в треугольнике  (нет шаблона «треугольник»)

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

6. МКСЭ позволяет использовать сетки с большим  шагом по пространству,  уменьшая объемы вычислений, но в то же время учитывает сложное поведение решения внутри ячейки.

7. Метод хорошо распараллеливается: расчет базиса в ячейке может осуществляться на отдельном процессоре.