Содержание к диссертации
Введение 4
Глава 1. Итерационное решение систем линейных алгебраических уравнений 18
§ 1 Постановка задачи 18
§2 Методы решения систем линейных алгебраических уравнений с незнакоопределенной почти особенной матрицей 28
2.1 Итерационная схема неполной аппроксимации 28
2.2 Асимптотическое свойство схемы 30
2.3 Ускорение сходимости 33
2.4 Многопараметрическая оптимизация 34
§3 Движение идеальной несжимаемой стратифицированной жидкости 41
3.1 Уравнение Гельмгольца 41
3.2 Разностная схема 41
3.3 Тестовые расчеты 43
Глава 2. Итерационное решение систем нелинейных алгебраических уравнений 54
§1 Постановка задачи 54
§2 Итерационная градиентная схема решения систем нелинейных алгебраических уравнений 56
2.1 Явная итерационная схема 56
2.2 Многопараметрическая оптимизация 57
2.3 Ускорение сходимости , 62
§3 Модельное уравнение переноса 65
3.1 Уравнение Бюргерса 3.2 Разностная схема 65
3.3 Сравнительные расчеты 67
Глава 3. Решение систем нелинейных дифференциальных уравнений Навье-Стокса 69
§1 Двухмерная задача движения вязкой несжимаемой жидкости в формулировке «функция тока»-«вихрь» 69
1.1 Постановка задачи 69
1.2 Разностная схема 70
1.3 Тестовые расчеты 75
§2 Двухмерная задача движения вязкой несжимаемой жидкости в формулировке «скорость»-«давление» 86
2.1 Постановка задачи 86
2.2 Разностная схема 87
2.3 Тестовые расчеты 89
§3 Трехмерная задача движения вязкой несжимаемой жидкости в фор
мулировке «скоро сть»-«давление» 94
3.1 Постановка задачи 94
3.2 Разностная схема 95
3.3 Тестовые расчеты Заключение 124
Введение к работе
В соответствии с изречением греческого философа Гераклита, который говорил «все течет...», в повседневной жизни мы постоянно сталкиваемся с течением обычных жидкостей: воды, воздуха, крови и т.д., в очень простых ситуациях таких как дыхание, мытье рук и т.п.
Большинство течений имеет природное (океаны, моря, ветер) и техногенное происхождение (самолеты, машины, биоинженерия). Существует потребность в моделировании проблем течения жидкостей, с целью лучшего понимания сложных явлений и повышения качества технологий. С появлением новых вычислительных инструментов моделирование становится все более подходящим для проведения экспериментов. При некоторых обстоятельствах, прямой эксперимент может быть слишком дорог и даже привести к разрушению исследуемого объекта, таким образом, моделирование является единственным способом изучить параметры исследуемой системы.
Появление компьютеров привело к развитию вычислительной гидродинамики. В начале, новая дисциплина принесла надежду, что вычисления могут открыть дорогу к моделированию трехмерных течений с высоким числом Рей-ноль дса. Реальность разрушила эти наивные ожидания. До сих пор остаются актуальными многие проблемы, ограничивающие использование существующих численных методов решения.
Большинство течений, с которыми человеку приходится сталкиваться, имеют нелинейную природу. Таковыми, например, являются течения вязкой несжимаемой жидкости, движение которой описывается нелинейными дифференциальными уравнениями Навье-Стокса. Используя различные упрощения моделей гидродинамики можно нелинейную задачу свести к линейной. Например, течение идеальной несжимаемой стратифицированной жидкости в приближении Буссинеска описывается одним линейным дифференциальным уравнением Гельмгольца относительно функции тока. Упрощенные модели движения жидкости отвечают многим процессам, происходящим в природе, поэтому их исследование также имеет огромное значение.
Как правило, при решении задач гидродинамики используют конечно-разностные аппроксимации исходных дифференциальных уравнений. Если аппроксимируется линейная задача, то получающиеся равенства являются системой линейных алгебраических уравнений вида
An = f, (1)
где A = ipLij) — матрица размерности m, f, и — неизвестный и известный m-мерные векторы из конечномерного гильбертового пространства Нт.
Первый итерационный метод для решения совместной системы линейных алгебраических уравнений (1) был предложен Гауссом [161] в 1828 году. Таким образом, история развития итерационных методов решения (1) насчитывает уже 175 лет. Поскольку работ посвященных данной теме достаточно много [35], остановимся лишь на ключевых моментах.
Большинство итерационных методов решения (1) можно разделить на два типа. Первый — методы, использующие информацию о спектральных свойствах оператора А. Операторы перехода от одного начального приближения к другому в алгоритмах этого типа зависят от некоторых параметров, которые выбираются из условия наиболее быстрой сходимости. Для применения итерационного процесса, построенного исходя из данного принципа, нужно располагать возможно большей информацией о матрице системы, в частности о границах ее спектра. Более подробно об этих методах см. [7,27,49,100,127,153,154,194,195,201,232,258-260,264]. Их недостатком является трудоемкость получения информации о спектре матрицы. В некоторых случаях, например, при решении уравнения Гельмгольца движения идеальной стратифицированной жидкости в сложной области с негладкой границей или с условиями на бесконечности, расчет границ спектра представляет собой не менее сложную задачу чем решение исходной системы уравнений.
Методы второго типа — это вариационные или градиентные методы, в основе которых лежит принцип минимизации некоторого функционала, минимум которого достигается на искомом решении системы. Подробнее см. [121,125,127,133,150-152,154,155,161].
При решении (1) вариационными методами большую роль играют такие свойства оператора А как знакоопределенность, самосопряженность, число обусловленности.
Если матрица А незнакоопределена, возникает необходимость преобразования системы (1). Умножая (1) на А (1-ая трансформация Гаусса), мы по-лучим новую систему с положительно определенной матрицей А = А А и правой частью Л і , которую можно решать обычными методами. Но когда detA « 0, такой путь неприемлем. Покажем это на примере. Пусть А = А , тогда получим систему А2и = АЇ.
Собственные числа д; матрицы А2 будут равны X2. Если какое-либо собственное число Afc w 0, например, пусть А = 10 5, то соответствующее собственное число новой матрицы ц будет равно 10 10, а это означает, что придется решать СЛАУ с фактически особенной матрицей. Таким образом, использование трансформации Гаусса в случае незнакоопределеннои, почти особенной матрицы приводит к существенному изменению свойств системы и полученного приближенного решения, т.к. в случае особенной матрицы решение неединственно и итерационный процесс может сойтись к любому из решений с точностью до вектора из ядра оператора.
Итерационные методы решения СЛАУ с несимметричной матрицей рассматривались в следующих работах [8,21,38,135,156,160,174,196,224,247]. В последнее время большое распространение получил метод обобщенных минимальных невязок (GMRES) решения систем линейных алгебраических уравнений с незнакоопределеннои матрицей [244,246,255,262]. Однако, сходимость
этого метода не доказана и его реализация предъявляет большие требования к объему памяти и количеству операций, т.к. смысл этого алгоритма заключается в ортогонализации большого набора векторов. В связи с этими ограничениями, как правило, проводят перезапуск метода после ортогонализации заданного числа векторов, которое меньше размерности системы. Однако, исследования показывают, что такое упрощение может привести к отсутствию какой-либо сходимости.
Скорость сходимости итерационных методов зависит от спектральных свойств матрицы коэффициентов. Таким образом, можно сделать попытку трансформировать исходную систему уравнений в эквивалентную ей, но с матрицей коэффициентов, имеющей лучший спектр, Предобуславливатель — матрица, проводящая данное преобразование. Например, если матрица М аппроксимирует матрицу А в некоторым смысле, система вида
М 1Ли = М Н
имеет то же решение, что и исходная, но спектральные свойства матрицы М 1А могут быть лучше.
Использование предобуславливателя в итерационном методе оборачивается дополнительными вычислительными затратами при его построении и использовании на каждой итерации.
Работы [109,175,177,182,189,197,199,245,263] посвящены предобуслав-ливанию систем линейных алгебраических уравнений (СЛАУ).
При решении СЛАУ хорошо зарекомендовали себя схемы неполной аппроксимации (НА), которые не являются инвариантными относительно решения (1). Эти схемы, впервые были предложенные Н.Н. Яненко, позволяют решать СЛАУ с незнакоопределенной матрицей без использования 1-й трансформации Гаусса. Существует теорема, доказывающая их сходимость. В ходе итерационного процесса, проводимого схемами НА, возможно применение ускорения сходимости, которое существенно сокращает время расчета, даже в
незнакоопределенном случае.
В работах [58,61,64,68,69,71,72,75,78,84,88,89,94-96] схемы неполной аппроксимации используются для решения линейных задач гидродинамики.
Необходимость решения СЛАУ возникает, например, при аппроксимации линейных дифференциальных уравнений гидродинамики (уравнение Гельм-гольца движения идеальной стратифицированной жидкости), однако, линейными уравнениями описываются значительно упрощенные модели течений жидкости. Например, отказ от допущения, что жидкость является идеальной, приводит к необходимости решения уже нелинейных дифференциальных уравнений Навье-Стокса, описывающих движение вязкой жидкости.
В двухмерном случае при численном исследовании течений вязкой несжимаемой жидкости используются уравнения Навье-Стокса, записанные как в переменных скорость-давление (v, Р) так и в переменных функция тока-вихрь (ф, со) (см. обзор [202]).
В последнем случае при численном расчете в первую очередь встает проблема постановки краевых условий для вихря на твердых стенках, которые не присутствуют явно в постановке задачи. Условия первого, второго порядка, получающиеся путем разложения условия прилипания дф/дп вблизи границы [147], с успехом использовались в численных расчетах многими авторами (см., например, [23,29,204]). В некоторых работах отмечается, что использование условий более высокого порядка аппроксимации при расчете методами стационирования может привести к зарождению численной неустойчивости. Кроме того, эти условия при решении явными схемами существенно замедляют сходимость, а при решении неявными схемами требуют организации отдельного итерационного процесса.
Система Навье-Стокса в физических переменных скорость-давление до недавнего времени обладала меньшей популярностью вследствие того, что в ней отсутствует уравнение для давления и краевые условия для него.
Обычно недостающее уравнение для Р получают либо путем дифферен цирования остальных уравнений системы с учетом соленоидальности поля скоростей, либо путем введения в уравнение неразрывности некоторой функции давления по времени. Граничные условия для полученного уравнения обычно являются следствием теории пограничного слоя и имеют, как правило, градиентный вид. Многие авторы, решающие систему Навье-Стокса в физических переменных, успешно используют условия градиентного типа для давления [144,229], однако, встречаются работы [235,250], в которых показано, что градиентные условия для давления при расчете нестационарной задачи являются некорректными.
Трехмерные задачи движения вязкой жидкости в большинстве своем решаются в переменных скорость-давление, однако, существуют некоторые работы, исследующие систему в переменных векторный потенциал-вихрь. Проблемы задания краевых условий для давления или векторного потенциала сохраняются и в пространственных задачах.
И двухмерная и трехмерная задача может включать в себя условия на бесконечности. В этом случае появляется проблема их переноса с бесконечности на границу конечной области. Использование на неотражающих границах каких-либо сложных равенств затрудняет использование существующих схем решения, замедляет сходимость и может привести к неустойчивости. Краевые условия на бесконечности а также неотражающие условия рассматривались в работах [42,48,99,106,142,143,145,149,176,203,217,228,231,251]. Отметим работу [98], содержащую обзор по неотражающим условиям на границах области расчета.
Краевые условия для уравнений Навье-Стокса исследуются также в работах [25,26,36,56,57,130-132,148,158,227,233,243,254].
В последнее время большое распространение получили методы модификации исходной системы уравнений путем введения искусственной вязкости [40], псевдо- [192,223,225,226,240] или слабосжимаемости (s-аппроксимация) [44,52,102,103], впервые предложенные Н.Н. Яненко.
В существующих итерационных схемах решения стационарных уравнений Навье-Стокса решение, в основном, ищется путем стационирования разностных уравнений, аппроксимирующих нестационарную систему Навье-Стокса. Преимуществом явных схем является экономичность, однако их условия устойчивости очень ограничены. Таким образом, шаг по времени в некоторых задачах должен быть достаточно мал. Неявные схемы, наоборот, имеют более широкую область устойчивости, требуя взамен большего числа операций для обращения блочно-диагональных матриц. Еще одним недостатком неявных схем является бесконечная скорость распространения ошибки по области решения.
В неявных схемах можно избавиться от нелинейности, если коэффициенты матрицы оператора, зависящего от решения, рассчитывать используя приближенное решение, полученное на предыдущем временном слое.
Наиболее эффективными неявными методами являются методы переменных направлений [3,4,180,210], расщепления (в иностранной литературе алгоритм SIMPLE и его модификации) [33,33,46,103,115,120,159,183,184,186, 191,206,210,211,257].
Лишь в очень немногих работах нелинейные уравнения, аппроксимирующие стационарную систему Навье-Стокса, решаются как система нелинейных алгебраических уравнений [51,179,208,213,242]. Это связано с тем фактом, что методы решения нелинейных уравнений накладывают сильные ограничения на оператор и начальное приближение. Например, в работе [178] отмечена большая чувствительность метода Ньютона к начальному приближению. Решение нелинейных операторных уравнений рассматривалось также в [11,14,20,28,37, ПО, 111,205,219,220].
При решении задач гидродинамики достаточно эффективными оказались многосеточные методы (см. обзор [253]). Принцип их работы может быть объяснен с использованием Фурье-анализа вектора ошибки. Оказывается, при получении следующего приближенного решения с помощью какого-либо итерационного метода уменьшаются те Фурье-составляющие вектора ошибки, длинна волны которых сравнима с размером шага сетки. Таким образом, чередование точных и грубых сеток позволяет уменьшать все Фурье-составляющие вектора ошибки. Многосеточные методы также применялись в работах [141,162,181,185,193,207,212,221,222,237-239,248,249,261,265].
Для решения систем нелинейных алгебраических уравнений также применимы схемы неполной аппроксимации, которые оказались достаточно эффективными при решении широкого класса задач. Применение схем НА позволяет использовать сложные краевые условия при решении задач гидродинамики. При решении стационарных задач течения вязкой жидкости схемами НА отпадает необходимость проведения внутренних итераций для решения промежуточных линейных задач.
Схемы неполной аппроксимации на примере различных нелинейных задач гидродинамики рассматривались в работах [12,13,50,59,60,62,63,65-67,73, 74,76,79-81,83,85-87,91-93,97,128,136-138,146,171,173,215,216].
Для получения более полной информации о численных методах решения уравнений Навье-Стокса см. также обзоры [17,147,170,202] и работы [1,2,9, 15,24,39,41,43,45,47,52,53,55,101,112,116,117,124,126,169,187,188,234,256].
Из последних работ по численному исследованию уравнений Навье-Стокса отметим следующие [6,10,22,30,31,34,54,104,107,108, ИЗ, 114,118, 122,123,139,140,157,163-167,172].
В диссертации исследуются схемы неполной аппроксимации решения как линейных, так и нелинейных задач. Строятся алгоритмы, позволяющие проводить в многопараметрических схемах неполной аппроксимации оптимизацию по нескольким параметрам одновременно, а также алгоритмы существенно ускоряющие сходимость процесса решения нелинейных задач. Исследуются способы задания краевых условий для вихря. Строится метод задания недостающих условий для давления, в случае решения задачи в физических переменных. Проводятся численные эксперименты по расчету двухмерных и трех
мерных внутренних задач и задач протекания в различных областях.
Первая глава состоит из трех параграфов и посвящена построению итерационного алгоритма неполной аппроксимации с многокомпонентной групповой оптимизацией для решения систем линейных алгебраических уравнений (СЛАУ).
В параграфе 1 приводится постановка задачи решения системы линейных алгебраических уравнений, анализируются основные проблемы, связанные с решением этой задачи, дается обзор существующих методов решения.
Параграф 2 посвящен построению итерационного алгоритма неполной аппроксимации с многокомпонентной групповой оптимизацией. Приводится описание итерационной схемы неполной аппроксимации для решения системы линейных алгебраических уравнений и процедуры ускорения сходимости итерационного процесса [90].
Рассмотрим итерационную градиентную схему НА
un+i/2 = цп _ Тп+1в(Лип - Л, (2а)
un+i = un+i/2 + Qn+lZ«) п = 0,1,2,..., (2Ь)
где В — неособенная квадратная матрица, zn Є Hm — произвольный вектор, Тп+и &п+1 — итерационные параметры, и0 — произвольный начальный вектор. Параметры итерации тп+ь ct +i выбираются исходя из минимума нормы векторов vn+1 2 = D 2(un+1 2 — и) и v"+1 = Dl/2(un+1 — u) соответственно {D — самосопряженный положительноопределенный оператор). Заметим, что ап+\ может быть как константой, так и матрицей.
Доказано, что существуют такие вектора zn+1, что последовательность {0H+i — v"+1 /v"2} является монотонно возрастающей и ограниченной сверху числом 0 1. Таким образом, схема НА является сходящейся и обладает асимптотическим свойством. Благодаря этим свойствам, процедура ускорения сходимости схемы неполной аппроксимации позволяет существенно сократить время расчета.
Можно проводить оптимизацию по нескольким параметрам одновременно. Пусть ап+і — диагональная матрица размерности т и zn — вектор, все элементы которого равны 1. Если находить элементы матрицы {( } из условия min vn+1, это приведет к необходимости решения системы линейных алгебраических уравнений, матрица которой, в случае блочно-диагональной матрицы Л, также будет блочно-диагональной, что упрощает ее решение. Кроме того, в этом случае можно так проводить минимизацию, что для нахождения параметров итерации a„j нужно будет решать систему уравнений с диагональной матрицей.
Параграф 3 гл. 1 содержит результаты расчетов некоторых тестовых задач.
Для проверки эффективности схем неполной аппроксимации по сравнению с обычным градиентным методом были проведены расчеты задач движения идеальной несжимаемой равномерно стратифицированной жидкости, описываемого уравнением Гельмгольца.
С помощью рассмотренного итерационного метода были получены приближенные решения тестовых задач в областях с различными формами и размерами. В случае расчетов краевых задач с условиями на бесконечности увеличение размеров области расчета и уменьшение сеточного шага практически не оказывало влияния на получаемые картины течения.
Сравнительные расчеты задач с незнакоопределенной несимметричной матрицей показали превосходство схемы неполной аппроксимации над обычным градиентным методом.
Вторая глава состоит из трех параграфов и посвящена построению итерационного алгоритма многокомпонентной групповой оптимизации и алгоритма ускорения сходимости итерационного процесса решения системы билинейных уравнений.
В параграфе 1 излагается общая физическая постановка задачи движения вязкой несжимаемой жидкости, дается краткий обзор литературы и результатов, полученных ранее другими авторами. Приводится постановка задачи
решения системы билинейных алгебраических уравнений.
Движение вязкой несжимаемой жидкости описывается системой нелиней
ных дифференциальных уравнений Навье-Стокса. Конечно-разностная ап
проксимации уравнений Навье-Стокса является системой нелинейных алгеб т раических уравнении вида
Л(и, u) = f, (3)
где и, f — соответственно неизвестный и известный вектора из Ят, Л(u, v) = A\(u,v) + A$v, Ai — линейный оператор, Лі(и, u) — билинейный оператор, такой что
А ащЫ + o2u(2), М(1) + b2v 2)) =
= 0 , v 1 ) + aiMi(u(1\v(2 )+o2b1 1(u(2 v(1 ) + o2Mi(u(2 v(2 (4)
где u , v J Є Hm, a,-, b( — произвольные постоянные.
Существующие в настоящее время методы решения нелинейных уравнений, например, метод Ньютона, накладывают серьезные ограничения на опе 4 ратор и начальное приближение, которое должно быть достаточно близким к
точному решению [178].
Параграф 2 посвящен построению итерационного алгоритма многокомпонентной групповой оптимизации и алгоритма ускорения сходимости итерационного процесса решения системы билинейных уравнений. Приводится описание итерационного градиентного метода решения систем билинейных
щ алгебраических уравнений [129]
un+l/2 = ип __ Tn+i [А(ип un) _ f ] (5а)
un+i = un+i/2 + an+lZnj п = о, 1, 2, . . ., (5Ь)
где z" Є Hm — некоторый вектор, u° — произвольное начальное приближение из области определения оператора Л, rn+i = const, an+i — итерационный
параметр, который может быть как константой, так и матрицей.
Параметр rrt+i находится из условия минимума нормы г"+1/2 с помощью свойства билинейности оператора и формулы Кардано решения кубических уравнений. Если параметр an+i — константа, будем выбирать его из условия минимума jrn+1, также используя свойство билинейности оператора и формулы Кардано.
Выбирая так тп+1 и ап+и получим гп+1[2 г"+1/22 г»2.
Пусть an+i = {а +1} — диагональная матрица.
(г)
Можно находить параметры c +i из условия глобального минимума [[г71"1"1!]. В этом случае, используя некоторые предположения относительно оператора Лі и векторов z", задача минимизации сводится к задаче решения си (І) стемы линейных уравнении относительно „+!•
Для указанного итерационного алгоритма (5) построена процедура ускорения сходимости, которая, как показали тестовые расчеты, существенно сокращает время, необходимое для получения хорошего приближенного решения.
Параграф 3 гл. 2 посвящен численному решению нелинейного дифференциального уравнения Бюргерса.
Были проведены сравнительные расчеты обычной градиентной схемой (5а), схемой (5а)-(5Ъ) с последовательной оптимизацией, схемой (5а)-(5Ь) с групповой оптимизацией. Имеются расчеты, показывающие эффективность применения ускоряющей процедуры. Кроме того, использование схемы (5 а)-(5Ь) позволило получить приближенное решение задачи с условием на бесконечности, практически совпадающее с точным.
Третья глава состоит из трех параграфов и посвящена решению системы нелинейных дифференциальных уравнений Навье-Стокса. Рассматривалось как двухмерное, так и трехмерное течение.
Параграф 1 посвящен численному решению системы нелинейных уравнений Навье-Стокса для «функции тока» и «вихря», описывающей двухмерное течение вязкой несжимаемой жидкости- Рассматривались задачи течения жид
кости в квадратной каверне с движущейся верхней крышкой, течения жидкости над каверной, течения жидкости в канале с обратным уступом.
Использование предложенного итерационного градиентного метода с многопараметрической оптимизацией, а также метода ускорения сходимости позволило провести расчеты в широком диапазоне коэффициента вязкости. Кроме того, алгоритм расчета позволяет использовать сложные краевые условия, которые возникают, например, при решении задач с условиями на бесконечности.
Полученные приближенные решения хорошо согласуются с результатами, полученными другими авторами [204,241].
Параграф 2 гл. 3 посвящен численному решению двумерной системы нелинейных уравнений Навье-Стокса, записанной для компонент вектора скорости и давления.
Данная задача представляет интерес прежде всего потому, что в ходе ее решения возникают те же самые проблемы, что и при решении трехмерной системы уравнений Навье-Стокса. А именно, отсутствие в постановке задачи краевых условий для давления на твердой стенке, проблема переноса краевых условий с бесконечности на границу конечной области, проблема эффективного решения разностных схем, аппроксимирующих дифференциальную систему уравнений.
В качестве тестовых были рассмотрены задачи течения жидкости в квадратной каверне с движущейся верхней крышкой, течения жидкости в канале с обратным уступом.
На примере этих задач были продемонстрированы методика задания недостающих краевых условий для давления на твердых стенках, способ переноса краевых условий с бесконечности на границу конечной области, итерационный градиентный метод решения уравнений Навье-Стокса в примитивных переменных.
Результаты расчетов всех задач практически совпадают с результатами, по
лученным при решении системы уравнений, записанной относительно «функции тока» и «вихря», которая была рассмотрена в первом параграфе третьей главы.
Поскольку предложенные алгоритмы оказались достаточно эффективными, они были применены при решении трехмерной задачи движения вязкой несжимаемой жидкости, описываемой системой уравнений Навье-Стокса, рассмотренной в параграфе 3 гл. 3.
В качестве тестовых рассматривались задачи движения жидкости в кубической каверне с движущейся верхней крышкой, движения жидкости в канале с обратным уступом, движения кубического тела в канале с жидкостью, результаты расчетов которых хорошо согласуются с численными [16,105,119,134,168,190,209,214,218,236] и натурными [252] экспериментами других авторов.
Полученные в диссертации теоретические результаты и численные расчеты показывают, что использование многопараметрической оптимизации в совокупности с ускорением сходимости и многосеточной методикой значительно сокращает время расчета как линейных, так и нелинейных задач. При этом, применяемые методы расчета позволяют использовать достаточно сложные краевые условия, возникающие при решении задач с условиями на бесконечности. Использование при решении двухмерных и трехмерных задач течения вязкой жидкости в формулировке скорость-давление предложенных методов задания на твердых стенках краевых условий для давления позволяет отказаться от введения искусственных условий, отсутствующих в исходной постановке задачи.
