Решение многих современных прикладных задач приводит к системам линейных алгебраических уравнений (СЛАУ) с седловой точкой. Характерной особенностью таких систем является знаконеопределенность. В симметричном случае имеются как положительные, так и отрицательные собственные значения. Подобные системы возникают, например, в теории динамических систем [37], в теории упругости [22, 28], в экономике [16, 32, 42], в теории оптимального управления [17, 18], при моделировании электромагнитных явлений [21, 54] и во многих других областях. Важной областью, требующей решения задач с седловой точкой, является численное решение линеаризованных уравнений Навье-Стокса, описывающих течение несжимаемой вязкой жидкости [38, 39, 55, 62, 63, 67]. Уравнения Навье-Стокса являются основными уравнениями гидродинамики и, соответственно, играют важную роль в современной науке. Отметим, что в ходе их решения, как правило, возникает необходимость проводить вычисления на мелких сетках, а значит решаемые системы имеют большую размерность. Знаконеопределенность делает процедуру выбора метода решения таких систем нетривиальной, так как многие эффективные методы решения СЛАУ, см., например, [58], пригодны только для систем с положительно определенной матрицей.
Уравнения Навье-Стокса во многих случаях хорошо описывают поведение жидкостей и газов. Однако, многие вещества в природе описываются моделями с переменным коэффициентом вязкости, зависящим от внешних факторов. Примером могут служить биологические жидкости (например, кровь), нефть, зубная паста, кетчуп, крахмал, разведенный в воде и многие другие вещества. В некоторых случаях, например, кетчуп, наблюдается уменьшение вязкости при возрастании внутренних напряжений. В других, таких как крахмал, поведение обратное — при возрастании напряжений вязкость резко увеличивается. Для моделирования подобных подобных веществ можно рассматривать модифицированные уравнения Навье-Стокса, при этом вязкость является не постоянным параметром среды, а функцией, зависящей от динамических, кинематических или других характеристик среды в данной точке пространства, например, тензора скоростей деформации, давления, температуры, и т. д. Возникающие при дискретизации линейные системы характеризуются тем, что их число обусловленности зависит от отношения максимальной вязкости к минимальной. Данное обстоятельство предъявляет дополнительное требование к методам решения СЛАУ, а именно — независимость числа итераций от отношения максимального значения коэффициента вязкости к минимальному и от градиента коэффициента вязкости как функции пространственной переменной. Эффективное решение уравнений Навье-Стокса с переменной вязкостью является темой ряда современных исследований. Упомянем некоторые из них. В статье [70] рассматривается моделирование мантийной конвекции в планетарном масштабе. Рассматривается влияние разных моделей вязкости, в том числе — распределенной по слоям, а также вязкости, зависящей от температуры. В [49] получены оценки сходимости конечноэлементного решения задачи Стокса с интерфейсом к решению непрерывной задачи. Вязкость при этом кусочно-постоянна и имеет скачок. Работа [56] посвящена построению эффективных итерационных методов для несжимаемой задачи Стокса с переменной вязкостью, в ней рассматриваются задачи, связанные с моделированием процесса экструзии и с геологией. В статье [68] представлено построение равномерно устойчивых по скачку коэффициента вязкости конечноэлементных аппроксимаций уравнений Дарси-Стокса-Бринкмана. В работе [52] получено тЛ^эир неравенство для задачи Стокса с интерфейсом, которое равномерно по скачку вязкости между подобластями. В [51] авторами предложен метод анализа переобуславливателей для дополнения по Шуру для абстрактных задач с седловой точкой. Применение этого метода к задаче Стокса с интерфейсом приводит к построению переобуславливателя, являющегося устойчивым по отношению к ряду параметров.
В численном анализе методов решения задачи Стокса важную роль играет условие ЬВВ (Ладыженской-Бабушки-Брецци) [23] и его непрерывный аналог, неравенство Нечаса [46]. В диссертации неравенство Нечаса обобщается на случай переменной вязкости и на основе этого обобщения получены оценки эффективности предлагаемого итерационного метода. Некоторые другие обобщения неравенства Нечаса получены в работах [1] и [6].
Одной из задач, где возникают уравнения Навье-Стокса с переменной вязкостью, является моделирование среды Бингама. Среда Бинга-ма является вязко-пластичной средой, которая при напряжениях ниже порогового значения ведет себя как твердое тело, а при превышающих пороговое значение — как вязкая жидкость. Примерами веществ, описываемых моделью Бингама являются зубная паста, различные суспензии, например грязь или незастывший бетон, геоматериалы. Впервые модель независимо была предложена Шведовым [60] и Бингамом [20] для описания движения суспензий в условиях чистого сдвига (одномерная модель). Позднее Генки [2] и Ильюшин [4] предложили пространственное обобщение уравнения состояния Шведова-Бингама и решили ряд задач для случая плоских течений вязко-пластичной среды. Течению среды.
Бингама посвящено большое количество литературы, среди отечественной можно отметить труды П. М. Огибалова и А. Х. Мирзаджанадзе [10], П. П. Мосолова и В. П. Мясникова [9], а также Д. М. Климова, А. Г. Петрова, Д. В. Георгиевского. [5]. Задача Бингама является более сложной, чем задача Навье-Стокса. Одним из подходов к ее численному решению является регуляризация — модель, когда среда Бингама рассматривается как жидкость с переменной вязкостью.
Задача Бингама трудно поддается математическому анализу и ее точные решения найдены только для узкого круга модельных задач, например, для течения среды между двумя параллельными пластинами [5]. По этой причине численные методы, зачастую, являются единственным способом анализа многих процессов. Численному решению задачи Бингама посвящены, например, следующие публикации: [19, 29, 30, 33, 44, 48, 53, 57, 65].
Задачи с переменной вязкостью появляются и во многих других научных областях, например, в геологии. В мантии Земли температура неоднородна, а вязкость магмы напрямую зависит от температуры. Численные методы позволяют проводить моделирование внутреннего строения Земли и эта тема рассматривается сейчас многими учеными. Например, в [45, 69] обсуждаются вопросы численного моделирования мантийной конвекции, в [43] — вопросы численного моделирования миграции магмы.
Цель работы.
Работа преследует следующие цели.
1. Построение эффективного метода решения систем линейных алгебраических уравнений, возникающих при дискретизации задачи Стокса с переменной вязкостью. Поскольку в реальных ириложениях отношение максимальной вязкости к минимальной может быть очень большим, к такому методу решения СЛАУ в настоящей работе предъявляется требование независимости (или слабой зависимости) количества итераций от этого отношения, а также от шага сетки.
2. Теоретический анализ эффективности предлагаемого метода решения систем линейных алгебраических уравнений. В частности, получение оценок скорости сходимости в терминах экстремальных значений коэффициента вязкости и других параметров систем уравнений.
3. Проверка эффективности предлагаемого метода на модельных задачах: регуляризованной задаче моделирования течения среды Бингама в канале и каверне, а также на линейной задаче, возникающей при моделировании всплытия раскаленного пузыря в мантии Земли.
Используемые обозначения.
• Обычным шрифтом (например, р) обозначаются скалярные функции, коэффициенты разложения по конечноэлементным базисам и вектора в матричной записи линейных систем.
• Полужирным шрифтом (например, V) обозначаются векторные и тензорные функции и операторы. Используются обозначения.
Уи = дих дих дх ду диу диу дх ду дих диг дъ диг дх ду сИУБ = а^Б1 ^ сИУ Г>2 у } где Ю — тензорная функция, действующая из К3 в М3×3, и Т>к =.
• Индекс Н у функции или оператора указывает на дискретную (ко-нечноразностную или конечноэлементную) аппроксимацию данной функции или данного оператора.
• Через 1Р обозначается пространство Лебега с индексом суммирования р, через || • \ьр — норма в I/. Запись || • || означает норму в Ь2. Скалярное произведение в Ь2 и евклидово скалярное произведение обозначаются как (•,•) и (•,¦), соответственно.
Используется — пространство всех функций / из таких, что /сЮ, = 0 и — пространство функций / из Ь2{0), для которых выполняется условие /ь>~1сЮ, = 0 с некоторой функцией и > 0 такой, что V" 1 € Ь°°.
• Символом Нц обозначено пространство интегрируемых вектор-функций, градиент которых, понимаемый в обобщенном смысле, принадлежит Ь2(Жах<�с/ = 2,3 и имеющих нулевой след на границе расчетной области.
Краткое содержание работы.
В первой главе приводится формулировка задачи Стокса с переменной вязкостью, доказывается обобщение неравенства Нечаса на случай весовой нормы. Приводится матричная постановка дискретной задачи, предлагается переобуславливатель для дополнения по Шуру и и при помощи обобщенного неравенства Нечаса получаются оценки его эффективности.
В первом разделе дается формулировка задачи Стокса с переменной вязкостью. Основные уравнения, рассматриваемые в первой главе, имеют вид.
— divi/(x)Du + Vp = f в ft —div u = 0 в Q ii = U () на dQ,.
Через Du обозначен тензор скоростей деформации |(Vu + VTu). Вязкость v в первой главе считается зависящей только от пространственных координат х.
Во втором разделе получен основной теоретический результат диссертации — неравенство Нечаса в весовой норме. Оно является обобщением хорошо известного неравенства Нечаса (g, divv) О < с0 < mi sup и «» и,veHjIMI IIVvH' которое, в свою очередь, является непрерывным аналогом условия Ладыженской — Бабушки — Брецци (LBB): если константа qh, divv/j) со, ft = mf sup т—?г=—?7 чьеQhVheYh Ш llVvftll равномерно по шагу сетки отделена от нуля, то пара конечноэлементных пространств Уд и Q/j (для скорости и давления, соответственно) приводит к устойчивой дискретной задаче. Под равномерным отделением от нуля понимается существование такой константы с* > 0, что при всех h выполняется co, ft > с*. Результат сформулирован в виде Теоремы 1, суть утверждения которой заключается в выполнении неравенства для функций д? Ь2(С1), на которые наложены условия ортогональности.
Приводится важное следствие Теоремы 1 — обобщение на случай, когда расчетная область Г2 представлена в виде объединения конечного числа подобластей Г^. В этом случае константа си для всей области О, может быть оценена снизу минимумом из всех констант сДГ^), вычисленных на подобластях.
Иногда при помощи удачного разбиения удается добиться того, что пищ с"(П). Ценой этого улучшения являются два условия ортогональности на каждой подобласти: (д, — (<7> и -½)щпг) = о.
Разделы 3 и 4 посвящены дискретизации задачи и построению СЛАУ. Дискретизация проводилась двумя методами, а именно — методом конечных элементов (?боР2-Р1) и методом конечных разностей (МАС-схема). Дается описание построения дискретной задачи при помощи обоих методов. Излагается вариационная постановка, дается определение конечно-элементных пространств для метода конечных элементов, вводятся сетки и задаются сеточные аналоги для непрерывных операторов при использовании метода конечных разностей. Оба метода приводят к решению системы линейных алгебраических уравнений с седловой точкой с матрицей вида где, А = Ат > 0.
В пятом разделе описывается метод решения полученной дискретной 0. Получены оценки на константу с1 V задачи. Основная идея заключается в использовании метода на подпространствах Крылова (MINR.ES или СМНЕЯ) со специальным переобу-славливателем вида.
Р=1 Л ° I (3).
О 5 или.
7 А О.
4).
Здесь, А — переобуславливатель для блока А, 5' — переобуславлива-тель для дополнения по Шуру 5 = ВА~гВт. Для дополнения по Шуру предлагается применять переобуславливатель.
М"}ц := (5) при использовании метода конечных элементов (через ф{ обозначены базисные функции для давления) и (6) при использовании метода конечных разностей. Переобуславливатель Ми при построени учитывает функцию что важно при большом отношении максимального значения вязкости к минимальному. Сформулированы дискретные аналоги Теорем 1 и 2, обозначенные в тексте как Предположение 1 и Теорема 4, соответственно. Предположение 1 сводится к предположению справедливости дискретного аналога обобщенного неравенства Нечаса для пары конечноэлементных пространств, удовлетворяющих ЬВВ условию устойчивости. Далее доказывается, что.
Зд, д) <(1{Мид, д) для всех qh 6 Qh и, если Предположение 1 справедливо, то выполняется соотношение cl, h (Mvq, q) < (Sq, q) для всех qh? Qh таких, что (qh, v~l)щ^) = (qh, = 0. Здесь d — размерность пространства, т. е. d = 2 или d = 3, через q обозначен вектор коэффициентов разложения по базису функции qh? Q/г, а с^л — константа из Теоремы 4. Таким образом устанавливается связь между неравенством Нечаса в весовых нормах и эффективностью переобуслав-ливателя Ми для дополнения по Шуру.
В шестом разделе первой главы дается оценка собственных значений для переобусловленной матрицы дополнения по Шуру. Она выводится из представления Куранта-Фишера (Sg, д) ли = max min ——г,.
CeVfc 1 qeic± (Mvq, q) где символом Vfci обозначены все (к — 1)-мерные подпространства Rm. Для (M~lS) в диссертации получены оценки.
0 = Ai< А2лг+ь (V где d — размерность пространства, N — количество подобластей, на которые разбита область, а с^д — константа из Теоремы 4.
В этом разделе также приводится оценка скорости сходимости метода Узавы-сопряженных градиентов [34] для систем линейных алгебраических уравнений с матрицей (2) и переобуславливателем для дополнения по Шуру Ми и оценка количества дополнительных итераций, вызванных возможным присутствием 2N — 1 малого собственного значения.
Во второй главе описывается применение теоретических результатов первой главы к решению двух модельных задач: задачи Бингама о течении вязкопластичной жидкости в канале и задаче моделирования мантийной конвекции.
В первом разделе даются определяющие соотношения модели Бинга-ма, вводятся понятия жестких зон и эффективной вязкости в жидких зонах. Задача Бингама нелинейна, и ее эффективная вязкость зависит от тензора скоростей деформации, а значит, и от поля скоростей и.
Соотношения, связывающие тензор скоростей деформации Du и де-виатор тензора напряжений т, имеют вид.
Dijи rlJ = 2(iDlJu + если |Du| > 0;
M < т8, если |Du| = 0, где т3 — предел пластичности, а ?1 — коэффициент вязкости. В жидкой зоне коэффициент пропорциональности между девиатором тензора напряжений и тензором скоростей деформаций часто называют эффективной вязкостью. При этом в жестких зонах эффективную вязкость формально можно считать бесконечной.
В диссертации рассматривается подход к численному решению задачи Бингама, основанный на регуляризации, последующей дискретизации на каждом шаге нелинейного процесса и применении нового переобусловленного итерационного метода для решения линейных вспомогательных задач. Важным свойством предлагаемого в настоящей диссертации метода является отсутствие параметров, требующих тонкой настройки.
Во втором разделе дается описание регуляризованной задачи. Регуляризация — подход к решению задачи Бингама, при котором среда рассматривается как жидкая во всей расчетной области, а в жестких зонах эффективная вязкость считается конечной. Величина вязкости в жестких зонах существенным образом определяется параметром регуляризации е. Рассматриваются два варианта регуляризации: г, г/е (|Ви|) = 2[1 Н— я = — Берковьера-Энгельмана;
Т>и2 + е2 |.
1 ехр (Шн1) г/е (|Ви|) = 2/л + т3-|—-— — Папанастасио.
В третьем разделе определяется итерационный метод, используемый для решения регуляризованной задачи Бингама. Нелинейные уравнения решаются методом Пикара:
7 и" -1 I { и71 (и" -1, Р'" / рп~г и" -1) О рп~1 / о где через Р обозначен оператор, линейный при фиксированной функции, а ч (-а^(|Ба|)В V 24 а) =.
УсНУ О, а через ^(ип-1)-1 — приближенное решение линейной задачи с матрицей Л (2), возникающей при дискретизации оператора Р.
В четвертом разделе выводится оценка собственных значений матрицы М~1Б для задачи о течении среды Бингама в канале. Это одна из немногих задач, где известно аналитическое решение и = (и, у):
1 — 2т3)2,. если — т3 < у < + та, и = <1 [(1 — 2т3)2 — (1 — 2т3 — 2у)2], если 0 < у < - т&bdquo- [(1 — 2т3)2 — (2у — 2та — I)2], если 1 > у > § + г",.
V = 0, р = —х.
С помощью разбиения области на три подобласти, соответствующие ядру течения в центре и жидким зонам по краям, непосредственно оценена константа в обобщенном неравенстве Нечаса (1) и, используя соотношения (7), получена оценка на собственные значения переобусловленного матрицей Mv дополнения по Шуру: для любого s G (0,1] существует не зависящая от? константа С2, что.
С£ < C2? s < Л7 < • • • < Am < d = 2, где константа с не зависит от е.
Последний раздел второй главы посвящен применению результатов первой главы к задаче моделирования мантийной конвекции. Рассматривается линейная задача, возникающая при моделировании всплытия раскаленного пузыря в магме. В диссертации не ставилась цель решать полную физическую задачу, для простоты считалось, что вязкость является функцией только от пространственных координат. Для данной задачи приводится алгоритм разбиения области на подобласти и вычисленные значения констант в обобщенном неравенстве Нечаса (1).
Таким образом, во второй главе описано применение метода, предлагаемого в первой главе, к модельным задачам и получены оценки его эффективности.
В третьей главе представлены результаты численного решения модельных задач, расмотренных во второй главе и задачи течения среды Бингама в каверне.
В первых двух разделах представлены и проанализированы результаты решения задач течения среды Бингама в канале и в каверне, соответственно. Приводятся графики, на которых изображены все собственные значения переобусловленного дополнения по Шуру с использованием как матрицы масс М, так и предлагаемого в работе нереобуславливателя Mv. Приводится количество нелинейных и линейных итераций для разных конфигураций метода решения СЛАУ. Расчеты проводились с обеими ре-гуляризациями с разными значениями регуляризационного параметра е. Дано также сравнение блочно — диагонального (3) и блочно — треугольного (4) переобуславливателей для матрицы Л. Вычисления показали, что для данной модельной задачи численное решение хорошо приближает аналитическое и выбор регуляризации слабо влияет на эффективность метода. При этом, использование блочно — треугольного переобуславли-вателя (4) приводит к заметно более высокой скорости сходимости, чем использование блочно — диагонального (3).
Для задачи о каверне представлены графики с приближениями жестких зон для различных значений предела пластичности т3 и регуляризационного параметра е. Эти результаты (вид и размеры жестких зон) хорошо согласуются с результатами расчетов в других работах. В [44], например, рассматривается регуляризованная задача, а в работах [30] и [50] — нерегуляризованная. Это свидетельствует о том, что выбранные численные параметры позволяют воспроизводить важные физические эффекты.
Из численных экспериментов по решению регуляризованной задачи Бингама в канале и в каверне можно сделать следующий вывод: предлагаемый в диссертации итерационный метод решения систем линейных алгебраических уравнений показывает практически полное отсутствие зависимости числа итераций как от шага сетки, так и от параметра регуляризации (соответственно, от отношения максимального значения вязкости к минимальному), что хорошо согласуется с теоретическими оценками (в задаче о течении в канале).
В третьем разделе третьей главы обсуждаются результаты численного решения линейной задачи, возникающей при моделировании всплытия раскаленного пузыря в магме. Сравнивается эффективность пере-обуславливателей М и Mv. Численные эксперименты показали, что число линейных итераций в задаче моделирования мантийной конвекции с увеличением отношения максимальной вязкости к минимальной растет незначительно. Более того, предлагаемый итерационный метод оказывается эффективнее, чем предсказывают теоретические оценки.
Апробация работы.
1. Международная конференция «Х-я Белорусская математическая конференция». Белоруссия, Минск, 2008.
2. Доклад на семинаре Кафедры вычислительной математики Математического факультета Технического университета Дортмунда под руководством Ш. Турека. Германия, Дортмунд, 2009.
3. XVI Международная конференция студентов, аспирантов и молодых ученых «Ломоносов». Москва, 2009.
4. International Workshop on «Computational Mathematics and Applications». Финляндия, Тампере, 2009.
5. XVII Международная конференция студентов, аспирантов и молодых ученых «Ломоносов». Москва, 2010.
6. Доклад на семинаре Кафедры вычислительной математики, Механико-математический факультет МГУ под руководством Г. М. Ко-белькова. Москва, 2010.
7. Доклад на семинаре Кафедры механики композитов под руководством В. И. Горбачева, Механико-математический факультет МГУ. Москва, 2010.
8. Доклад на семинаре «Технологии математического моделирования течений со свободной границей» под руководством Ю. В. Василевского и М. А. Ольшанского, ИВМ РАН. Москва, 2010.
9. Доклад на семинаре Кафедры математического моделирования МЭИ (ТУ) под руководством А. А. Амосова и Ю. А. Дубинского. Москва, 2010.
На защиту выносятся следующие основные результаты и положения:
1. Доказано обобщенное неравенство Нечаса для случая переменной вязкости. Доказано обобщение неравенства для случая, когда область представлена в виде объединения непересекающихся подобластей.
2. Предложен иереобуславливатель для дополнения по Шуру для дискретной задачи Стокса, учитывающий переменную вязкости. Получена оценка на собственные значения переобусловленного дополнения по Шуру. Получена оценка скорости сходимости метода Узавы-сопряженных градиентов.
3. Предлагаемый переобуславливатель применен для численного решения регуляризованной задачи Бингама и линейной задачи, возникающей при моделировании всплытия раскаленного пузыря в магме. Для задачи о течении среды Бингама в канале получены оценки собственных значений переобусловленного дополнения по Шуру и теоретически доказано, что зависимость от параметра регуляризации оценок собственных значений для предлагаемого метода меньше, чем для обычно используемого переобуславливателя на основе матрицы масс.
4. Проведены численные эксперименты, хорошо согласующиеся с теоретическими результатами задачи течения среды Бингама в канале и показавшие практическую независимость числа линейных итераций как от шага сетки, так и от отношения максимального значения вязкости к минимальному для задачи течения среды Бингама в каверне и линейной задачи, возникающей при моделировании всплытия раскаленного пузыря в магме при использовании предлагаемого переобуславливателя.
Публикации.
1. P.P. Grinevich, М.А. Olshanskii. An iterative method for the Stokes type problem with variable viscosity // SIAM Journal on Scientific Computing, Vol.31, Iss.5, 2009, pp.3959−3978.
2. U.U. Гриневич, М. А. Ольшанский. Итерационный метод решения регуляризованной задачи Бингама // Вычислительные методы и программирование, Том 10, 2010, стр. 78−87.
3. П. П. Гриневич. Об итерационном методе решения задачи Стокса с переменной вязкостью // Вестник МГУ, № 3, 2010, стр. 38−41.
Структура и объем диссертации
.
Диссертация состоит из введения, трех глав, заключения и списка литературы, содержащего 70 наименований. Общий объем работы — 105 страниц, работа включает 22 иллюстрации и 18 таблиц.
Основные результаты, полученные в диссертации, могут быть сформулированы следующим образом:
1. Доказано обобщенное неравенство Нечаса для случая переменной вязкости. Доказано обобщение неравенства для случая, когда область представлена в виде объединения непересекающихся подобластей.
2. Предложен переобуславливатель для дополнения по Шуру для дискретной задачи Стокса, учитывающий переменную вязкость. Получена оценка на собственные значения переобусловленного дополнения по Шуру. Получена оценка скорости сходимости метода Узавы-сопряженных градиентов.
3. Предлагаемый переобуславливатель применен для численного решения регуляризованной задачи Бингама и линейной задачи, возникающей при моделировании всплытия раскаленного пузыря в магме. Для задачи о течении среды Бингама в канале получены оценки собственных значений нереобусловленного дополнения по Шуру и теоретически доказано, что зависимость от параметра регуляризации оценок собственных значений для предлагаемого метода меньше, чем для обычно используемого переобуславливателя на основе матрицы масс.
4. Проведены численные эксперименты, хорошо согласующиеся с теоретическими результатами задачи течения среды Бингама в канале и показавшие практическую независимость числа линейных итераций как от шага сетки, так и от отношения максимального значения вязкости к минимальному для задачи течения среды Бинга-ма в каверне и линейной задачи, возникающей при моделировании всплытия раскаленного пузыря в магме при использовании предлагаемого переобуславливателя.
Заключение
.
