Содержание к диссертации
Введение
Глава 1. Постановка задачи 27
1.1. Математическая модель движения вязкой несжимаемой жидкости: уравнения Навье - Стокса, уравнения Стокса 27
1.2. Математическая модель движения идеальной жидкости: уравнения Эйлера 29
1.3. Вариационные постановки 29
Глава 2. CLASS Построение дискретных аналогов вариационных постановок CLASS 32
2.1. Триангуляция расчетной области, построение двойственных сеток 32
2.2. Дискретные аналоги вариационных постановок 34
2.3. Пространства интерполяционных функций в смешанных постановках 38
2.4. Алгоритмы расчета поля скоростей и давления
2.4.1. Смешанный метод конечных элементов 41
2.4.2. Алгоритм Удзавы 45
2.4.3. Стабилизированные методы конечных элементов 46
2.4.4. Модификация конвективных членов в противопотоко-вых схемах 49
2.4.5. Противопотоковая схема с вычислением неизвестных в узлах конечных элементов 51
2.4.6. Противопотоковая схема с вычислением неизвестных в серединах ребер конечных элементов
Глава 3. CLASS Реализация алгоритмов моделирования, результаты вычислительных экспериментов CLASS 59
3.1. Структура комплекса программ 59
3.1.1. Структура данных комплекса программ 59
3.1.2. Основные модули комплекса программ 63
3.2. Тестирование методов 72
3.2.1. Численное решение уравнений конвективно-диффузионного типа в двумерной области. Проти-вопотоковые схемы 73
3.2.2. Численное решение уравнений конвективно-диффузионного типа в трехмерной области 77
3.2.3. Численное решение нестационарных уравнений конвективно-диффузионного типа в трехмерной области 79
3.2.4. Численное решение уравнений Навье - Стокса и Стокса в двумерной области 82
3.2.5. Численное решение уравнений Навье - Стокса, Стокса и Эйлера в трехмерной области 87
3.2.6. Численное решение задачи о каверне 89
3.3. Моделирование течений в криволинейных каналах и каналах переменного сечения 92
3.3.1. Моделирование течения в двумерном криволинейном канале с разворотом потока на 270 92
3.3.2. Моделирование течения в криволинейном канале с разворотом потока на 180 93
3.3.3. Моделирование течений вязкой несжимаемой жидкости в каналах переменного сечения
Заключение
- Математическая модель движения идеальной жидкости: уравнения Эйлера
- Дискретные аналоги вариационных постановок
- Стабилизированные методы конечных элементов
- Численное решение уравнений конвективно-диффузионного типа в двумерной области. Проти-вопотоковые схемы
Введение к работе
В современной науке для изучения самых разнообразных природных явлений, технологических процессов и экологических проблем широко используется математическое моделирование [36]. Исходный объект заменяется математической моделью, которая исследуется средствами вычислительной математики. Математическое моделирование сочетает в себе многие достоинства как теории, так и эксперимента [36]. Работа с моделью объекта дает возможность относительно быстро и без существенных затрат исследовать его свойства и поведение в различных ситуациях (преимущества теории). В то же время вычислительные эксперименты с моделями объектов позволяют подробно и глубоко изучать объекты в достаточной полноте, недоступной теоретическим подходам (преимущества эксперимента). Совместное использование физического и вычислительного экспериментов в исследовании объекта позволяет, с одной стороны, уменьшить количество натурных дорогостоящих измерений, а с другой стороны - провести верификацию и усовершенствование математических моделей.
Широкий круг проблем, стоящих перед современной наукой и техникой, связан с изучением движения жидкости. Для описания явлений используются модели стационарных и нестационарных течений, вязкой и невязкой жидкости, пограничных слоев и другие [26, 42].
В рамках выбранной математической модели, состоящей из системы дифференциальных уравнений и краевых условий (и начальных - для нестационарной задачи), определяется поведение течения жидкости в тех или иных условиях. Эффективным способом решения дифференциальных уравнений являются численные методы. В отличие от аналитических методов, где практически для каждой задачи разрабатываются свои самостоятельные приемы решения [33], численные методы отличаются большей универсальностью и
применимы для исследования широкого класса задач.
При подборе методов решения необходимо учитывать особенности исследуемых течений и соответствующих математических моделей. В диссертационной работе рассматриваются стационарные, внутренние, однофазные, химически однородные течения. Движение вязкой несжимаемой жидкости описывается уравнениями Навье - Стокса. При создании комплекса программ необходимо наличие не только программных реализаций сложных моделей, но и реализации более простых приближенных моделей, с помощью которых можно получать предварительное представление о характере течения. В случае больших значений кинематической вязкости или при малых значениях вектора скорости решение уравнений Навье - Стокса может быть заменено решением уравнений Стокса. Решения уравнений Эйлера, описывающие течение идеальной несжимаемой жидкости, могут служить начальным приближением для итерационных методов решения уравнений Навье - Стокса. Поэтому в диссертационной работе рассмотрен и включен в комплекс программ набор моделей, что помогает более полно исследовать явление.
Существуют следующие математические формулировки перечисленных уравнений, описывающих движение несжимаемой жидкости [40, 125]:
в естественных переменных вектор скорости - давление (и - р формулировка);
в переменных векторный потенциал - вектор вихря (ф - Со формулировка);
в переменных вектор скорости - вектор вихря (и- ш формулировка).
Каждая постановка имеет свои преимущества и недостатки. Исторически описание в переменных функция вихря - функция тока было весьма популярно при расчете двумерных течений вязкой несжимаемой жидкости [28, 30, 39, 40, 94, 125]. Трудности обобщения этого подхода на трехмерный случай связаны с установлением граничных условий для трехмерного анало-
га функции тока - векторного потенциала [30, 40, 125]. Это является одной
—*
из причин, по которой ф - ш формулировки использовались относительно редко при численном моделировании трехмерных течений. Однако, переход к зависимым переменным фиш имеет ряд преимуществ по сравнению с и - р подходом. Так, для ф - ш формулировок уравнение неразрывности на дифференциальном уровне выполняется автоматически, а на разностном -при подходящей аппроксимации компонент вектора скорости и использовании разнесенных сеток. С введением векторного потенциала (функции тока в двумерном случае) давление не входит в число зависимых переменных, поэтому на промежуточных итерациях давление исключается из вычислительного процесса и при необходимости может быть восстановлено после сходимости
—>
итераций для ф и ш. Примеры вычислительных алгоритмов, основанных на ф - ш формулировках, представлены в [40, 121, 127].
В работах [74, 104, 105] используется и — w формулировка, в которой поле скоростей вычисляется путем решения уравнений Пуассона для каждой компоненты вектора скорости. В статье [77] на каждом итерационном шаге компоненты вектора скорости находятся из решения уравнения неразрывности и соотношений, задающих связь между вектором скорости и вектором вихря. В трехмерном случае u — w подход менее эффективен, чем описание в естественных переменных, если только движение вихря, в особенности нестационарного, не представляет особого интереса. Более подробный обзор данных подходов можно найти в книге [40].
В диссертационной работе рассмотрены математические модели в естественных переменных вектор скорости - давление.
Система уравнений Навье - Стокса состоит из уравнений неразрывности и движения. Преобладание конвективного или диффузионного переноса в течении жидкости определяется числом Рейнольдса. Уравнения движения и неразрывности образуют систему связанных уравнений с неизвестными ве-
личинами - компонентами вектора скорости и давления. Для несжимаемых течений важно влияние перепадов давления на поле скоростей, способ же расчета давления не очевиден, поскольку давление не является термодинамической переменной и для его определения не существует отдельного уравнения. Отметим, что при использовании моделей в естественных переменных на уровне вычислительной схемы возникают трудности корректного учета условия, накладываемого на дивергенцию аппроксимируемого поля скоростей, так называемого дивергентного условия.
Численное решение уравнений Навье - Стокса, Стокса и Эйлера в естественных переменных предполагает решение следующих задач, определяющих структуру комплекса программ [14, 17, 28]: пространственной дискретизации уравнений; выбора базисных функций для вектора скорости и и давления р, гарантирующих существование решения (и,р); учета нелинейности, связанной с конвективными членами в уравнении движения; выбора алгоритма расчета неизвестных вектора скорости и давления; решения системы линейных алгебраических уравнений (СЛАУ).
На сегодняшний день ведущие позиции при численном решении дифференциальных уравнений занимают сеточные методы. В зависимости от способа построения приближенного решения можно выделить следующие классы сеточных методов, часто используемые в вычислительной гидродинамике [18, 24, 28, 39, 40, 41, 46]: дифференциальные и вариационные.
Отличительной характеристикой дифференциальных методов является непосредственная аппроксимация исходного дифференциального уравнения, либо соответствующих интегробалансных соотношений, с помощью разностных аналогов производных. В этот класс входят метод конечных разностей (МКР) и метод конечных объемов (МКО).
Идея МКР состоит в замене расчетной области множеством дискретных точек (сеткой) и построении дискретизации дифференциальных операторов
в узлах этой сетки. Существует более десятка конечно-разностных схем для решения уравнений Навье - Стокса, отличающихся различными признаками [18, 24, 28, 35, 39, 40, 41, 46, 92]. Особенностью этих схем является наличие различных сеточных параметров, рациональный выбор которых трудоемок [28]. К недостаткам метода относятся также сложность построения конечно-разностных операторов, аппроксимирующих уравнения Навье - Стокса на произвольном сеточном шаблоне, учет краевых условий и адаптация регулярных сеток к областям с произвольной криволинейной границей [28]. Среди преимуществ МКР - низкие вычислительные затраты и возможность применения к широкому спектру решаемых задач. Применение в контексте метода конечных разностей разнесенных („шахматных") сеток для вычисления компонент вектора скорости и давления, впервые предложенное в 1965 г. Харлоу и Велч [81], позволяет связать значения вычисляемых неизвестных в соседних точках. Осцилляции, как правило, возрастают при увеличении числа Рейнольдса, поскольку диссипативные члены, посредством которых осуществляется связь значений компонент вектора скорости в соседних точках, в этом случае малы. Случай разнесенных расстановок переменных эквивалентен применению разных порядков базисных функций для обеспечения условия существования и единственности решения (и,р) [30, 91].
В отличие от МКР, метод конечных объемов [34] связан с построением дискретизации интегральных формулировок на заданной сетке. Преимуществом конечно-объемных аппроксимаций являются простота реализации и возможность работы с неструктурированными сетками. Для построения устойчивых конечно-объемных аппроксимаций в ряде работ используются пространственные „шахматные" сетки. В работах [ПО, 116] была предложена разнесенная дискретизация для неструктурированных сеток. В работе [110] для компонент вектора скорости используются одни и те же контрольные объемы, а вклад градиента давления в уравнение движения определяется с помощью
специальной реконструкции давления на разнесенной сетке.
Необходимо отметить, что использование „шахматных" сеток в МКР и МКО позволяет повысить устойчивость схем и на уровне вычислительной схемы корректно учесть дивергентное условие.
Другой класс сеточных методов - методы, основанные на использовании вариационных постановок, эквивалентных исходной дифференциальной краевой задаче. В этот класс входят различные варианты метода конечных элементов (МКЭ) [28, 29, 37, 38, 45, 51, 57, 58, 59, 60]. Одними из важных достоинств МКЭ для задач гидродинамики является эффективность конструирования конечно-элементных аппроксимаций на неструктурированных сетках и возможность практически полной автоматизации всех этапов решения задачи - от построения сетки до сборки глобальной матрицы результирующей СЛАУ [45, 51]. Высокую вычислительную эффективность МКЭ обеспечивает выбор функций кусочно-полиномиального вида [29, 38].
В последние годы создано множество методов моделирования течений несжимаемой жидкости с использованием МКЭ. Важным критерием эффективности таких схем является точность учета дивергентного условия. Наиболее распространенным подходом, позволяющим учесть дивергентное условие, является введение смешанных вариационных постановок [15, 128, 129]. Под термином „смешанные постановки" понимаются такие постановки, в которых аппроксимация происходит на паре специально подобранных пространств. Пары пространств, используемых в смешанных вариационных постановках, необходимо подбирать в соответствие с условием Ладыженской - Бабушки - Брецци (ЛББ), которое является одними из необходимых условий существования и единственности решения [15, 27, 39]. Интерполяционные функции в смешанных постановках, удовлетворяющие условию ЛББ, называются div-устойчивыми (равномерно устойчивыми) [15].
Имеются и другие подходы [40] к удовлетворению дивергентного условия.
Например, появление методов штрафа дало начало новому направлению в области моделирования течений вязкой несжимаемой жидкости методом конечных элементов - численным алгоритмам, допускающим произвольные сочетания интерполяционных полиномов решения. Соответствующие методы называются методами, „обманывающими" условие ЛББ. Применение этих методов [39, 56, 57, 64] позволяет исключить явное наличие давления в уравнении движения.
К классу методов с произвольным сочетанием интерполяционных полиномов решения также относятся стабилизированные схемы Галеркина [61, 83, 84, 97, 98, 120]. Назначение параметров стабилизации в этих методах - регулировать величину искусственной вязкости, добавляемой в уравнения для подавления больших нефизических осцилляции. Способы задания параметров стабилизации, зависящих от пространственных и временных характеристик, приводятся в [120]. Получение оптимального значения численной диффузии в направлении линий тока также исключает возможность возникновения осциллирующего поля давления.
Отметим также, что хорошие результаты дает алгоритм Удзавы [51, 52, 64], традиционно используемый при решении задачи о седловой точке. Для улучшения свойств сходимости алгоритма Удзавы в уравнение неразрывности вводится слагаемое, содержащее давление и штрафной параметр [52, 64, 117, 118, 120].
Противопотоковые схемы. При больших числах Рейнольдса Re преобладающим является конвективный перенос, что требует устойчивых алгоритмов учета конвективных членов в уравнении движения. Существуют различные итерационные процедуры, используемые для линеаризации конвективных членов в зависимости от значений числа Рейнольдса.
Для течений, в которых доминирует конвекция, использование центральных разностей приводит к появлению в решении сильных нефизических ОС-
цилляций [30]. При попытке стабилизировать решение путем дискретизации конвективных членов двухточечными разностями против потока или суммами с весами центрально-разностных выражений и разностей против потока получается неточное решение, особенно в случаях совпадения локального направления вектора скорости с ориентацией (левой, правой) сетки и больших локальных градиентов вектора скорости.
Более точное решение получается, если для аппроксимации конвективных членов используются разности высокого порядка точности против потока, подобные четырехточечным формулам. В комбинированной схеме [30] в зависимости от абсолютной величины сеточного числа Пекле сочетаются центрально-разностная и противопоточная схемы аппроксимации конвективных членов уравнения. Схема квадратичной интерполяции против потока [96] (схема Леонарда) имеет второй порядок точности, но менее устойчива при расчетах, чем комбинированная схема. В работе Дэвиса и Мура [65] для решения нестационарных уравнений Навье - Стокса используются многомерные разности третьего порядка точности против потока, которые являются обобщением одномерной квадратичной интерполяции против потока Леонарда [96]. Схема с косыми разностями против потока Рейсби позволяет учитывать направление вектора скорости потока по отношению к грани контрольного объема [110].
Далее перечислим противопотоковые схемы на базе методов конечных элементов и конечных объемов. Наиболее исследованными и часто применяемыми в случае неструктурированных сеток являются противопотоковые схемы, основанные на методе конечных объемов с расчетом неизвестных в центрах ячеек [53, 69]. Несмотря на то, что грани элементов разбиения не параллельны координатным осям, данные схемы в большинстве случаев являются одномерными [111]. Расчеты по подобным схемам не воспроизводят многомерную структуру потока и обладают чрезмерной численной диффузией [69]. Для
построения противопотоковых схем второго порядка необходимо расширение шаблона, что для неструктурированных разбиений приводит к усложнению структур данных [53, 69].
Противопотоковые схемы с расчетом неизвестных в узлах сетки и серединах ее ребер немногочисленны [99, 108, 111]. В ряде случаев противопото-ковый принцип аппроксимации сводится к использованию одного значения скалярной субстанции (в узле симплекса, лежащего против потока) [93], либо двух взвешенных значений (на концах ребра симплекса, лежащих против потока) [99, 111].
Схема, учитывающая направления потока (Flow Oriented Upwind Scheme), использует преимущество расчета неизвестных в узлах - возможность построения асимметричных функций профиля. Но расчеты по данной схеме признаны неудовлетворительными и итерационные процессы часто расходятся [99], поскольку схема не обладает свойством положительности.
Конечно-элементный аналог схем против потока, предложенный Кристи [62], состоит в использовании функций формы, модифицированных с тем, чтобы вклады каждого из элементов распределялись в зависимости от доминантного направления течения на элементе, а не равномерно по его узлам. Однако, предложенные модификации вносят чрезмерную численную диффузию в направлении перпендикулярном вектору средней скорости на элементе.
В работах Келли и соавторов [89] в методе Петрова - Галеркина со стабилизацией оператора конвекции в направлении линии тока [83, 84, 122] предложено использовать функции формы, искривленные в направлении средней скорости на элементе. Другим подходом является метод стабилизирующей тензорной вязкости [78], позволяющий работать с конечно-элементными пространствами младших порядков. Выражения для стабилизирующей тензорной вязкости учитывают направления средней скорости на элементе. Анализ этого метода показывает [78], что в нем не отражается многомерная струк-
тура потока, поэтому привносимая численная диффузия велика.
Другим эффективным подходом к решению задач с преобладанием конвекции является смешанный метод конечных элементов/объемов, в котором используются метод Галеркина с симметричными тестовыми функциями для самосопряженной части дифференциальных операторов и противопотоковые схемы МКО - для несимметричной их части [67, 69].
Метод конечных суперэлементов (МКСЭ) [19, 20, 21] был предложен как метод дискретизации дифференциальных уравнений, учитывающих мелкомасштабные (относительно шага сетки) неоднородности задачи (стержни в ядерном реакторе, инородные включения в композитных материалах и т. п.) или, если явных неоднородностей нет, то присутствие различных пространственных масштабов решения. Основные вычислительные затраты МКСЭ связаны с построением базиса на элементе разбиения. С другой стороны, при использовании схем высокого порядка МКСЭ включает в себя процедуру локального исключения неизвестных, ассоциированных с внутренними узлами элемента. Такое исключение приводит к уменьшению размерности результирующей линейной системы уравнений, улучшению ее обусловленности и уменьшению вычислительных затрат на ее решение.
В разрывном методе Галеркина [50, 63, 66, 70] решение аппроксимируется на конечно-элементной сетке кусочно-полиномиальными функциями, разрывными на межэлементной границе. Главной особенностью этих методов является то, что решение удовлетворяет закону сохранения локально на каждом элементе (поэлементно), степени свободы решения ассоциируются с элементами пространственной дискретизации, а не с вершинами элементов как в классическом методе Галеркина. Несмотря на большой интерес, проявляемый к методу и его модификациям, для него еще не создано единой вычислительной технологии [70].
Основная особенность класса противопотоковых конечно-элементных
схем, использующих аналог конечно-объемной двойственной сетки, состоит в построении ячеек двойственной сетки, на которых используется кусочно-постоянный оператор модификации конвективных членов. В данном классе выделяются противопотоковые схемы с расчетом неизвестных в вершинах симплексов [80] и в серединах ребер симплексов [103] первичного разбиения расчетной области.
Методы аппроксимации по времени. Несмотря на то, что диссертационная работа посвящена стационарным течениям, проведено исследование разработанных вычислительных схем на нестационарных задачах. По этой причине рассмотрим схемы аппроксимации по времени, используемые при численном решении дифференциальных уравнений. В литературе (см., например, [16, 35]) схемы аппроксимации по времени классифицируются на явные, неявные и полунеявные.
К классу явных схем относятся, например, схемы Адамса - Башфорта и Рунге - Кутта [122]. В явных схемах эллиптическая часть исходного уравнения вычисляется на предыдущих временных слоях, и матрица массы имеет ненулевые элементы только на главной диагонали. В данном случае нет необходимости решать СЛАУ на каждом временном слое, так как компоненты вектора решения явным образом выражаются через компоненты вектора решения на предыдущем временном слое. Диагональную матрицу массы можно получить при использовании подходящей дискретизации. В работе [127] предлагается набор базисных функций, которые позволяют получить диагональную матрицу массы. Показано, что при использовании таких базисных функций скорость расчета увеличивается в три раза, однако количество базисных функций и размерность СЛАУ при этом также увеличивается в три раза.
Если матрица массы, получаемая при дискретизации, не является диагональной, то выполняется ее диагонализация [113], состоящая в поэлементном
преобразовании матрицы массы, при котором диагональным элементам присваиваются значения суммы всех элементов в строке, а внедиагональные элементы полагаются нулевыми. В диссертационной работе поэлементное преобразование матрицы массы с целью ее диагонализации применяется для обращения матрицы массы в вычислительных схемах расчета вектора скорости и давления.
Явные схемы характеризуются низкими вычислительными затратами. Общим недостатком всех явных схем являются строгие ограничения на соотношение шагов по пространству и времени: с уменьшением размеров ячейки дискретизации максимально допустимый для устойчивого счета шаг по времени должен быть уменьшен пропорционально квадрату размера ячеек дискретизации. Поэтому часто используются неявные и полунеявные схемы.
Неявные схемы характеризуются вычислением матрицы жесткости на текущем временном слое, следовательно, на каждом временном слое требуется решение СЛАУ. Однако неявные схемы более устойчивы, чем явные, и порой оказываются более экономичными по времени счета, поскольку для устойчивости явных схем требуется большее число временных слоев. К классу неявных схем относятся, например, схемы с обратным дифференцированием и схемы типа Адамса - Моултона, среди которых можно выделить неявную схему Эйлера и схему Кранка - Николсон [25, 122].
Полунеявные схемы сочетают в себе достоинства явных и неявных схем. К этим схемам можно отнести методы типа „предиктор-корректор", в которых в качестве начального приближения на каждом временном шаге неявной схемы служит решение на этом же шаге, полученное с помощью явной схемы. В работе [55] предложен новый метод подобного класса, использующий для решения СЛАУ только некоторое фиксированное число итераций метода обобщенных минимальных невязок. Число итераций обычно выбирается небольшим, что значительно сокращает вычислительные затраты на реше-
ниє СЛАУ. Полученная схема может рассматриваться как явная, поэтому приходится накладывать ограничения на шаг по времени.
Одним из способов повышения эффективности схем аппроксимации по времени является метод Тейлора - Галеркина [48, 71]. Эта схема основана на разложении искомой функции в ряд Тейлора по времени, последующей замене временных производных пространственными, и применении к полученной системе метода Галеркина. На основе метода Тейлора - Галеркина предложен ряд явных и неявных схем [71]. Изначально метод предлагался для задач конвективного переноса, но в последствии был применен и в других приложениях [71].
Несмотря на разнообразие подходов к аппроксимации дифференциальных задач, четких рекомендаций по выбору оптимальной схемы для конкретной задачи не существует. Этот выбор должен делаться на основании многих факторов, таких как физические и геометрические особенности решаемой задачи, метод пространственной дискретизации, вычислительные ресурсы и другие.
Решение систем линейных алгебраических уравнений. Решение систем линейных алгебраических уравнений, получаемых в результате дискретизации нелинейных задач, является важным этапом в математическом моделировании. Несамосопряженность оператора исходной дифференциальной постановки порождает ряд сложностей при применении многих численных методов. В сеточных методах несамосопряженным операторам, как правило, соответствуют несимметричные матрицы СЛАУ, что затрудняет решение таких систем.
Применение обобщенных методов сопряженных градиентов [23] или методов минимальных невязок [23] для несимметричных систем уравнений фактически эквивалентно трансформации Гаусса (то есть решению системы с матрицей ААТ), что в общем случае приводит к понижению эффективности
итерационного процесса (по этим вопросам в работе [23] приводятся ссылки на обширную литературу).
Более эффективным и устойчивым среди итерационных методов решения таких систем уравнений является класс проекционных методов с использование проектирования на подпространства Крылова [23, 100, 113, 114, 123]. Общий подход к построению проекционных методов заключается в применении условия Петрова - Галеркина, которое означает ортогональность вектора невязки некоторому n-мерному подпространству. Наиболее известными являются следующие методы [23, 100, ИЗ, 114, 123]: метод сопряженных градиентов и метод полной ортогонализации, применяемые для решения СЛАУ с симметричной матрицей, а также метод обобщенных минимальных невязок, метод бисопряженных градиентов и метод бисопряженных градиентов со стабилизацией - для случая несимметричных матриц.
О предмете и содержании диссертации. Диссертация состоит из введения, трех глав, заключения, списка использованных источников (131 наименование), приложения. Работа изложена на 138 страницах, включая 23 таблицы и 40 иллюстраций.
Цель работы состоит в моделировании внутренних течений несжимаемой жидкости с преобладанием конвективного переноса в каналах сложной геометрии, являющихся элементами различных технических устройств.
Методы исследования: методы вычислительной математики; сравнительный анализ результатов моделирования и имеющегося аналитического решения; расчеты на последовательности сгущающихся сеток с последующим анализом сходимости к аналитическим решениям.
В первой главе сформулированы постановки задач о внутренних стационарных течениях вязкой и идеальной несжимаемой жидкости в естественных переменных вектор скорости - давление. Глава состоит из трех параграфов.
В параграфе 1.1 приведена математическая постановка задачи об устано-
вившимся движении вязкой несжимаемой жидкости через ограниченную односвязную область П, которая заключается в определении вектора скорости и и давления р, удовлетворяющих в области Q уравнениям Навье - Стокса и краевым условиям на границе области: заданному вектору скорости на входе в область, условиям непротекания на твердых стенках и заданной нормальной компоненте тензора напряжения на выходе.
В случаях больших значений кинематической вязкости или при малых значениях вектора скорости решение уравнений Навье - Стокса может быть заменено решением уравнений Стокса с соответствующими краевыми условиями на границе области.
В параграфе 1.2 приведена математическая постановка задачи об установившимся протекании идеальной несжимаемой жидкости через ограниченную односвязную область О, которая заключается в определении вектора скорости и и давления р, удовлетворяющих в Q, уравнениям Эйлера и краевым условиям на границе области: заданному вектору скорости на входе в область, условиям непротекания на твердых стенках и заданной нормальной компоненте вектора скорости на выходе.
В параграфе 1.3 для поставленных задач приведены вариационные постановки в форме Галеркина, эквивалентные исходным [39, 57, 79].
Вторая глава посвящена построению дискретных аналогов вариационных постановок в форме Галеркина, полученных в первой главе. Разработана вычислительная схема расчета неизвестных компонент вектора скорости и давления в гипервекторе, позволяющая единообразно моделировать стационарные внутренние течения вязкой несжимаемой жидкости (уравнения Навье - Стокса, Стокса) и идеальной несжимаемой жидкости (уравнения Эйлера). Реализованы стабилизированные схемы Галеркина для моделирования трехмерных течений на тетраэдральном разбиении. Линеаризация конвективных членов выполнена противопотоковыми схемами с расчетом неизвестных в уз-
лах конечных элементов и серединах их ребер. Для аппроксимации потоков через ребра двойственной сетки предложено использовать формулы интегрирования базисных функций по отрезкам медиан конечного элемента.
Моделирование течения идеальной жидкости требует задание нормальной компоненты вектора скорости на границе области. В вариационных постановках уравнений Эйлера при использовании базисных функций пространств Тейлора - Худа в явном виде краевое условие для нормальной компоненты вектора скорости не задается. В диссертационной работе для учета краевого условия предложено построение матрицы преобразования декартовых компонент вектора скорости в тангенциальную и нормальную компоненты скорости. Глава состоит из пяти параграфов.
В параграфе 2.1 введены (для реализации противопотоковых схем учета конвективных членов) двойственные по отношению к первичной сетке разбиения: 1) разбиение на ячейки, построенные вокруг вершин симплексов первичной сетки; 2) разбиение на ячейки, построенные вокруг середин сторон симплексов первичной сетки.
В параграфе 2.2 построены дискретные аналоги полученных в первой главе вариационных постановок. Важным критерием вычислительных схем для моделирования течений вязкой несжимаемой жидкости является эффективность учета дивергентного условия. В связи с этим требуется введение специальных постановок и выбор функциональных пространств в вариационных постановках. Для корректного учета уравнения неразрывности построены смешанные [57] и стабилизированные конечно-элементные постановки [84, 120]. В стабилизированных конечно-элементных постановках для компонент вектора скорости и давления используются базисные функции одного порядка - кусочно-линейные. В смешанных конечно-элементных постановках базисные функции одного порядка приводят к неустойчивости решения, поэтому выбору базисных функций для смешанных конечно-элементных по-
становок посвящен следующий параграф второй главы.
В параграфе 2.3 отмечены допустимые дискретные пространства базисных функций для смешанных конечно-элементных постановок. Пространства выбираются в соответствии с условием Ладыженской - Бабушки - Брецци [27], которое является одним из необходимых условий существования и единственности решения вариационной постановки. В работе [15] показано, что пары пространств Тейлора - Худа и Крузея - Равьяра соответствуют условию Ладыженской - Бабушки - Брецци.
В параграфе 2.4 представлены алгоритмы расчета поля скоростей и давления, основанные на смешанных и стабилизированных конечно-элементных постановках.
В смешанных конечно-элементных постановках несжимаемое слагаемое штрафуется пропорционально давлению, что дает возможность отдельно получить уравнения для расчета компонент вектора скорости и уравнение для вычисления давления. В алгоритмах для обращения матрицы массы выполнена ее диагонализация - преобразование элементов матрицы, при котором диагональным элементам присваиваются значения суммы всех элементов в строке, а внедиагональные элементы полагаются нулевыми.
В качестве стабилизирующего оператора в стабилизированных конечно-элементных постановках используется конвективная часть исходного оператора. В диссертационной работе параметр стабилизации определяется в каждой ячейке ее размерами и коэффициентами уравнений. Назначение этого параметра - регулировать величину искусственной вязкости, добавляемой в уравнения для подавления больших нефизических осцилляции.
В дискретных вариационных постановках уравнений Эйлера для задания нормальной компоненты вектора скорости на границе области и п|хєГ = h, (n - вектор нормали к границе Г, h - заданная функция, если h — О, то задано условие непротекания на твердых стенках) введена матрица преобразования
декартовых компонент вектора скорости в тангенциальную и нормальную компоненты.
Для учета конвективных членов реализованы конечно-элементные проти-вопотоковые схемы с расчетом неизвестных в узлах и серединах ребер конечных элементов. При построении противопотоковых схем выполнена модификация дискретного аналога конвективных членов с помощью кусочно-постоянного объединяющего оператора. Вклады в матрицу конвективных членов вносят значения в узлах, которые лежат вверх по потоку. Для аппроксимации потоков через границы ячеек двойственной сетки предложено использовать формулы интегрирования базисных функций по отрезкам медиан конечного элемента.
Таким образом, в диссертационной работе для корректного учета уравнения неразрывности используются смешанные и стабилизированные постановки и соответствующие им вариационные задачи. Учет конвективных членов выполнен противопотоковыми схемами с расчетом неизвестных в узлах и серединах ребер симплексов сетки и методом Петрова - Галеркина со стабилизацией оператора конвекции в направлении линии тока.
Третья глава посвящена описанию разработанного комплекса программ, вычислительным экспериментам и анализу их результатов. Глава состоит из трех параграфов.
В параграфе 3.1 описана структура комплекса программ, структура данных и основные модули комплекса программ. Разработана структура данных для хранения элементов матриц СЛАУ, в которых неизвестный вектор является гипервектором, содержащим компоненты вектора скорости и давления. Проведено исследование структуры глобальных матриц и их размерности для разных вариационных постановок. Процесс решения задачи разделен на этапы, что позволяет в комплекс программ без особых затрат включать новые вычислительные алгоритмы и расширять круг задач.
Параграф 3.2 посвящен тестированию реализованных методов на ряде модельных задач с гладким решением на последовательности вложенных пространственных сеток. На двумерной задаче конвекции-диффузии с доминированием конвекции выполнено сравнение следующих методов: метода Петрова - Галеркина стабилизации давления со стабилизацией оператора конвекции в направлении линии тока, метода конечных элементов/конечных объемов, МКЭ и конечно-элементных противопотоковых схем с расчетом неизвестных в узлах и серединах ребер треугольников.
Выполнено сравнение точности решения уравнений конвективно-диффузионного типа для трехмерной задачи, полученного методом конечных элементов и стабилизированным методом конечных элементов. Для этой же тестовой задачи решено нестационарное уравнение конвективно-диффузионного типа. Аппроксимация по времени выполнена схемой Кранка - Никол сон.
Программная реализация алгоритмов расчета вектора скорости и давления верифицирована на модельных двумерных и трехмерных уравнениях На-вье - Стокса (в стационарном и нестационарном случаях), Стокса и Эйлера. Для методов, использующих функции штрафа, определены оптимальные значения штрафного параметра.
Численно решена задача о каверне, наблюдается хорошее совпадение полученных результатов с опубликованными данными [73].
В параграфе 3.3 приведены результаты моделирования течений в двумерных и трехмерных криволинейных каналах с разворотом потока на 180 и 270.
Расчеты в криволинейном канале с разворотом потока на 180 выполнены при задании сдвига скорости на входе. При увеличении сдвига скорости на входе характер течения существенно не меняется, увеличивается интенсивность вторичного течения. Сравнение результатов расчетов, полученных в
переменных вектор скорости - давление методом конечных элементов и в переменных векторный потенциал - вектор вихря методом конечных разностей [12, 47], показало хорошее согласование. Полученные результаты качественно согласуются с выводами работ [47, 95] и заключаются в следующем: перенос частиц жидкости в окрестности верхней стенки трубы происходит в сторону внешней стенки, а около нижней - в сторону внутренней.
Трубопроводы, используемые в промышленности, представляют собой систему, состоящую из трубы („русла") и различных конструктивных включений (заслонки, резкие сужения, расширения и другое). Проведена серия расчетов для каналов с расширением и сужением русла. Показано влияние соотношения узкого и широкого сечений на погрешность вычислений при измельчении сетки.
В приложениях 1-3 приведены компоненты локальных матриц и векторов правых частей, вычисленные на треугольных конечных элементах для кусочно-линейных и кусочно-квадратичных базисных функций, а также на тетраэдральных конечных элементах для кусочно-линейных базисных функций. Для полученных в работе вариационных постановок представлены портреты генерируемых матриц СЛАУ.
Основные результаты, выносимые на защиту:
Разработана технология реализации противопотоковых схем на неструктурированных сетках с расчетом неизвестных в узлах конечных элементов и серединах их ребер для моделирования течений с преобладанием конвективного переноса. В контексте метода конечных элементов предложен новый способ аппроксимации потоков через границы ячеек двойственной сетки.
Разработана и программно реализована вычислительная схема на основе смешанных постановок с введением функции штрафа, позволяющая единообразно моделировать стационарные внутренние течения вязкой
несжимаемой жидкости (уравнения Навье - Стокса, Стокса) и идеальной несжимаемой жидкости (уравнения Эйлера). Предложен способ задания нормальной компоненты вектора скорости и на границе области в вариационных постановках уравнений Эйлера (и п|хеГ = h, n - вектор нормали к границе Г, h - заданная функция). Разработана структура данных для хранения элементов матрицы системы линейных алгебраических уравнений с расчетом неизвестных компонент вектора скорости и давления в гипервекторе.
Реализованы стабилизированные схемы Галеркина для моделирования трехмерных течений на тетраэдральном разбиении, в которых используется метод стабилизации давления и в качестве стабилизирующего оператора выбирается конвективная часть исходного оператора.
Показана зависимость интенсивности вторичного течения от величины сдвига скорости на входе. Определено влияние соотношения узкого и широкого сечений на погрешность вычислений при измельчении сетки в трубопроводах, состоящих из различных конструктивных включений (элементов сужения и расширения).
Обоснованность и достоверность результатов, полученных в диссертационной работе, обеспечивается корректностью постановок рассматриваемых задач и методов их решения, основывается на расчетах широко известных и рекомендуемых тестовых задач и сопоставлении результатов численных расчетов с результатами, полученными другими авторами.
Апробация работы. Созданный комплекс программ внедрен в расчетную практику в Сибирском центре научно-технического обеспечения АНО "Пром-безопасность - Сибирь". Результаты решения задач с доминированием конвекции применяются в учебном процессе на факультете летательных аппаратов Новосибирского государственного технического университета в курсе лекций "Экологические проблемы энергетики "и "Численные методы в зада-
чах экологии".
Основные результаты диссертации докладывались на: международной конференции молодых ученых по математическому моделированию и информационным технологиям, октябрь 2002 г., г. Новосибирск; региональной научной конференции студентов, аспирантов и молодых ученых „Наука. Техника. Инновации" (НТИ-2002), декабрь 2002 г., г. Новосибирск; международной конференции „Вычислительные и информационные технологии в науке, технике и образовании" (ВИТ-2003), сентябрь 2003 г., г. Усть-Каменогорск, Казахстан; IV всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям, ноябрь 2003 г., г. Красноярск; всероссийской научной конференции студентов, аспирантов и молодых ученых „Наука. Технологии. Инновации" (НТИ-2003), декабрь 2003 г., г. Новосибирск; международной конференции по вычислительной математике (МКВМ-2004), июнь 2004 г., Академгородок, г. Новосибирск; V всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям с участием иностранных ученых, ноябрь 2004 г., г. Новосибирск; всероссийской научной конференции студентов, аспирантов и молодых ученых „Наука. Технологии. Инновации" (НТИ-2004), декабрь 2004 г., г. Новосибирск; семинарах Новосибирского государственного технического университета, Института вычислительных технологий СО РАН, Института вычислительной математики и математической геофизики СО РАН.
Результаты выполненных исследований опубликованы в работах [1] - [12].
Автор выражает искреннюю признательность и глубокую благодарность научному руководителю д.т.н., профессору Элле Петровне Шуриной и к.ф.-м.н., с.н.с. ИВТ СО РАН Нине Юрьевне Шокиной за помощь и поддержку при работе над диссертацией.
Математическая модель движения идеальной жидкости: уравнения Эйлера
Трубопроводы, используемые в промышленности, представляют собой систему, состоящую из трубы („русла") и различных конструктивных включений (заслонки, резкие сужения, расширения и другое). Проведена серия расчетов для каналов с расширением и сужением русла. Показано влияние соотношения узкого и широкого сечений на погрешность вычислений при измельчении сетки.
В приложениях 1-3 приведены компоненты локальных матриц и векторов правых частей, вычисленные на треугольных конечных элементах для кусочно-линейных и кусочно-квадратичных базисных функций, а также на тетраэдральных конечных элементах для кусочно-линейных базисных функций. Для полученных в работе вариационных постановок представлены портреты генерируемых матриц СЛАУ.
Основные результаты, выносимые на защиту: 1. Разработана технология реализации противопотоковых схем на неструктурированных сетках с расчетом неизвестных в узлах конечных элементов и серединах их ребер для моделирования течений с преобладанием конвективного переноса. В контексте метода конечных элементов предложен новый способ аппроксимации потоков через границы ячеек двойственной сетки. 2. Разработана и программно реализована вычислительная схема на основе смешанных постановок с введением функции штрафа, позволяющая единообразно моделировать стационарные внутренние течения вязкой несжимаемой жидкости (уравнения Навье - Стокса, Стокса) и идеальной несжимаемой жидкости (уравнения Эйлера). Предложен способ задания нормальной компоненты вектора скорости и на границе области в вариационных постановках уравнений Эйлера (и пхеГ = h, n - вектор нормали к границе Г, h - заданная функция). Разработана структура данных для хранения элементов матрицы системы линейных алгебраических уравнений с расчетом неизвестных компонент вектора скорости и давления в гипервекторе. 3. Реализованы стабилизированные схемы Галеркина для моделирования трехмерных течений на тетраэдральном разбиении, в которых используется метод стабилизации давления и в качестве стабилизирующего оператора выбирается конвективная часть исходного оператора. 4. Показана зависимость интенсивности вторичного течения от величины сдвига скорости на входе. Определено влияние соотношения узкого и широкого сечений на погрешность вычислений при измельчении сетки в трубопроводах, состоящих из различных конструктивных включений (элементов сужения и расширения).
Обоснованность и достоверность результатов, полученных в диссертационной работе, обеспечивается корректностью постановок рассматриваемых задач и методов их решения, основывается на расчетах широко известных и рекомендуемых тестовых задач и сопоставлении результатов численных расчетов с результатами, полученными другими авторами.
Апробация работы. Созданный комплекс программ внедрен в расчетную практику в Сибирском центре научно-технического обеспечения АНО "Пром-безопасность - Сибирь". Результаты решения задач с доминированием конвекции применяются в учебном процессе на факультете летательных аппаратов Новосибирского государственного технического университета в курсе лекций "Экологические проблемы энергетики "и "Численные методы в задачах экологии".
Основные результаты диссертации докладывались на: международной конференции молодых ученых по математическому моделированию и информационным технологиям, октябрь 2002 г., г. Новосибирск; региональной научной конференции студентов, аспирантов и молодых ученых „Наука. Техника. Инновации" (НТИ-2002), декабрь 2002 г., г. Новосибирск; международной конференции „Вычислительные и информационные технологии в науке, технике и образовании" (ВИТ-2003), сентябрь 2003 г., г. Усть-Каменогорск, Казахстан; IV всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям, ноябрь 2003 г., г. Красноярск; всероссийской научной конференции студентов, аспирантов и молодых ученых „Наука. Технологии. Инновации" (НТИ-2003), декабрь 2003 г., г. Новосибирск; международной конференции по вычислительной математике (МКВМ-2004), июнь 2004 г., Академгородок, г. Новосибирск; V всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям с участием иностранных ученых, ноябрь 2004 г., г. Новосибирск; всероссийской научной конференции студентов, аспирантов и молодых ученых „Наука. Технологии. Инновации" (НТИ-2004), декабрь 2004 г., г. Новосибирск; семинарах Новосибирского государственного технического университета, Института вычислительных технологий СО РАН, Института вычислительной математики и математической геофизики СО РАН.
Дискретные аналоги вариационных постановок
Построим двойственные сетки, которые будут использоваться при аппроксимации конвективных членов противопотоковыми схемами [80, 103]. Подобласти, которые образуются при построении двойственной сетки, являются ячейками двойственной сетки (аналогами конечных объемов). Рассмотрим два типа ячеек для аппроксимации конвективных членов противопотоковыми схемами: 1) ячейки, построенные вокруг вершин симплексов первичной сетки; 2) ячейки, построенные вокруг середин сторон симплексов первичной сетки.
Ячейки, построенные вокруг середин сторон симплексов, используются в противопотоковых схемах с расчетом неизвестных в серединах ребер симплексов. Вершинами ячейки fJj, построенной вокруг середины стороны г симплексов первичной сетки, будут барицентры и вершины треугольников, а ребрами - отрезки медиан треугольников, соединяющие вершины с барицентрами (рис. 2, а). Если сторона треугольника і принадлежит границе расчетной области, то ей соответствует ,;неполная" ячейка П; (рис. 2, б).
Первичная и двойственная сетки: а) ячейка S\, построенная вокруг середины j не принадлежащего границе расчетной области; б) "неполная" ячейка Q;, построенная вокруг середины р принадлежащего границе расчетной области
При равенстве выражения в левой части (2.12) нулю может не существовать uh и ph, удовлетворяющих (2.4), (2.6) или (2.8). Если inf-sup-выражение из (2.12) положительно, но стремится к 0 при h —» 0, то единственное решение хотя и существует, но может не сходиться к непрерывному решению дифференциальной задачи. О том, как проверять выполнение условия ЛББ-устойчивости можно найти в [15, 64]. При выборе интерполяционных функций из пространств, не удовлетворяющих условию ЛББ, в поле давления появляется шахматный эффект. Например, аппроксимация с использованием 12 степеней свободы для компонент вектора скорости и трех для давления (в узлах треугольника) может проявлять осцилляции поля давления. Важным следствием условия ЛББ является то, что скорость сходимости для давления и компонент вектора скорости можно точно оценить в пространстве Vg, тогда как, если условие ЛББ не выполняется, скорость сходимости можно вычислить лишь в пространстве дивергентно-свободных функций. Несмотря на то, что ЛББ гарантирует оптимальную сходимость конечно-элементного решения, это только достаточное, но не необходимое условие, то есть сходимость может быть достигнута в некоторых случаях без выполнения условия ЛББ.
Численные результаты показывают [86], что элементы Лагранжа (девять расчетных узлов) дают наибольшую точность при вычислении компонент вектора скорости и давления, элементы серендипова семейства (восемь расчетных узлов) показывают наименьшую точность при расчете давления. В [86] отмечено, что на треугольных элементах (шесть расчетных узлов) результаты расчетов зависят от шаблона триангуляции и на прямоугольных элементах (четыре расчетных узла) генерируются ошибочные решения для давления, если не используются специального типа граничные условия.
Рассмотрим часто используемые пространства интерполяционных функций для компонент вектора скорости и давления, удовлетворяющих условию ЛББ [15, 57]: Pf-Pl (пространство MINI-элементов); р +1_р ; к 0 (пространство Тейлора - Худа); пространство Равьяра - Томаса; пространство Крузея - Равьяра. Применение пространства Равьяра - Томаса для аппроксимации компонент вектора скорости сохраняет непрерывной нормальную компоненту скорости, давление исключается из системы уравнений, что приводит к положительно определенной дискретной системе со значительно меньшим числом неизвестных [56, 102].
Идея MINI-элементов заключается в пополнении пространства кусочно-непрерывных функций, аппроксимирующих компоненты вектора скорости, функциями-шапочками (bubble function) [58, 75, 76, 106], которые равны нулю на границе каждого элемента разбиения области. Такой способ построения решения включает процедуру исключения внутренних переменных и применяется в стабилизированных методах конечных элементов, таких, как метод RFB (от англ. residual-free bubbles) [58].
Стабилизированные методы конечных элементов
Представим функции uh Є Vft, ph Є Ph в виде разложения по базису (2.1)-(2.2). Подставив (2.1)-(2.2) в стабилизированные дискретные вариационные постановки (2.9)-(2.11) уравнений Навье - Стокса, Стокса и Эйлера, перейдем к вычислению элементов локальных матриц.
Определим локальные базисные функции ф{ компонент вектора скорости и давления, равные единице в расчетном узле і и нулю в остальных узлах конечного элемента. Получим вклады конечного элемента ек в глобальную матрицу СЛАУ.
Введем следующие обозначения для матриц, входящих в глобальную конечно-элементную матрицу СЛАУ, и глобального вектора правой части: D - диффузионная матрица, составленная из локальных матриц dk; С -конвективная матрица, сгенерированная из локальных матриц Ck; Q - матрица, соответствующая дискретному дивергентному ограничению, получена из локальных матриц qk; М - матрица массы, определяемая локальными матрицами т ; S, К - матрицы стабилизации, составленные из локальных матриц Sk, kk; f - вектор правой части, состоящий из вкладов локальных векторов fk. Матрицы с нижним индексом SUPG (SSUPG, SUPG, SUPG) получены поэлементной сборкой соответствующих локальных матриц Sjfc, kfc, Cfc, элементы которых умножаются на параметр стабилизации TSUPG [120]. Матрицы с нижним индексом PSPG (DPSPG, KPSPG) генерируются из соответствующих локальных матриц djt, kfc, элементы которых умножаются на параметр стабилизации TPSPG [120].
Введем на конечном элементе еп локальную нумерацию вершин, определим части ячеек О; = 17( П еп ( - ячейка, построенная вокруг вершины г симплекса) и границу между частями ячеек Sk = dftiPidUj, где г, j, к = 1,2,3 и j ф г ф к (рис. 4). Конечный элемент с узлами г, j", к содержит части трех ячеек fij, lj, Q,k- Вклад треугольного конечного элемента с глобальными номерами г, j, к расчетных узлов в глобальную матрицу представлен в виде локальной матрицы 3x3.
Для реализации противопотоковой схемы с вычислением неизвестных в серединах ребер симплексов рассмотрим конечный элемент еп, содержащий части трех ячеек Oj = Пеп и границу между частями ячеек Sk = d&iddUj, Qi - ячейка, построенная вокруг середин ребер МІ симплекса, г, j, к = 1,2,3, j Ф і ф к (рис. 6). Элементы Cjk локальной конвективной матрицы" конечного элемента еп имеют вид (2.25), где для противопотоковой схемы с расчетом неизвестных в серединах ребер симплекса введены следующие обозначения: Sk = 5 П dty - граница между частями ячеек конечного элемента, образованная медианой, опущенной в середину стороны с локальным номером Mfc; iijk единичная нормаль, внешняя по отношению к ячейке с барицентром j и восстановленная к границе Sk, k,j = 1,2,3. Для нормалей, входящих в выражения (2.25) имеют место соотношения: пгі = — Пзь піз = — Пгз, П32 = — Пі2 (рис. 6).
Элементы локальной конвективной матрицы конечного элемента имеют вид (2.30), где ті - выражения из (2.33), полученные интегрированием базисных функций по границам ячеек, построенных вокруг середин ребер симплексов.
Численное решение уравнений конвективно-диффузионного типа в двумерной области. Проти-вопотоковые схемы
Исследование интерполяционных свойств противопотоковых схем учета конвективных членов проведено на уравнениях конвективно-диффузионного типа. Решение задачи (3.2) выполнено следующими вычислительными схемами: конечно-элементной противопотоковой схемой с расчетом неизвестных в вершинах симплексов (п. 2.4.5); стандартным методом конечных элементов; конечно-элементной противопотоковой схемой с расчетом неизвестных в серединах сторон симплексов (п. 2.4.6); методом Петрова - Галеркина со стабилизацией оператора конвекции в направлении линии тока (SUPG); методом конечных элементов/конечных объемов (МКЭ/МКО), в котором диффузионные члены уравнения аппроксимируются с использованием метода конечных элементов, а конвективные - метода конечных объемов. Для аппроксимации скалярной функции используются линейные базисные функции.
Расчеты проведены для значений VT = 1, Ю-1, Ю-2, 10 3. Равномерное разбиение расчетной области выполнено сетками с левоориентированным шаблоном. Сетка размерностью 41 х41 включает 3200 треугольников, размерность генерируемой на этой сетке матрицы СЛАУ равна 1681 и заполненность матрицы составляет 4880 ненулевых элементов. Сетка размерностью 81 х 81 состоит из 12800 треугольников, размерность матрицы СЛАУ, генерируемой на этой сетке, равна 6561, заполненность матрицы - 19360 ненулевых элементов. В табл. 1-5 приведены значения гт = \\Т — Th\\/\\T\\ для задачи конвекции-диффузии (3.2). При сопоставлении полученных значений гт из табл. 1-5 видно, что для всех рассмотренных методов наблюдается сходимость на измельчающихся сетках.
Уменьшение значения VT приводит к уменьшению значения гт в методе SUPG, для остальных методов наблюдается увеличение значений гт- При доминировании конвективных членов классический МКЭ дает наихудший результат.
В табл. 7-8 представлены результаты решения краевой задачи методами SUPG и МКЭ на линейных базисных функциях при различных значениях VT на двух пространственных сетках с шагом по времени At = 1.25--10-1 и At = 6.25 Ю-2, где nmax - максимальное количество итераций метода решения СЛАУ. При доминировании конвекции (VT ОКОЛО Ю-1, 10 2) метод SUPG показывает лучшие результаты на пространственных и временных сетках по сравнению с МКЭ.
Проведен ряд тестовых расчетов на функциях пространства Тейлора Худа, вид элементов локальных матриц представлен в приложениях 1, 2. На рис. 10 отображено влияние значений штрафного параметра на значение евклидовой нормы \\р — ph\\ для алгоритма Удзавы и смешанного метода конечных элементов. Изменение параметра в большей мере оказывает влияние на точность вычисления давления, чем на точность вычисления компонент вектора скорости. Значения штрафного параметра, при которых ошибка вычисления давления наименьшая, лежат в пределах от 2.5 до 10 для алгоритма Удзавы и от 0.5 до 2.5 для смешанного МКЭ.
Для этой же задачи с известным аналитическим решением (3.4) исследована точность решения уравнений Навье - Стокса стабилизированным МКЭ. В стабилизированном методе конечных элементов параметр стабилизации определяется на каждом конечном элементе по формуле [120].
С измельчением сетки увеличивается число итераций решения СЛАУ (табл. 11), вычисляемое значение параметра стабилизации уменьшается и примерно в 1.8 раза уменьшаются ошибки при вычислении компонент вектора скорости и давления.
В начале поворотного участка проекция вектора скорости направлена к внутренней стенке канала, в конце этого участка - к внешней стенки канала. Вблизи плоскости симметрии вектор скорости перпендикулярен поперечному сечению.
В табл. 21 приведены значения дискретной дивергенции (3.1). Для решения задачи использован стабилизированный метод SUPG/PSPG, в котором параметр стабилизации вычисляем по формуле (3.5). При расчете стабилизированным методом конечных элементов значения дискретной дивергенции меньше по сравнению со значениями, полученными смешанным МКЭ.
Выполнено сравнение полученных результатов с опубликованными данными других авторов. На рис. 19 представлены расчеты (значения компоненты вектора скорости иу в правом поворотном изгибе канала), полученные стабилизированным МКЭ, смешанным МКЭ и МКР, реализованном в переменных векторный потенциал - вектор вихря [43, 47, 95]. Значения компоненты вектора скорости, вычисленные этими методами, отличаются друг от друга на величину 10_3.
