Содержание к диссертации
Введение
1 Аналитические решения для уравнений Стокса в цилиндрических и сферических координатах 13
1.1 Уравнения Стокса в цилиндрических координатах 13
1.2 Построение решения в осесимметричном случае 15
1.3 Уравнения Стокса в сферических координатах 22
1.4 Построение эталонных решений в сферических координатах 23
1.5 Анализ решений и метод тестирования численных алгоритмов 27
1.6 Выводы 30
2 Решение уравнений Стокса многосеточным методом в цилиндрической и сферической системах координат 31
2.1 Многосеточный метод в цилиндрической системе координат 31
2.1.1 Метод Гаусса-Зейделя для уравнений Стокса 34
2.1.2 Дискретизация уравнений Стокса и неразрывности 35
2.1.3 Проверка численной схемы на эталонных решениях 40
2.2 Многосеточный метод в сферической системе координат 45
2.2.1 Дискретизация уравнений для метода Гаусса-Зейделя 46
2.2.2 Проверка численной схемы на эталонных решениях 53
2.3 Выводы 61
3 Аналитические решения для плоских течений 62
3.1 Аналитические решения для Стоксова течения над полостью 62
3.1.1 Постановка задачи 62
3.1.2 Построение решений для течений над полостью 63
3.1.3 Картина линии тока для течений над полостью 70
3.2 Течение в полости, имеющей форму эллиптической подковы 71
3.2.1 Постановка задачи 71
3.2.2 Построение решений 72
3.2.3 Геометрическая картина течений
3.3 Проверка конечно-разностной схемы на эталонных решениях 78
3.4 Выводы 81
Заключение 82
Литература 83
- Уравнения Стокса в сферических координатах
- Анализ решений и метод тестирования численных алгоритмов
- Дискретизация уравнений Стокса и неразрывности
- Картина линии тока для течений над полостью
Введение к работе
Актуальность темы.
Математическое моделирование течений в вязких средах является
актуальной задачей для различных областей науки. В геодинамике при
построении моделей течения Земной мантии традиционно применяются система уравнений Стокса. В ряде исследований система уравнений Стокса дополняется уравнением переноса тепла и используется для описания процесса конвекции. Уравнения Стокса используются при моделировании течения крови, а также течений в наноструктурах.
Для некоторых задач геодинамики возникает необходимость находить решения уравнений Стокса с учетом неоднородной вязкости среды. Вязкость горных пород, как правило, зависит от их температуры. Если температура верхних и нижних слоев среды значительно отличается, то моделирование течений для таких систем приводит к численному решению уравнений при высоких контрастах вязкости.
Важным этапом разработки программ численного решения уравнений, является процесс тестирования. Результатом тестирования может стать оценка корректности работы алгоритма. Одним из способов тестирования является сравнение результата работы численного алгоритма с аналитическим решением, полученным для некоторой тестовой задачи. Аналитическое решение в таком случае выступает в качестве эталона. Тестовая задача может быть специализирована для проверки корректной работы программы при варьировании той или иной группы параметров, входящих в уравнения. Например, программа для моделирования процесса субдукции в геофизике, должна в первую очередь выдавать корректное решение при расчете течений в системах с высокими контрастами вязкости и наличием гравитационного поля.
Естественным шагом к построению решений уравнений Стокса в областях, граница которых имеет сферическую или осевую симметрию, часто становится переход из декартовой в сферическую или цилиндрическую систему координат. Ряд исследований был посвящен нахождению аналитических решений уравнений Стокса в декартовой системе координат, ставших основой для создания тестовых задач, используемых для проверки численных алгоритмов. Значительно меньшее количество работ было посвящено нахождению аналитических эталонных решений в криволинейных системах координат. В то же время, построение численного решения уравнений в этих системах зачастую оказывается более сложным процессом, чем в декартовой системе координат.
Таким образом, тема диссертационной работы отвечает современным проблемам рассматриваемой области науки и поэтому является актуальной.
Целью данной диссертационной работы является развитие подходов к тестированию численных алгоритмов решения уравнений Стокса. Первой задачей является вывод аналитических решений уравнений Стокса с переменными вязкостью и плотностью в цилиндрической и сферической системах координат, позволяющих создавать на их основе тестовые задачи для
оценки качества работы численных методов и комплексов программ. Вторая задача состоит в нахождении аналитических решений краевых задач для уравнений Стокса, задающих двумерные течения в областях с эллиптическими границами.
Научная новизна результатов, полученных в работе:
1. Впервые получены частные аналитические решения уравнений Стокса с
переменными вязкостью и плотностью в цилиндрической и сферической
системах координат, для таких параметров системы, при которых вязкость,
плотность и напряженность гравитационного поля являются функциями
координаты г : г/ = г/(г), р = p(r), G = G(r).
2. Предложены и протестированы численные схемы для многосеточного
метода решения уравнений Стокса с переменной вязкостью на ортогональных
смещенных сетках в цилиндрической и сферической системах координат.
-
Выведено явное решение краевых задач для уравнений Стокса, задающих двумерные течения над полостью.
-
Выведено явное решение краевых задач для уравнений Стокса, задающих течения в области, имеющей форму эллиптической подковы.
Теоретическая значимость результатов данной работы связана с расширением класса задач для тестирования численных методов. Важным для теории является построение частных аналитических решений уравнений Стокса с переменными вязкостью и плотностью в цилиндрических и сферических координатах, нахождение аналитических решений краевых задач для уравнений Стокса для двумерного течения над полостью и для области, имеющей форму эллиптической подковы, вывод дискретизации для уравнений Стокса с переменной вязкостью на ортогональных смещенных сетках в цилиндрической и сферической системах координат.
Практическая значимость результатов подтверждается возможностью применения найденных эталонных решений при создании тестовых задач, используемых для тестирования численных алгоритмов, использованием программного комплекса для моделирования геофизических течений с переменной вязкостью.
Методы исследования. В работе используются методы математической физики решения дифференциальных уравнений в частных производных, численные методы решения систем дифференциальных уравнений в частных производных и методы программной инженерии.
На защиту выносятся:
1. Метод тестирования численных алгоритмов на основе эталонных решений
уравнений Стокса, моделирующих течения при малых числах Рейнольдса с
переменной вязкостью и плотностью в цилиндрической и сферической
системах координат.
2. Эталонные решения краевых задач для уравнений Стокса с постоянной
вязкостью, моделирующих течения вязкой жидкости в областях с
эллиптическими границами.
Достоверность результатов обеспечивает надежность используемых теоретических методов и согласованность теоретических расчетов с численными экспериментами. Результаты прошли многократную апробацию на научных конференциях.
Апробация работы. Основные результаты работы докладывались на международных конференциях «Mathematical Challenge of Quantum Transport in Nanosystems» (Санкт-Петербург, 2014 г. и 2015 г.), III Всероссийском конгрессе молодых ученых (Санкт-Петербург, 2014 г.), IV Всероссийском конгрессе молодых ученых (Санкт-Петербург, 2015 г.).
Публикации. Основные результаты по теме диссертации изложены в 8 печатных работах, из них 3 в журналах из Перечня ВАК, 1 в журнале, индексируемом в Web of Science и Scopus.
Объем и структура работы. Диссертация состоит из оглавления, введения, трех глав, заключения, списка сокращений и условных обозначений, списка литературы и приложений. Полный объем диссертации составляет 105 страниц с 48 рисунками. Список литературы содержит 88 наименований.
Построенные в диссертационной работе аналитические решения уравнений
и разработанный комплекс программ внедрены и используются для
тестирования алгоритмов расчета стоксовых течений в компании
ООО «Корнинг СНГ».
Уравнения Стокса в сферических координатах
Для оценки качества работы численных алгоритмов может быть использован следующий метод, основанный на сравнении численного решения с аналитическим. Выберем некоторое частное решение из множества аналитических решений, построенных в разделах 1.2 и 1.4. На основе данного решения поставим тестовую краевую задачу. В соответствующих криволинейных координатах рассмотрим область в форме прямоугольного параллелепипеда. Для компонент скорости на границе данной области зададим значения, равные значениям выбранного аналитического решения. В самой области решим уравнения Стокса и неразрывности численно, используя тестируемый численный алгоритм. В силу единственности решения с данными граничными условиями, численное решение должно совпадать с аналитическим. Выполнив сравнение этих двух решений оценим качество численной схемы.
Заметим, что из-за наличия особенностей у эталонных решений возникает необходимость исключать ось г = 0 в цилиндрической системе координат и точку г = 0 в сферической системе координат из области, в которой рассматривается течение. 1.6 Выводы
В главе построены частные аналитические решения уравнений Стокса и неразрывности с переменной вязкостью и плотностью в цилиндрической и сферической системах координат для моделей геофизических течений в мантии Земли и нанотечений в тубулярных наноструктурах. Рассмотрен случай, когда вязкость и плотность являются функциями радиуса. В сферической системе координат решения найдены для вязкости, являющейся степенной функцией радиуса. Описан метод, позволяющий на основе данных решений составлять тестовые краевые задачи для оценки качества работы численных алгоритмов решения уравнений Стокса и неразрывности с переменной вязкостью и плотностью. Глава 2
Решение уравнений Стокса многосеточным методом в цилиндрической и сферической системах координат Во многих задачах, особенно геофизических, стоксовы течения с переменными параметрами (вязкость, плотность) рассматриваются в областях, имеющих осевую или сферическую симметрию границы. Численные схемы, основанные на декартовых координатах, в этой ситуации работают плохо. Естественно использовать уравнения в цилиндрических или сферических координатах и соответствующие разностные схемы. Наличие высокого контраста вязкости в моделируемых системах создает дополнительные вычислительные сложности [46-48]. В данной главе построены численные схемы для многосеточного метода решения уравнений Стокса и неразрывности в цилиндрических и сферических координатах. Производится тестирование данных схем на эталонных решениях, выведенных в главе 1.
Численная схема для решения уравнений Стокса с переменной вязкостью в декартовой системе координат описана в [3]. Построим численную схему для решения уравнений Стокса с переменной вязкостью в цилиндрической системе координат. Численную схему будем рассматривать применимо к многосеточному методу. Решение уравнений Стокса и неразрывности строится на последовательности иерархически вложенных сеток с разным разрешением (Рисунок 2.1). При переходе на каждый следующий уровень разрешение сетки увеличивается. Как правило, итерации цикла построения решения многосеточных методов включают три основные операции [3]: smoothing (сглаживание), restriction (проекция), prolongation (интерполяция). Операция smoothing осуществляет уменьшение невязки численного решения. Во время операции restriction осуществляется проекция невязки на следующий уровень более грубой сетки. Во время операции prolongation осуществляется интерполяция коррекций для численного решения на более мелкую сетку, находящуюся уровнем выше.
При реализации многосеточного метода для задания структуры итерационного процесса построения численного решения будем использовать V-цикл [3], схема которого изображена на Рисунке 2.2. Белым кругом обозначена операция сглаживания, нисходящими стрелками обозначена операция restriction, восходящими стрелками обозначена операция prolongation. Процесс построения решения начинается с выполнения операции smoothing на сетке первого уровня с самым высоким разрешением. Затем выполняется операция restriction, переход на сетку следующего уровня и операция smoothing на данной сетке. Процесс продолжается до тех пор, пока не будет достигнута сетка последнего уровня. Затем выполняется операция prolongation, переход на сетку уровнем выше, коррекция значения для скорости и давления и операция smoothing на данной сетке. Последовательность переходов заканчивается на сетке первого уровня.
Будем использовать трехмерную смещенную сетку для давления, вязкости, компонент скорости и плотности жидкости. Ячейка сетки (Рисунок 2.3) представляет собой прямоугольный параллелепипед с ребрами Ar,A ?,Az вдоль координатных осей r, p,z. Для нумерации узлов сетки будем использовать индексы (j,j,l). Значения вязкости задаются для узлов, совпадающих с серединами ребер ячейки. Значения плотности задаются для узлов, находящихся в вершинах ячейки. Узел, для которого вычисляется значение давления, находится в центре ячейки. Узлы, для которых вычисляется значение скорости, совпадают с
Анализ решений и метод тестирования численных алгоритмов
Построим численное решение многосеточным методом в этой области. Значения компонент скорости на границе области приравниваются к значениям, определяемым аналитическим решением.
На Рисунке 2.17 представлен результат тестирования программы Stokes3D-spher, реализующей описанный многосеточный метод. В верхнем ряду расположены распределения для трех компонент скорости и давления в сечении = 0.5 для численного решения. В среднем ряду расположены аналогичные распределения для аналитического решения. В нижнем ряду представлены распределения разности между численным решением и аналитическим для компонент скорости и давления. На Рисунке 2.18 показана зависимость относительной ошибки многосеточного метода от величины шага сетки d в логарифмическом масштабе. Положительный тангенс угла наклона линий свидетельствует о сходимости численного решения к аналитическому при измельчении сетки.
Одним из самых эффективных подходов к численному моделированию геофизических процессов в мантии Земли является использование многосеточных методов [3]. В главе детально описан способ построения численных схем для многосеточного метода решения уравнений Стокса с переменной вязкостью в цилиндрической и сферической системах координат. Качество схем было проверено путем сравнения численного решения, полученного в результате работы программы, с аналитическим решением. Вычислительный процесс сходится при измельчении сетки и при малом, и при большом контрасте вязкости.
Тестирование алгоритма при различных вариантах поведения вязкости показало его хорошее качество даже для случая высокого контраста вязкости. Основные погрешности наблюдаются в окрестности изломов границы. Наибольшая погрешность наблюдается при расчете давления. Наибольшая относительная ошибка - Lx. Это связано с неравномерным распределением ошибки по области. Наибольшая ошибка имеет место в углах области, что и приводит к росту ошибки типа Lro. Разумеется, проведенное тестирование, как и любое тестирование вообще, не дает абсолютно точного ответа на вопрос о качестве алгоритма, а показывает его хорошую работу лишь в сходных с тестовыми ситуациях. Глава 3 Аналитические решения для плоских течений
В данной главе рассматриваются два вида плоских течений жидкости, наблюдаемых в процессах химического травления [10] и стеклообразования. Такие течения моделируются уравнениями Стокса с постоянной вязкостью. Ставится задача получить аналитические решения уравнений, записанных для функции тока. Первая задача состоит в нахождении в явном виде функции тока, описывающей течение над полостью специального вида. Второй задачей является получение функции тока, описывающей течение в области, ограниченной двумя половинками эллипса и двумя отрезками прямой. При решении краевой задачи для функции тока используется метод конформных отображений и метод разделения переменных.
Рассмотрим двумерное стоксово течение жидкости над полостью, представляющей собой выемку в форме половинки эллипса в безграничной полуплоскости. Область ограничена двумя лучами L 1,L2 и половинкой эллипса S (Рисунок 3.1). Стенки полости продолжены до фокусов эллипса. Жидкость находится в выделенной части плоскости. При рассмотрении течений будем считать, что на лучах L 1,L2 выполняются условия свободного проскальзывания [3]. На границе S выполнены условия прилипания жидкости. Граница S движется, мгновенная скорость ее точек направлена по касательной к границе. Таким образом, движение точек границы на Sзадается условием для тангенциальной компоненты скорости:
Дискретизация уравнений Стокса и неразрывности
Для определения решения остается найти коэффициенты ах,а2,а3,а4 для функции тока в представлении (3.13) и ЪХ,Ъ2,Ъ3,Ъ4 для функции тока в представлении (3.16) такие, чтобы {х,у) удовлетворяла граничным условиям и условиям склейки (3.18).
Условия (3.7), (3.10) выполняются при любых значениях коэффициентов. Подставив (3.13) в условие (3.8) и приравняв значения функций при членах ряда, содержащих sin(Jt&), к нулю, имеем п линейных уравнений для коэффициентов в представлении (3.15):
Значения коэффициентов с, могут быть найдены численно. Приравняв коэффициенты при функциях sin(Jt) в выражении (3.20), получаем п линейных уравнений для коэффициентов в представлении (3.15):
Далее, коэффициенты о1,а2,а3,«45 Ьг,Ь2,Ь3,Ь4 находятся путем решения системы линейных уравнений, составленной из наборов уравнений (3.19), (3.21), (3.22). Таким образом, для заданной функции v(2), определяющей скорость движения границы S, могут быть в явном виде получены выражения (3.13), (3.16), задающие функцию тока. 3.1.3 Картина линий тока для течений над полостью
В предыдущем разделе были получены выражения (3.13), (3.16), задающие функцию тока. Изолинии функции тока являются линиями тока жидкости. На Рисунке 3.4 представлена картина линий тока жидкости, для случая, когда скорость на границе S задается функцией от координаты 2 в виде F(2) = sin(2).
При вычислении функции тока в формулах (3.13), (3.16) выбирается параметр п = 5; а,Ъ - большая и малая полуоси эллипса.
Рассмотрим двумерное стоксово течение жидкости в области, ограниченной половинками S1 ,S2 двух софокусных эллипсов и двумя отрезками AB , CD, лежащими на прямой, проходящей через фокусы этих эллипсов (рис. 3.6). Введем декартовы (x, y) координаты с нулем в середине отрезка, соединяющего фокусы эллипсов (Рисунок 3.6). Остается найти коэффициенты а„а2,а3,а4 для функции тока в представлении (3.29) такие, чтобы W&y) удовлетворяла граничным условиям данной задачи. Условия (3.23),(3.27),(3.28) выполняются при любых значениях коэффициентов. Подставив (3.29) в условие (3.24) и приравняв значения функций при членах ряда, содержащих sin(fc), к нулю, имеем 2п линейных уравнений для коэффициентов в представлении (3.29):
Коэффициенты а1,а2,а3,а4, могут быть найдены путем решения линейной системы уравнений. В эту систему входят уравнения (3.30), (3.31), обеспечивающие выполнение условия (3.24), и часть уравнений из множеств (3.33), (3.34), задающие значение тангенциальной скорости на подвижной части границы. В систему можно включить, например, все уравнения из множества (3.33) и часть уравнений из (3.34) или все уравнения из множества (3.34) и часть уравнений из (3.33). Таким образом, при выводе решений (3.29) в данной задаче можно задать произвольную скорость только на одной из границ S1,S2, а на противоположной границе удается поставить ограничения на часть коэффициентов в выражениях (3.15).
В предыдущем разделе описан вывод выражений (3.29), задающих функцию тока. На рисунках ниже представлена картина линий тока жидкости, для различных видов зависимости скорости на участках границы S1, S2 от координаты 2 на промежутке [я-, 2л-]. Вычисления производились при следующих значениях параметров:а = 0.4;Ъ = 0.12 для S1. а = 0.86; 6 = 0.76 для S2, где а,Ъ - большая и малая полуоси эллипса.
Картина линии тока для течений над полостью
В предыдущем разделе описан вывод выражений (3.29), задающих функцию тока. На рисунках ниже представлена картина линий тока жидкости, для различных видов зависимости скорости на участках границы S1, S2 от координаты на промежутке [я-, 2л-]. Вычисления производились при следующих значениях параметров:а = 0.4;Ъ = 0.12 для S1. а = 0.86; 6 = 0.76 для S2, где а,Ъ - большая и малая полуоси эллипса.
Приведем пример тестирования конечно-разностной схемы на полученных аналитических решениях (3.13),(3.16) и (3.29). Построим численное решение уравнения (3.1) методом конечных разностей для течения над полостью и для течения в области, имеющей форму эллиптической подковы. Затем выполним сравнение численного решения с аналитическим.
При построении численного решения будем использовать прямоугольную сетку с квадратными ячейками dxd. Для записи уравнения (3.1) в конечно-разностной форме воспользуемся шаблоном, включающим 13 узлов сетки (Рисунок 3.18).
При построении численного решения задачи о течении над полостью ограничим область в верхней полуплоскости некоторым прямоугольником R, заданным неравенствами:-/! х 1{, 0 у 12. На верхней и на боковых частях границы R зададим значения для численного решения, равные значениям аналитического решения (3.16).
На Рисунке 3.19 приведен график зависимости относительной ошибки численного решения є, вычисленной с нормой LY, от величины ребра ячейки сетки d для течения над полостью. Вычисления производились для конфигурации системы, приведенной к Рисунку 3.4. (V() = sin(), a = 0.9,
На Рисунке 3.20 приведен график зависимости относительной ошибки с нормой от величины ребра ячейки сетки d для течения в области, имеющей L1 форму эллиптической подковы. Вычисления производились для конфигурации системы, описанной в примере №2 предыдущего раздела. Рисунок 3.20. Зависимость относительной ошибки от шага сетки для течения, описаного в примере №2. Положительный тангенс угла наклона линий графиков на Рисунках 3.19, 3.20 свидетельствует о сходимости численного решения к аналитическому при измельчении сетки. 3.4 Выводы
В данной главе строятся эталонные решения для моделей процессов химического травления и стеклообразования. А именно, описан способ построения аналитических решений уравнений Стокса для двумерного течения жидкости в областях двух типов. В первом случае область представляет собой безграничную полуплоскость с выемкой в форме половины эллипса. Во втором случае область замкнута и имеет форму эллиптической подковы. Аналитические решения могут быть построены для различных видов функций, задающих скорость подвижной части границы. В главе представлены картины линий тока для разных конфигураций системы.
Полученные решения когут бать использованы в качестве эталона при тестировании численных алгоритмов решения уравнений Стокса. Подобное тестирование, основанное на сравнении с аналитическими решениями, применялось в работах [10,40]. В главе приведен пример тестирования численного алгоритма решения бигармонического уравнения методом конечних разностей. Построенные решения также позволяют моделировать вихревую структуру течения и выбирать параметры для обеспечения требуемого характера течений. Заключение
Основные результаты работы заключаются в следующем: 1. Для моделей, описывающих геофизические течения в мантии Земли и нанотечения в тубулярных наноструктурах, построены частные аналитические решения уравнений Стокса и неразрывности с переменной вязкостью и плотностью в цилиндрической и сферической системах координат в случае, когда вязкость и плотность являются функциями радиуса. В сферической системе координат решения найдены для вязкости, являющейся степенной функцией радиуса. Предложен метод, позволяющий на основе аналитических решений составлять тестовые краевые задачи для оценки качества работы численных алгоритмов решения уравнений Стокса с переменной вязкостью и плотностью. 2. Предложены и протестированы конечно-разностные схемы решения уравнений Стокса и неразрывности с переменной вязкостью в цилиндрической и сферической системах координат для моделей геофизических течений в мантии Земли и нанотечений в тубулярных наноструктурах. Решения уравнений строятся многосеточным методом на последовательности ортогональных смещенных сеток. Качество работы схем проверено путем сравнения численного решения, полученного в результате работы программы, с эталонными аналитическими решениями. 3. Для моделей, описывающих процесс химического травления и стеклообразования, предложен способ построения аналитических решений уравнений Стокса, задающих двумерное течение над полостью и в области, имеющей форму эллиптической подковы. Представлена картина линий тока для различных конфигураций системы. Продемонстрирована возможность использования полученных решений при тестировании численных алгоритмов. Предложенный подход позволяет моделировать вихревую структуру течений при различных параметрах системы.