Содержание к диссертации
Введение
Глава 1. Построение неконформного метода конечных элементов для двумерной задачи Стокса с разрывным коэффициентом в эллиптическом операторе 17
1.1. Классическая постановка задачи Стокса 17
1.2. Основные обозначения и определения пространств 18
1.3. Обобщенная постановка задачи Стокса 23
1.4. Существование и единственность обобщенного решения задачи Стокса 26
1.5. Схема метода конечных элементов 35
1.5.1. Триангуляция исходной области 35
1.5.2. Определение конечно-элементных пространств на подобластях 36
1.5.3. Определение мортарного конечно-элементного пространства на интерфейсе между подобластями 37
1.5.4. Свойства конечно-элементных пространств 38
1.5.5. Определение приближенного решения задачи Стокса . .43
Глава 2. Получение оценок скорости сходимости задачи Стокса с разрывным коэффициентом 44
2.1. Аналог второй леммы Стрэнга 44
2.1.1. Существование и единственность приближенного решения задачи Стокса 44
2.1.2. Формулировка и доказательство аналога второй леммы Стрэнга 46
2.2. Оценка скорости сходимости в норме Vft(ftft) и {p-ph) в норме 49
2.2.1. Вспомогательные утверждения 49
2.2.2. Получение оценки скорости сходимости 59
2.3. Оценка нормы (w — w/J в L2(fi/j) 60
Глава 3. Численная реализация неконформного метода конечных элементов для задачи Стокса с разрывным коэффициентом 68
3.1. Постановка дифференциальной задачи 68
3.2. Конечно-элементное представление 69
3.2.1. Функциональный вид 69
3.2.2. Алгебраический вид 71
3.2.3. Конечно-элементная аппроксимация 73
3.2.4. Преобразование полученной системы линейных уравнений 93
3.3. Построение итерационного процесса 114
3.4. Численный эксперимент и анализ результатов 127
Литература
- Основные обозначения и определения пространств
- Определение мортарного конечно-элементного пространства на интерфейсе между подобластями
- Формулировка и доказательство аналога второй леммы Стрэнга
- Конечно-элементное представление
Введение к работе
Настоящая диссертация посвящена разработке и исследованию неконформного метода конечных элементов для решения двумерной стационарной задачи Стокса с кусочно-постоянным коэффициентом в эллиптическом операторе с использованием мортарных склеек на линии его разрыва.
Изучению стационарной задачи Стокса уделено особое внимание математиками как у нас в стране, так и за рубежом. Интерес к ней связан, прежде всего, с возможностью эффективного применения предлагаемых алгоритмов и подходов для нахождения решения нелинейных уравнений Навье-Стокса.
Первоначально задача Стокса, в ее классической постановке, была решена методами теории потенциалов. Независимо друг от друга Ф. Одквист [91] и Л. Лихтенштейн [86], благодаря построению гидродинамических потенциалов, нашли классическое решение задачи.
Позднее, в работах Д. Лерэ [81]-[83] была дана вариационная формулировка для уравнений Стокса (в общем случае нелинейных) при изучении слабых решений задачи. Существование решения в такой интерпретации вытекает из классической проекционной теоремы (см., например, [28]), при этом само решение называется обобщенным, а вариационная формулировка — обобщенной постановкой задачи Стокса.
В работах Л. Катабриги [55], С. Агмона, А. Даглиса и Л. Нирен-берга [29], В.А. Солонникова [22]-[24], а также в [25], [2], [6], [76], [110], [87] и [7] разработаны различные подходы к исследованию регулярности как классического, так и обобщенного решений стационарной задачи.
Отметим, что при определенных условиях на исходные данные классическое решение задачи Стокса может не существовать, в то время как существование обобщенного решения имеет место. В диссертации рассмотрен один из таких случаев — разрывность коэффициента кинематической вязкости в дивергентно-градиентной части уравнения. Неоднородность этого характера возникает, например, при моделировании физического процесса, протекающего в химическом реакторе.
В связи с этим в настоящей работе (см. также работы автора [13], [14]) для задачи Стокса предложена-новая обобщенная постанов-
ка, основанная на разбиении области на подобласти, где коэффициент вязкости постоянен, и условиях согласования решения на интерфейсе между ними. Для того, чтобы найти обобщенное решение в такой вариационной формулировке, потребовалось разработать эффективный подход к дискретизации задачи с последующей ее численной реализацией.
Созданию и исследованию численных методов решения краевых задач посвящено большое количество работ, которые могут быть условно поделены на две части. К первой из них относятся публикации по аппроксимации задач методом конечных разностей, а ко второй — методом конечных элементов (М К Э).
В своей простейшей форме метод конечных элементов есть процесс построения конечномерных пространств, называемых конечно-элементными пространствами. Эти пространства состоят из.кусочно-полиномиальных функций, определенных на разбиении исходной области.
Впервые М К Э был предложен Р. Курантом в [61], но это важное исследование тогда осталось незамечанным. Затем в начале 50-х годов прошлого столетия метод независимо был переоткрыт инженерами. Наиболее ранние ссылки, широко встречающиеся в литературе, относятся к работам Д. Аргириса [31], [32], М. Тернера, Р. Клафа, X. Мартина и Л. Топпа [100]. Название метода было предложено Р. Клафом [60].
Позднее, в 60-ые годы математиками, начиная с работ С.Г. Мих-лина [10], [11], была показана значимость анализа метода Галеркина и метода Ритца при построении подхода нахождения приближенного обобщенного решения краевых задач. Хотя математикам в то время не были известны достижения инженеров, изучавшиеся ими приближенные методы, как показывают работы Р. Варги [101], Ж. Сеа [56] в одномерном случае и Д. Биркгофа, М. Шульца и Р. Варги [45] в многомерном случае, напоминали М К Э.
Более подробно теоретические результаты М К Э изложены в монографиях Ж. Деклу [3], Г. Стрэнга и Г. Фикса [26], Ф. Сьярле [27], Э. Митчелла и Р. Уэйта [9], Г.И. Марчука и В.И. Агошкова [8], в работах И. Бабушки и Э. Азиза [36], П. Равьяра [93], М. Зламала [106]. В [4], [20], [21] описаны применения М К Э к решению задач, возникающих в различных областях физики и техники.
По своей природе методы конечных элементов разделяются на два типа. Первый тип — согласованные (конформные) М К Э (см. [3], [27], [9], [8], [36], [93]) в том смысле, что
сеточное пространство Sh является подпространством исходного пространства S;
используемые в определении дискретной задачи билинейные и линейные формы тождественно соответствуют формам исходной задачи.
Второй тип М К Э включает в себя методы, нарушающие хотя бы одно из условий согласованности. Из практики вычислений можно выделить три основных способа таких нарушений:
Численное интегрирование (численная квадратура, приближенное интегрирование, приближенная квадратура). При этом подходе все еще имеет место вложение Sh С S, но для вычисления коэффициентов результирующей системы линейных алгебраических уравнений используются квадратурные формулы. (Общее введение в предмет численного интегрирования представлено в работах С. Хабера [77], Ф. Дэ-виса и Ф. Рабиновича [64], дальнейшее его изучение проводилось Л. Мэнсфилдом [89], [90], П. Равьяром [93], Ф. Сьярле [58]).
Метод криволинейных конечных элементов. При этом подходе сеточное пространство Sh не содержится в исходном пространстве S. Такое нарушение встречается при аппроксимации краевой задачи, поставленной в области с криволинейной границей. Внутри области аппроксимация проводится с помощью прямолинейных конечных элементов (элементов с плоскими гранями), а вблизи границы — с помощью конечных элементов, имеющих по крайней мере одну криволинейную грань. Криволинейные конечные элементы используются для того, чтобы как можно лучше аппроксимировать границу области. Преимущественно применяются изопараметрические конечные элементы, впервые предложенные Д. Аргирисом и И. Фридом в [33], И. Ерга-тудисом, Б. Айроном и О. Зенкевичем в [66] и получившие развитие в работе [59]. Наряду с ними, рассматриваются и другие, отличные от изопараметрических, криволинейные конечные элементы. В связи с этим, прежде всего, необходимо отметить работы М. Зламала [107] -[109].
Неконформный метод конечных элементов. При этом подходе также не выполняется вложение Sh С 5, такое нарушение включения, например, для задач второго порядка вызвано тем, что конечно - эле-
ментные функции из пространства Sh при переходе через границы конечных элементов терпят разрыв, то есть не являются непрерывными. Метод был разработан и проанализирован Г. Стрэнгом в работах [98], [99].
Что касается задачи Стокса, то при ее конечно-элементной аппроксимации возникает дополнительная трудность — правильный учет условия несжимаемости. Существуют два способа для ее преодоления. Первый подход состоит в том, чтобы использовать такое конечно - элементное пространство, в котором предполагается, что дискретный аналог уравнения несжимаемости выполнен точно. Однако данный процесс, чаще всего, приводит к усложнению элементов. Методы этого подхода были изучены М. Фортэном в [70]. Второй подход, предложенный М. Крузеем и П. Равьяром, основан на аппроксимации условия несжимаемости. В работе [62] этими авторами были построены и изучены как конформные, так и неконформные конечно-элементные пространства. Дальнейшему исследованию методов конечно-элементной аппроксимации задачи Стокса посвящены также публикации [34], [43], [67] - [69], [72], [73], [102]. Одними из первых вопросами получения оценок погрешности аппроксимации, зависящих от гладкости обобщенного решения, занимались Р. Келлог и Д. Осборн [78], Д. Осборн [92] и Р. Темам [28]. Общий обзор по данной тематике изложен в монографии Ф. Брецци и М. Фортэна [53].
Первыми численными методами для решения стационарной задачи Стокса в переменных "скорость-давление"были итерационные методы К. Эрроу, Л. Гурвица и X. Удзавы, предложенные без обоснования на дифференциальном уровне в работе [35]. Несмотря на то, что с момента появления этих подходов прошло уже более сорока лет, они остаются основой для создания новых методов. Вопросы их применения для задач гидродинамики изучали Ж. Сеа, Р. Гловински и Д. Неделек [57], М. Фортэн [70], М. Фортэн, Р. Пейрэ и Р. Темам [71]. Одни из первых экспериментальных исследований проблемы оптимального выбора итерационных параметров р (или р и а ) были проведены Д. Бежи в [38] и М. Фортэном в [70]. Теоретическое решение этой проблемы для одного частного случая дано М. Крузеем в работе [63]. Отметим ряд публикаций [37], [50], [51], [65], [74], [95], [96], [103], [88] и [79], в которых предложены и изучены модернизированные под-
ходы к построению итерационных методов решения системы линейных алгебраических уравнений, возникающей в результате дискретизации исходной задачи, как с использованием переобусловливающих матриц, так и без них.
В настоящей работе обобщенное решение задачи Стокса определяется независимо на подобластях, а на границе между ними производится его согласование с помощью интегральных соотношений (условий слабой непрерывности). Таким образом, возможно проведение независимой дискретизации задачи на подобластях и, более того, при разбиении каждой из подобластей на конечные элементы не требуется следить за тем, чтобы узлы аппроксимации конечно-элементных пространств совпадали на общем интерфейсе, то есть допустимо использование нестыкующихся сеток.
Тематика применения нестыкующихся сеток для приближенного решения краевых задач привлекает внимание многих исследователей с конца 80-х - начала 90-х годов минувшего столетия. Использование таких сеток неразрывно связано с разбиением, или декомпозицией, области и учетом условий согласования на линиях ее раздела. В связи с этим К. Бернарди, И. Мадей и А. Патерой [44] был предложен новый подход - метод мортарных конечных элементов, который значительно расширяет область допустимых аппроксимаций. Во-первых, в разных подобластях возможно применение различных сеток, не связанных друг с другом, что позволяет:
использовать разные формы сеточных ячеек;
строить "скользящие"сетки, когда аппроксимация задачи с движущимися границами осуществляется за счет сдвига, или скольжения вдоль интерфейса, одной подобласти относительно другой, а не за счет трансформации всей сетки;
генерировать сетки в подобластях независимо от генерации сеток в соседних подобластях;
использовать сетки с сильным скачком шага при переходе через границу между подобластями, это необходимо в случае сильных скачков коэффициентов задачи;
строить легко параллелизуемые методы решения возникающих систем линейных алгебраических уравнений.
Во-вторых, в разных подобластях допустимо применение различных типов метода конечных элементов: базисные функции могут зависить
как от формы ячейки сетки, так и от степени полиномов и используемых в конечном элементе степеней свободы.
В-третьих, в разных подобластях возможно использование различных типов аппроксимаций, например, метода конечных элементов и метода конечных объемов.
В-четвертых, допустимо применение различных моделей в разных подобластях с минимальными условиями согласования решения на интерфейсе.
Метод мортарных элементов, основанный на выполнении условий слабой непрерывности функций на границах между соседними подобластями, оказывается неконформным. Важно отметить, что в первой версии метода, кроме этого, требовалось удовлетворение функций условиям сильной непрерывности в точках на концах интерфейса между подобластями. Во второй же версии мортарного метода, которая уже стала классической, Ф. Бен Велгасемом в [39] было предложено отказаться от этого требования. В работе [40] Ф. Бен Велгасемом и И. Мадей было показано, что использование второй версии метода дает ощутимое преимущество при исследовании трехмерных краевых задач эллиптического типа. Следует подчеркнуть, что на дискретном уровне подобные неконформные методы независимо от них изучали П. Ле Таллек [84], П. Ле Таллек и Т. Соши [85].
Существуют два подхода для построения метода мортарных элементов. Первый подход (см., например, [48], [39]) основан на использовании:
множителей Лагранжа в интегральном тождестве;
дополнительного уравнения в обобщенной постановке задачи — условия слабой непрерывности решения на интерфейсе между подобластями.
При построении приближенного решения задачи получается система линейных алгебраических уравнений, матрица которой обладает седловой структурой и имеет аналогичный алгебраический вид, что и матрица системы при дискретизации задачи типа Стокса. В качестве итерационных процедур для решения вышеупомянутой системы уравнений используют неточный алгоритм Удзавы [50], метод сопряженных градиентов в подпространстве, где седловая матрица положительно определена [52], метод минимальных итераций Ланцоша с блочно-диагональным переобусловливателем [80], [1]. Кроме этого, возможно
применение многосеточных методов непосредственно к задачам с сед-ловым оператором [105], [48].
Второй подход основан на отдельном учете условия слабой непрерывности решения на интерфейсе между подобластями. Итерационные методы решения задачи в такой интерпретации рассмотрены в [97], [75].
Отметим, что в случае разрывного множителя задачи второй подход, в отличие от первого, не позволяет качественно учитывать особенности решения на линиях его разрыва.
Методика получения оценок скорости сходимости приближенных по методу мортарных конечных элементов решений к точным решениям краевых задач эллиптического типа в нормах различных пространств разработана Д. Браессом в [49], Б. Волмус в [104], Ф. Бен Белгасемом в [39].
Решению задачи Стокса с постоянным коэффициентом кинематической вязкости на декомпозиционной области с нестыкующимися сетками на интерфейсе между подобластями посвящены публикации Ф. Бен Белгасема [41], [42], в которых применен второй подход, в качестве условия согласования на общей границе подобластей использовано условие слабой непрерывности вектора скоростей. Основная направленность первой статьи состоит в том, чтобы установить inf — sup неравенство (глобальную div - стабилизацию), в случае произвольного разбиения исходной области на непересекающиеся подобласти; во второй статье автором на границах между подобластями предложено использовать (для тетраэдральных сеток) метод с добавлением интерфейсных bubble функций [54].
В настоящей диссертации для двумерной стационарной задачи Стокса рассмотрен случай разрывного (кусочно-постоянного) коэффициента кинематической вязкости в эллиптической части уравнения. Исходная область разбивается на подобласти так, чтобы на каждой из них коэффициент вязкости был постоянен. На декомпозиционной области предложена новая вариационная постановка задачи, учитывающая условия согласования решения на интерфейсе между подобластями. Решение поставленной задачи определяется как обобщенное, что позволяет установить его существование и единственность в специальных пространствах. Вследствие специфики рассматриваемой задачи, а именно, разрыва коэффициента кинематической вязкости,
неконформный метод конечных элементов характеризуется следующими особенностями: схема М К Э строится независимо на подобластях на основе определения обобщенного решения соответствующей задачи; при этом конечно-элементное пространство для вектор-функции скорости, не являясь подпространством исходного пространства, содержит непрерывные на конечных элементах и кусочно-линейные на подобластях (разрывные при переходе через границу соседних конечных элементов) вектор-функции, а конечно-элементное пространство для скалярной функции давления, являясь подпространством исходного пространства, состоит из непрерывных на конечных элементах и кусочно-постоянных на подобластях функций; сетки на подобластях имеют постоянные, но отличные друг от друга шаги, вследствие чего они не стыкуются на интерфейсе между подобластями; на общей границе двух подобластей введено мортарное конечно-элементное пространство для выполнения сеточных аналогов условий согласования решения.
Благодаря использованию такого подхода и выбору специальных пространств, для задачи Стокса с разрывным множителем в эллиптической части уравнения установлены степенные (по шагу сетки) оценки скорости сходимости приближенного по методу мортарных конечных элементов решения к точному решению. Кроме этого, для решения системы линейных алгебраических уравнений М К Э предложен способ преобразования при помощи исключения ряда неизвестных, в результате чего удается существенно уменьшить их число и количество уравнений, а также построить эффективный итерационный метод решения поставленной задачи с переобуславливаем, используя обобщенный метод минимальных невязок.
Перейдем к более подробному описанию полученных результатов. Диссертация состоит из введения и трех глав.
В первой главе рассматривается двумерная стационарная задача Стокса с кусочно-постоянным коэффициентом кинематической вязкости ао(х, у) в дивергентно-градиентной части уравнения. Для данной задачи исходная область Г2 разбивается на подобласти таким образом, чтобы на каждой из них коэффициент вязкости был постоянен. Для такой декомпозиции области предложена новая вариационная постановка задачи, учитывающая следующие условия согласования на линии разрыва множителя: (1) условие слабой непрерывности компо-
нент вектор-функции скорости; (2) равенство потоков скоростей с давлением на функционалах. Изучен вопрос о существовании и единственности обобщенного решения; а также построена схема метода конечных элементов для нахождения приближенного решения поставленной задачи с использованием мортарных склеек на интерфейсе между подобластями.
В пункте 1.1 приведена классическая постановка задачи в прямоугольной области 1.
В пункте 1.2, благодаря определениям стандартных пространств С.Л. Соболева, введены специальные пространства с нормами, основанные на разбиении исходной области О, на подобласти Qk] изучены их свойства.
В следующем пункте с помощью рассмотренных в пункте 1.2 пространств определена обобщенная постановка задачи Стокса, в которой использованы условия согласования решения на интерфейсе Гі2 между подобластями Qk- Кроме этого, на Гі2 введена вспомогательная вектор-функция Л, совпадающая с суммой потока скоростей и давления на функционалах.
В пункте 1.4 установлено, что при выполнении определенных условий рассматриваемой задачи Стокса существует единственное обобщенное решение (теорема 1.1).
В пункте 1.5 на основе обобщенного решения поставленной задачи построена схема метода конечных элементов. С этой целью в подпункте 1.5.1 выполняется триангуляция Th области Q,. Сначала каждая из подобластей \ вертикальными и горизонтальными линиями разбивается на элементарные квадраты со стороной 2Д&, & = 1,2. Затем каждый из них диагоналями делится на четыре треугольника, которые обозначаются символом е и называются конечными элементами. Далее, вводится 2 ' — триангуляция подобласти Qk, а через &h = U є, Єікн — U е обозначаются разбиения П, 0,^ к = 1, 2,
соответственно. В качестве узлов аппроксимации {а~- } для каждой компоненты вектор-функции скорости на Qk> выбираются середины сторон конечных элементов е, совокупность которых обозначается через S . Подмножество узлов, принадлежащих О^иГ^, определяется как S(k\ к = 1, 2. Выбор узлов аппроксимации для скалярной функции давления не принципиален. В случае h — 2 h\ = 4 / подмножества
узлов аппроксимации для компонент вектор-функции скоро-
сти w на Г\2 ГШі и Гі2 П О2 не совпадают, то есть сетки не стыкуются на интерфейсе Гі2-
В подпункте 1.5.2 введены сеточные пространства на Clkh l.VJj (Q.kh) - подпространство кусочно-линейных непрерывных на е конечно-элементных функций из пространства 1^0) 5 в котором г-я базисная функция (<>f) принимает значение единицы в г-м узле аппроксимации (а\ Є S^) и значение нуль во всех остальных узлах. При этом базисная функция $ тождественно равна нулю вне конеч-
(к)
ных элементов, содержащих узел а\ .
2. ХІ (Q.kh) - подпространство кусочно-постоянных конечно - элементных функций из пространства І^О^/ь/і), в котором г-я базисная функция (ipf) принимает значение единицы на г-м конечном элементе и нуль на всех остальных элементах.
В следующем подпункте определено мортарное конечно - элементное пространство на 1. Отрезки Аі, г = 0, ...,iV, образующие разби-
ение Ti2h = U Аг, называются мортарными элементами. В качестве
г'=0
узлов аппроксимации выбирается подмножество узлов S^ на Гі2, и обозначается через Лд(Гі2л) подпространство линейных на Ду, j ~ 1,..., N ~ 1, и постоянных на До, Ajv, непрерывных конечно-элементных функций из пространства L2(Ti2h)-
В подпункте 1.5.4 определены конечно-элементные пространства со специальными нормами: для вектора скоростей w^ и функции давления ph на всей области О,, для вектор-функции склейки решения на интерфейсе Гі2- Изучены их свойства, в частности, доказано утверждение об эквивалентности норм || Ці^ и || ||уЛ на множестве вектор-функций пространства Y/^f^) (лемма 1.9).
В последнем подпункте пункта 1.5 определено приближенное обобщенное решение задачи Стокса на декомпозиции области Г2 по методу конечных элементов.
Во второй главе установлены априорные оценки (w — w^) в нормах пространств V/^O/i) и 1^(^), на множестве вектор-функций w^ пространства Yft(^), так и (р — рь) в норме пространства Xh(flh), на множестве функций рн пространства Xh{lh), удовлетворяющих интегральному условию постановки задачи.
В пункте 2.1 на основе вспомогательных утверждений получен
аналог второй леммы Стрэнга. В подпункте 2.1.1 сформулировано важное утверждение о том, что дискретизация в вариационной постановке удовлетворяет условию Ладыженской-Бабушки-Врецци формы 0^(-, ) на интерфейсе Гіг (утверждение 2.5). Благодаря утверждениям 2.1 -2.5, установлен факт существования единственного приближенного обобщенного решения задачи Стокса на декомпозиционной области по М КЭ (лемма 2.1).
В подпункте 2.1.2 сформулирован и доказан основной факт пункта
— аналог второй леммы Стрэнга (теорема 2.1).
В пункте 2.2 на основе вспомогательных утверждений установлена априорная оценка ошибок вектора скоростей и функции давления в нормах сеточных пространств.
Подпункт 2.2.1 посвящен доказательству вспомогательных утверждений - оценке слагаемых правой части неравенства теоремы 2.1 (леммы 2.2 и 2.4).
В подпункте 2.2.2 получена оценка ошибок (w — w/J в норме пространства Vh(fifc) и (p — ph) в норме пространства Х^Гі^) (теорема 2.2) с помощью доказанных соотношений лемм подпункта 2.2.1.
В пункте 2.3, благодаря использованию методов и подходов, изложенных в предыдущем пункте, установлена априорная оценка для (w — w/j) в норме пространства 1^() (теорема 2.3).
В третьей главе приведено описание численной реализации неконформного метода конечных элементов для двумерной задачи Стокса в прямоугольнике с кусочно-постоянным коэффициентом в дивергентно - градиентной части уравнения.
В пункте 3.1 дана постановка дифференциальной задачи, а в пункте 3.2 рассмотрено ее конечно-элементное представление.
В подпункте 3.2.1 приведен функциональный вид представления, а в следующем подпункте - подробный алгебраический.
В подпункте 3.2.3 описан процесс построения приближенного обобщенного решения задачи Стокса по М К Э, получена система линейных алгебраических уравнений.
В последнем подпункте пункта 3.2 показано, как в результате преобразований системы уравнений подпункта 3.2.3 можно существенно (почти вдвое) уменьшить число неизвестных и уравнений.
Множество узлов S(k\k = 1,2, разбивается на три типа:
- середины вертикальных сторон элементарных квадратов (первый
тип);
середины горизонтальных сторон элементарных квадратов (второй тип);
середины сторон элементарных треугольников, внутренних по отношению к сторонам элементарных квадратов (третий тип). Переменные компонент вектор - функции скорости в узлах S(k' таким же образом делятся на три типа. Поскольку давление постоянно на каждом элементе е, оно определяется в одной внутренней точке треугольника. В качестве такой точки выбирается середина высоты, проведенной из центра элементарного квадрата на его сторону. В результате множество всех треугольников разбивается на два вида.
Из полученной системы исключаются узлы (и соответствующие переменные) третьего типа:
для уравнений, имеющих дивергентную структуру, берется линейная комбинация на конечных элементах двух смежных элементарных квадратов (8 треугольников);
остальные уравнения преобразуются следующим образом: из уравнений, соответствующих узлам третьего типа, выражаются переменные того же типа и подставляются в уравнения, соответствующие соседним узлам первого и второго типов. (Узел является соседним по отношению к некоторому узлу третьего типа, если они оба принадлежат одному и тому же элементу триангуляции Tfr').
Пункт 3.3 посвящен построению и исследованию итерационного процесса для решения преобразованной системы уравнений, матрица которой имеет седловую структуру. Для нахождения решения системы построен итерационный процесс с переобуславливанием матрицы системы, состоящий из четырех этапов.
В ходе численного эксперимента (пункт 3.4) создана система программ на языке Си для нахождения приближенного решения с помощью неконформного метода конечных элементов, проведены расчеты серии модельных задач с различными разрывами коэффициента. Проделан численный анализ полученных результатов и сделаны выводы об аппроксимационных свойствах неконформного метода конечных элементов для решения задачи Стокса данного типа.
Основные положения диссертации опубликованы в работах [94], [12] - [19]. Полученные результаты докладывались на Второй Международной конференции по вычислительной математике (МКВМ-2004;
Новосибирск, 2004 г.), на Дальневосточной школе-семинаре по математическому моделированию и численному анализу (FESS-MMNA'01; Находка, 2001 г.), на Конкурсах-конференциях молодых ученых и аспирантов Хабаровского края (Хабаровск, 2004, 2005 г.г.), на Дальневосточных математических школах-семинарах имени академика Е. В. Золотова (2002-2004 г.г.), на совместных семинарах лаборатории математического моделирования в физике и технике ВЦ ДВО РАН, кафедр математического анализа Хабаровского государственного педагогического университета и систем автоматизированного проектирования Дальневосточного государственного университета путей сообщения а также на семинаре Института математического моделирования РАН (Москва).
В тексте диссертации принята двойная нумерация утверждений и формул: первая цифра указывает на номер главы, а последующие -на порядковый номер утверждения или формулы (при этом для нумерации однотипных формул используются также буквы русского алфавита, например, (3.14а) и (3-14^)).
Основные обозначения и определения пространств
Пусть QQ С R2 - произвольная выпуклая ограниченная область с достаточно гладкой границей TV Через 7о обозначим кусок границы Го, а через По замыкание области По, По = По U Го.
Через С00 (По) обозначим пространство бесконечно дифференцируемых функций на По, а через Со(По) его подпространство, содержащее функции с компактным носителем в П0.
Как обычно, H3(Q,o) - пространство функций из 1/2(По) (суммируемых с квадратом) таких, что их все производные, до s-ro порядка ВКЛЮЧИТеЛЬНО, фуНКЦИИ ИЗ 1/2() (H(QQ) = 1/2(По)), то есть Я5(П0) = {z є L2(n0);Vm, \т\ s,Dmz є 2(П0)}, с нормой IWU = (E {[D dxdyf r (1-2) где дИг т = a a m = ть m2 mi -0,m2 - " целые и m = ГП\ + 77.
Определим подпространство Яо(По) С ІЗ"1 (По) как замыкание по норме (1.2) множества функций из пространства Со(По); при этом норма на Яд (По) совпадает с нормой пространства Я1 (По).
Введем пространства Hs(П0) и Но(По) вектор-функций w = (w, г ), каждая компонента которых принадлежит Я3 (По) и Яд (По) соответственно. Норма в Hs(Ho) определяется как ми = (н;л + н\1л)1/г Через ) \3д0 ( \s,nQ) обозначим полунорму на пространстве Я3(По) (HS(HQ)) И определим ее следующим образом: Ы.Л = ( Е / ID ufdzdy)1 (wsA = («n„ + Мїд,)1/"). lml=J»J2o Далее, введем пространство Я1/2(70) как множество функций г из 1/2(7о); Для кажДй из которых существует такая функция z из пространства Я1 ), что оператор следа этой функции z совпадает с функцией г на 7о: # 1/2Ы = {г є Ыю); 3z Є ЯХ(ОД, k = г}, с нормой г70=г Через Н1/2(7о) обозначим пространство вектор - функций г = (гі, гг), каждая компонента которых принадлежит Я1/2(70), с нормой Irlv2,70 = (1Ы?/2,70 + Ы\1Рл)ф 1/2
Через Я0о (7о) обозначим множество таких функций ц(х, у) из пространства Я1 2(70); которые можно продолжить нулем на Го\7о ДО функций Д(ж, у) Є Я1/2 ): то есть #0 (70) = { є я1/2Ы; Д є #1/2(г0)}, с нормой IH#»(,) = ІІАІІ1АГ. = Д, ІИІіл гІ70=ЛгІг0\70=0 Определим (Я0о (7о)) как дуальное пространство к Н0д (70) (относительно 2(70)) с нормой Ли Ы)У= SUP 7 дєя#ам И Ия ы Обозначим через Н0о (70) пространство вектор - функций М = (мъ А ), каждая компонента которых принадлежит Я0д (70) j с нормой к2ы = (И 11я0 ы + Ы1я ы)1/2 Теперь введем определения и обозначения пространств для области О, из пункта 1.1. Обозначим через wk = (ик, vk) сужение вектор - функции w = (и, v) на подобласть ГІ/г, к = 1,2.
Далее, через [Z]\T12 обозначим функцию разности двух следов на общей границе двух подобластей, то есть если zl\r12 — след функции zl = 21 Є Н1{р,і) т. Г12, a z\12 — след функции z2 = z\n2 Є Hl(fk) на Гі2) то йг12 := zl\Tl2 - z2\Tl2 на Гі2, Гі2 = Пі ГіП2.
Аналогичную запись будем использовать и для вектор - функции z, при этом полагая, что вышеизложенное определение относится к каждой из компонент.
Функцию (вектор - функцию) Игі2([2]гі2) будем называть скачком функции (вектор - функции) z(z) на интерфейсе Гі2 Через и(П) определим множество вектор-функций w = (гл, v) на П, сужение каждой из компонент которых на к-ю подобласть принадлежит пространству Я1(П/с),/с = 1,2, удовлетворяющих условию w = 0 на Г, то есть U(ft) = {w= (г , «);« Є L2(fi), v Є Z fi), и = «flfc Є ЯХ(ВД, vk = v\ak Є ЯХ(ОД, w = 0 на Г, A; = 1,2}, с нормой llwllufl = (\M\U + \H\l,n2)1/2 Через У(П) обозначим множество вектор-функций из пространства U(О), скачок каждой компоненты которых на Гіг принадлежит пространству Я0о (Г12), то есть V(fi) = {w = (г , г;) Є U(ft); МГі2 Є Я010/2(Гі2), МГи Є #І/2(Г12)}, с нормой lwv,« = (HwH + ИМ /3(Гіг)) 2. Определим пространство Y(ft) = {w = (и, v) Є V(ft); / НГі2 - i/dy = 0, Гі2 с нормой пространства V(H).
Определение мортарного конечно-элементного пространства на интерфейсе между подобластями
В этом пункте построим схему МКЭ: зададим сетку на подобластях Q/., k = 1,2, а следовательно, и на О, определим конечно - элементные пространства, сначала на подобластях и интерфейсе между подобластями, затем на всей области, изучим их свойства и, наконец, определим приближенное обобщенное решение рассматриваемой задачи по МКЭ.
Выполним триангуляцию Тд области О. Сначала каждую из подобластей Пі и 2 вертикальными и горизонтальными линиями разбиваем на элементарные квадраты со стороной 2 hi и 2 / соответственно. Затем каждый из них диагоналями делим на четыре треугольника, которые будем обозначать символом е и называть конечными элементами.
Далее, введем Тд и Тд — триангуляции подобластей Пі и П2, а через Пд = U є, Під = U е и Пгд = U е обозначим разбиения еєТ(і) ee2f 0, 0\ и 2 соответственно. Отметим, что построенная триангуляция квазиравномерна на подобластях (см., например, [27]).
В качестве узлов аппроксимации {с - } для каждой компоненты вектор-функции скорости на 0 , выберем середины сторон конечных элементов е, совокупность которых обозначим через S . Подмножество узлов, принадлежащих OkijTu, определим как S k\ к = 1, 2.
Выбор узлов аппроксимации для скалярной функции давления не принципиален, исходя из введенного для нее ниже конечно - элементного пространства.
Рассмотрим случай h = 2 hi = 4/i2, тогда подмножества узлов аппроксимации 5(1) и 5(2) для компонент вектор-функции скорости на Гіг П 0,1 и Гіг П 2 не совпадают, то есть сетки не стыкуются на интерфейсе Гі2 Определение конечно-элементных пространств на подобластях
Введем конечно-элементные пространства на Okh l.Vh (Okh) - подпространство кусочно-линейных непрерывных на е конечно-элементных функций из пространства І П&д), базисные функции в котором определены следующим образом: V /? Є V (Okh) : rf(aT]) = 1 и tf(a ) = 0, у, ,« І,4Ч Є 5»), af є 3 , = 1,2. При этом базисная функция (р% тождественно равна нулю вне конеч ных элементов содержащих узел а) , к — 1,2. 2.Хд (Okh) - подпространство кусочно-постоянных конечно-элементных функций из пространства 1/2 ( ) , базисными функциями которого являются функции ф!? такие, что: Vi/jf Є Х (Okh) (еі) = 1 и $(е,-) - 0, Vj 7= і, є , Є,- Є if, fc - 1,2.
Сеточное решение системы (1.4), (1.5), на подобластях Г2&, /с = 1,2, будем искать в виде и (я?,з/)= Е viffay), {i;ofc)65№)} {i;afc)65( )}
Введем конечно-элементное пространство на интерфейсе Г и-Отрезки До - [(тг, 0), (тг, Лі)], Aj = [(тг, h + 2 (j - 1) /н), (тг, /її + 2 j /її)], j = 1,..., TV - 1, Д = [(тг, + 2( -1)/ ), (тг, 7T)},NGN: N= N образуют разбиение Гш = U Дг- и называются (см., например, [44]) г=0 мортарными конечными элементами. В качестве узлов аппроксимации выберем подмножество узлов S на Гі2 и обозначим его через S: S = (}1 = {fah1 + 2jhi),3 = 0,...,JV -1, iV = -}. Л/1(Гі2/1) — подпространство линейных на Д./, j = 1,..., TV — 1, и постоянных на AQ, AN, непрерывных конечно-элементных функций из пространства 1/2 (Г д), базисные функции которого определены следующим образом:
1) щфо) = 1, vo(bj) = 0 V? = 1,..., iV - 1, Ї/О(ТГ, з/) = 1 при г/ Є Д0, і/о(тг, 2/) при у Ді — линейная,г/0(тг, 2/) = 0 при у Є Д? Vj = 2,..., TV;
2) -() = 1, -(6,-) = 0 V j, г (7г, 2/) при у Є Aj, Vj(ir, у) при 2/ Є Aj+i — линейные, і/Дтг, 2/) - 0 при у Є Дг Мі ф j,j + 1 Vj = 1,..., TV - 2;
3) izjv-i(bjv-i) = 1, -1() = 0 Vj = 0,..., N - 2, I/JV-I(TT, 2/) = 1 при у Є ДАТ, VN-I{K, 2/) ПРИ 2/ Є Ддт-i — линейная, _і(тт, 2/) = 0 при 2/ Є До Vj = 0,..., iV - 2. Функции склейки в пространстве Л/г (ГІ2Л) будем искать в виде JV—1 г=0 / НЕ - У). г=0
Как и в работе [80], будем полагать, что на функции Л и /?, А; = 1,2, наложены условия Ад = Ад = — Л, / = / = -/ которые являются конечно-элементными аналогами компонент вектор-функции Л - (1.8) и определяют конечно-элементную вектор - функцию склейки Ад (Xh = (А , /? ), & = 1,2).
Сначала определим конечно-элементные пространства на Пд и на интерфейсе Г12/1, затем исследуем их свойства. Введем пространство УЛ(ОД = {ЛУД = (uh, vh);ukh = 1 Є vk\Qkh), 4 = vh\akevk\nkh),k = l,2}. В целях последующего анализа потребуется определить на пространстве Vh(&h) норму. Будем использовать для этого отображение wfc - Ык = (wfcJA + ІІК]ї/2ЛГі2)1/2, (1.35) где Шы = (Е llw,21AJ1/2 = (Е Е w e21)e)1/2, (1.36) Ы\\г/2Ат12 = h 1 2 \\[wh]\\o,r12. (1.37) Под функцией (вектор-функцией) [ ], в конечно-элементном случае, будем понимать разность значений функций (вектор-функций) на ГігППі иГі2ПП2.
Формулировка и доказательство аналога второй леммы Стрэнга
Теорема 2.1. Пусть w up - компоненты решения задачи (1.1), удовлетворяющие системе уравнений (1.4) и условию (1-5), и w Є Y(p)nH2(Q), р є І/ВДПЯ1 ), fL2(fi). Пусть (wb Ph, Ад) - приблиоюенное решение, удовлетворяющее системе уравнений (1.48) и условию (1-49), w Y/, /,), рд Є Хд(Т2д), Ад Є Мд(Гі2д). Тогда имеет место оценка \\-Wh\\vh + \\p-Ph\\xh С17 (_mf w - wAVfc + m p - qh\\xh+ t _ o (w, y?jQ + 6/,( , p) 4- dh(iph, Xh) f, 9?g К Г01Пч -і- sup —Qjj j, (/. IUJ где Ci7 - положительная константа. Доказательство. Пусть w — произвольная вектор-функция пространства Y/l(Q/l) такая, что w Ф w , а дд — произвольная скалярная функция пространства .ХД(ПД) такая, что q ф р и удовлетворяет условию (1.49).
Сложив правые и левые части неравенств (2.11), (2.12), соответственно, будем иметь (\\Ph - Qh\\xh + WA - WAVJ WA - WAvfc Cis(ah{w - WA, WA - wh) + 6A(WA - WA, p - qh) + ( f, WA - wh -( (w, WA - WA) - &A(WA - %, p) - dh(wh - wh, \))), (2.13) где Cis = max{Ci4, d5}.
Разделим обе части (2.13) на величину WA — WAHV И оценим правую часть неравенства своим модулем. Далее учитывая, что модуль суммы не больше суммы модулей ее слагаемых, устанавливаем Возьмем точную нижнюю грань, от обеих частей неравенства, по всем w/i Є Yh(Clh) и qh Є -Х /і(\), удовлетворяющим условию (1.49). В результате получим (2.10).
В этом пункте установим априорную оценку (w — w/j) в норме пространства V/ fift), на множестве вектор-функций пространства Y/ O/j), и (p—Ph) в норме пространства (0/,,), на множестве скалярных функций пространства Xh(flh), удовлетворяющих условию (1.49). Для этого предварительно докажем вспомогательные утверждения.
Лемма 2.2. Пусть w и р - компоненты решения задачи (1.1), удовлетворяющие системе уравнений (1-і) и условию (1.5), w Є Y(Q)nH2(Q), р Є X/R{Q,)f)H1(Q), и триангуляция на каждой подобласти П/-, k = 1,2, квазиравномерна. Тогда существует положительная постоянная Сч\, не зависящая от компонент решения w и р задачи (1.1) и h такая, что имеет место оценка Ы w-wVh + inf \\р-чЛхк С21К{Ы\1 Ы\и)- (2.17) Доказательство. 1. Пусть I/ w - интерполянт вектор - функции w на V/l(Q/l), a JhP - интерполянт скалярной функции р на Xhiflh)- Тогда, согласно работе [53], на подобластях fi& справедливы оценки (w - lh w)fiJififcfc + (р - Лр)кІІо,п C22 fc Л(w2 + bi,aj. (2.18)
2. Рассмотрим вектор-функцию разности (w — w ) в норме пространства Vh{ h)i где w/j - произвольная конечно-элементная вектор -функция этого пространства, и воспользуемся определением нормы (1.35) lw-whfrfc = w-% №+(w-wh21)Q2ft+[w- ll/2Ari2. (2.19)
Возьмем точную нижнюю грань от обеих частей (2.19) по всем элементам w/г, Є YhiP-h)- Имеем І& IW - »v» = «j& Hw - ИЇЯи+ +йШ w - w„]21A4 + Jnf I[w - Л]Ц5/2АГи. (2.20) Оценим каждое из слагаемых правой части (2.20): (1) для оценки первых двух слагаемых применим неравенство .inf. \\w-wh\\2lflkh (w-I,w)k2lfe, к = 1,2; (2.21) wftfc "ft (2) для оценки третьего слагаемого воспользуемся определением нор мы ЦІ/2,Л,ГІ2 (1-37) и неравенством (1.43) леммы 1.7; получим wSk[w ііії/w»= -14ll[w " %1 г- Г1 [w - 1Н11о,г12 h l (IKW - І )Піола+ +(w - І )п2о,г12)2 2/Г1 ((w - I w)ftl2 2+ +(w-Iftw)ta5iria) 2(C23)2/ 2(w )2. Таким образом, .in [w - %]21/2Л,Гіа 2 (С23)2 Л2 (Ы\ 2Д)\ (2.22) где положительная постоянная С2з = тах{Сзд, Сз]2}. Далее, сложив неравенства (2.21), (2.22), и подставив в (2.20), получим wjevfc "W " " И(w " ЬіИіА + (w - IfcwJIn.ll? + 2 (C23)2 u2 (w )2. (2.23) 3. Оценим норму разности (p — qh) в пространстве Х/ Пд), где g произвольная функция этого пространства, удовлетворяющая усло вию (1.49). Воспользуемся определением нормы lb - чь&ъ = lb - длііол + lb - gAV (2-24) Возьмем точную нижнюю грань от обеих частей (2.24) по всем элементам qh пространства Xh{Qh), удовлетворяющим условию (1.49):
Конечно-элементное представление
Использование данного алгоритма неприемлемо. Дело в том, что если первые строки матрицы X не содержат большого числа нулей, то она, в результате преобразований алгоритма (1)-(7), будет заполненной (значения элементов г-й строки передаются элементам j-й строки, если j г). В таком случае, во-первых, необходимо хранить в памяти компьютера огромное число элементов матрицы X, а во-вторых, потребуются большие временные затраты умножения вектора на эту матрицу, и как следствие - нецелесообразность использования такого подхода для построения переобуславливающей матрицы системы.
В монографии [96] предложен алгоритм построения матрицы X, который преодолевает отмеченные выше трудности. Через lev(rriij) определим уровень заполнения (level-of-fill) произвольного элемента mi j : foo, mitj = О, Тогда каждый раз в процедуре исключения Гаусса (для произвольной матрицы М), определенной по формуле TYli j = m i j TTli k TTlk ji уровень заполнения элемента rriij вычисляется следующим образом: lev(mij) = mm{lev(mij), lev{m k) + lev(mk,j) + 1}. Здесь {mij} представляет собой элементы нижнетреугольной матрицы, если i j и верхнетреугольной матрицы, если г j.
Таким образом, ILU(p) преобразование исходной матрицы состоит из процедуры исключения Гаусса, с одной поправкой: все элементы, чей уровень заполнения превосходит р - исключаются, то есть принимают значение нуль. Алгоритм для р = О рассматривался ранее при построении переобуславливателя А : Л = L U.
Далее приведем алгоритм построения матрицы X = L l ІЗ, когда уровень заполнения равен р : (1) для г = 1, ...,А4; (2) присваиваем г-й строке матрицы X значение г-й строки матрицы В : ХІ = ВІ\ (3) вычисляем уровень заполнения каждого элемента г-й строки матрицы X = {&i,j}: где bjj- элемент г-й строки, j-ro столбца матрицы В; (4) для & = 1,...,г-1, и ї к О; (5) выполняем преобразование над г-й строкой матрицы X : ХІ = Х{ — li:k Xk , (6) вычисляем уровень заполнения элементов г-й строки матрицы X (J = 1,...,) : lev(xitj) = mm{lev(xij), lev(xi k) + lev{xk,j) + 1}; (7) окончание цикла по к, (4); (8) если lev(xij) р, то Xij = 0; (9) в заключение поделим Х{ на диагональный элемент г-й строки матрицы L : Xi = Xi/kj, (10) окончание цикла по г, (1). Такой алгоритм в зарубежной литературе называется Row-oriented Level-p процедурой вычисления матрицы X.
Следует отметить (см. [88]), что если при построении переобуславливающей матрицы к Л был выбран алгоритм с уровнем заполнения равным р . то при построении переобуславливающей матрицы к S необходимо выбирать алгоритм с уровнем заполнения равным 2р-Ь 1. В качестве параметра р, для того чтобы матрица X, а значит и матрица S были разряженными и легко умножались на вектор, необходимо выбрать нуль.
Прямоугольная матрица X (У) полного ранга не является квадратной, размерность которой совпадает с размерностью матрицы В (С). Следует отметить, что эффективное решение системы с матрицей S невозможно, в отличие от эффективного умножения на эту матри-цу. Тем не менее, матрица S служит вспомогательной матрицей при построении переобуславливателя S. В работах [80], [1] изложен алгоритм нахождения решения системы, не используя непосредственного отыскания матрицы cS_1. В общем случае, пусть для решения системы уравнений
Ну построена переобуславливающая матрица Н, но она обладает неприятным свойством: эффективное решение системы с матрицей Н (нахождение обратной - Я-1) невозможно. Тогда следует построить такую переобуславливающую матрицу к Н - Н, для которой Н будет вспомогательной с тем свойством, что ее можно легко умножать на вектор. Процесс решения системы с матрицей Н будет выглядить как внутренняя итерационная процедура.
В настоящей работе в качестве такой процедуры, при получении результата умножения вектора пунктов (2) и (6) алгоритма (1)-(16) на матрицу S-1, предлагается использовать GMRES-метод с вспомогательной матрицей S. Внутренняя итерационная процедура содержит в качестве своих итераций также алгоритм (1)-(16) с той особенностью, что К = К г = Е.
С целью численного подтверждения теоретических результатов об аппроксимации, установленных в главе 2 (теоремы 2.2 и 2.3), были проведены расчеты задачи Стокса с кусочно-постоянным коэффициентом кинематической вязкости в дивергентно-градиентной части уравнения.
В модельных задачах рассматривалась система дифференциальных уравнений (3.1), для которой известно точное решение (w, р). Функция давления р выбиралась так, чтобы было выполнено условие j pdxdy = О, а. а вектор-функция скорости w = (и, v) удовлетворяла уравнению несжимаемости и краевым условиям задачи div w = 0 в Q, w п = 0 на Г, w т = qi на Г\, г = 1,..., 4, где пит — единичные нормальный и касательный векторы к границе. В численных экспериментах ао(х, у) имела вид V UJ I 02, (ж, У) ЄП2.
В таком случае, подставляя точное решение (w, р) в систему (3.1), при каждом значении положительного параметра a i, находили вектор-функцию f = (/і, /2) правой части. Для проведения расчетов был создан комплекс программ на языке Си, реализующий алгоритм нахождения приближенного обобщенного решения задачи по МКЭ, представленный в пункте 3.3. Разность между точным и приближенным решениями для компонент вектор-функции скорости определялась в сеточных нормах пространств Li и Н1, а для скалярной функции давления - в сеточной норме пространства L2, на подобластях П/ , к = 1,2. Одна из рассмотренных характеристик, для компонент вектор-функции скорости - поведение погрешности в норме пространства L i на полосах в окрестности интерфейса Гі2 между подобластями. Ширина полосы, примыкающей к интерфейсу Г\2, равна с каждой стороны. При этом полосу в к-й подобласти обозначили через Ek,k = 1, 2.
Для того, чтобы ввести норму сеточного пространства L на Сік потребовалось: (1) разбить Сік (квадрат) вертикальными и горизонтальными линиями на элементарные квадраты Е\ со стороной 2/, их совокупность обозначить через Е ; (2) множество точек, являющиеся серединами сторон квадратов Б\ , объединить в SQ , а подмножество его точек не лежащих на Г, которые согласно пункта 3.2 (подпункта 3.2.4) - узлы аппроксимации, объеДИНИТЬ В SQ , (3) разбить исходную подобласть Сік на непересекающиеся треуголь ники U (см. рис. 1) так, чтобы каждый треугольник в качестве своих —(к) (к) вершин имел представителей множества S0 .
