Содержание к диссертации
Введение
2 Квазифотоны для волн на поверхности тяжелой жидкости 17
2.1 Энергетическое соотношение 17
2.1.1 Постановка задачи 17
2.1.2 Энергетическое соотношение 19
2.1.3 Вектор Умова S как “вектор работы”
2.2 Пространственно-временной лучевой метод для волн на поверхности тяжелой
2.3 Уравнения переноса
2.3.1 О решении уравнений переноса 27
2.3.2 Расчет высших приближений 30
2.4 Комплексный вариант лучевого метода и квазифотоны 30
2.4.1 Общая схема построения 30
2.4.2 Нулевой и первый порядки 32
2.4.3 Второй порядок 33
2.4.4 Расчет высших приближений 35
2.4.5 Заключение 35
3 Волновые валы на поверхности тяжелой жидкости 37
3.2 Исходные уравнения 38
3.3 Пространственно-временной лучевой метод 38
3.4 Построение в 40
3.5 Построение Ф.
3.5.1 Построение Фо 45
3.5.2 Построение Фі 51
3.5.3 Построение Ф.,- 52
4 Метод эталонных решений для исследования магматических течений в мантии Земли 54
4.1 Построение решений в двумерном случае 54
4.1.1 Формулировка уравнений с переменной вязкостью в двумерном случае 54
4.1.2 Линейно меняющаяся вязкость 56
4.1.3 Экспоненциально меняющаяся вязкость 58
4.2 Построение решений в трехмерном случае 59
4.2.1 Формулировка уравнений с переменной вязкостью в трехмерном случае 59
4.2.2 Линейно меняющаяся вязкость. Трехмерный случай 61
4.2.3 Экспоненциально меняющаяся вязкость. Трехмерный случай 63
4.3 Схема метода и проверка корректности работы конкретного алгоритма 65
4.3.1 Двумерный случай. Пример для линейно меняющейся вязкости 66
4.3.2 Двумерный случай. Пример для экспоненциально меняющейся вязкости 69
4.3.3 Трехмерный случай. Пример для линейно меняющейся вязкости 71
4.3.4 Трехмерный случай. Пример для экспоненциально меняющейся вязкости 5 Заключение 79 Список литературы 81
6 Приложения
- Пространственно-временной лучевой метод для волн на поверхности тяжелой
- Комплексный вариант лучевого метода и квазифотоны
- Построение Фо
- Линейно меняющаяся вязкость. Трехмерный случай
Введение к работе
Актуальность. Проблема описания движений жидкости возникает во многих физических задачах. Большой интерес представляют геофизические задачи, актуальность которых связана с необходимостью прогнозирования природных явлений и обеспечения функционирования различных технических систем. Данные задачи чрезвычайно сложны и обычно допускают только численное решение. В то же время, прямые расчеты без предварительного математического анализа часто оказываются недостаточно эффективными. Кроме того, возникают вопросы верификации численных результатов, проверки корректности методов, качества работы численных алгоритмов. Для этой цели могут использоваться аналитические результаты, полученные строгими математическими методами. Разработке именно таких подходов и посвящена данная диссертационная работа. Рассматриваются две задачи о движениях жидкости, порожденные геофизическими явлениями: распространение волн на поверхности океана и течения в мантии Земли.
В задаче о поверхности океана показано существование волн, локализованных в окрестности движущейся по поверхности точки или линии. Использован асимптотический подход. В основе анализа – пространственно-временной лучевой метод, предложенный в теории дифракции В.М.Бабичем, В.С.Булдыревым, И.А.Молотковым. Рассмотрена задача об океане с медленно меняющейся глубиной. Скорость изменения глубины взята в качестве малого параметра задачи, по которому строится асимптотическое разложение. Анализируются решения двух типов: квазифотон (волна, сосредоточенная в окрестности движущейся точки) и волновой вал (волна, сосредоточенная в окрестности движущейся линии). В обоих случаях построено полное асимптотическое разложение соответствующего решения. Проведено доказательство корректности процедуры на всех ее этапах. Данное решение, полученное аналитически в рамках строгого математического подхода, можно использовать в качестве эталонного в данной задаче.
В задаче о течениях в мантии Земли предложен другой подход к построению эталонных решений. Данные течения описываются уравнениями Стокса и неразрывности с переменными вязкостью и плотностью. Подобного типа уравнения возникают и в других областях физики, в частности, в задачах о течениях в наноструктурах. Для выбранных конкретных типов зависимостей (линейной и экспоненциальной) вязкости и плотности от координат найдены аналитически явные решения уравнений в двумерном и трехмерном случаях. Данные эталонные решения очень удобно использовать для тестирования вычислительных алгоритмов. При этом имеется возможность достаточно просто исследовать характер сходимости, скорость работы алгоритма и т.п. при различных параметрах системы или вариациях вычислительной схемы, выбирая наилучший вариант для данной задачи.
Проверка точности работы численных алгоритмов с помощью эталонных решений широко используется в геофизике. В основном, они получены в двумерном случае для ряда задач. В частности, в [B. Blankenbach et al. - Vol. 98. Geophys. J. Int., 1989, P. 23–38] проведено сравнение различных эталонных решений для устойчивой конвекции с постоянной вязкостью, конвекции с переменной вязкостью и зависящей от времени конвекции с внутренним
разогревом. В [P.E. van Keken et al. - Vol. 171 (1-4). Phys. Earth Planet. Inter., 2008, P. 187– 197] разработаны эталонные решения для проверки численных алгоритмов, описывающих течения и распределение тепла в зонах субдукции. В [H. Schmeling et al. - Vol. 171 (1-4). Phys. Earth Planet. Inter., 2008, P. 198–223] предложены эталонные решения для проверки численных алгоритмов расчета динамики литосферы с субдукцией для различных вычислительных кодов (эйлеровых и лагранжевых, конечно-элементных и конечно-разностных, с маркерами и без маркеров). Для трехмерного случая следует отметить работы [F.H. Busse et al. - Vol. 75. Geophys. Astrophys. Fluid Dyn., 1993, P. 39–59], где предложены эталонные решения для сравнения с численными методами исследования конвекции с высоким числом Прандтля в трехмерной декартовой геометрии, и [S. Zhong et al. - Vol. 9. Geochem. Geophys. Geosyst., 2008], где приведены эталонные решения для проверки численных методов расчета конвекции мантии в сферической геометрии.
Объектом исследования являются модельные дифференциальные уравнения, возникающие в геофизических задачах, связанных с динамикой жидкости. Предмет исследования - математические модели, описывающие волны на поверхности океана переменной глубины и течения в мантии Земли.
Цель диссертационного исследования – исследование волн на поверхности океана переменной глубины, сосредоточенных в окрестности движущихся точек и линий, асимптотическими методами и построение необходимых для тестирования вычислительных схем и алгоритмов эталонных решений уравнений Стокса и неразрывности с переменными вязкостью и плотностью.
Для достижения этой цели в диссертации:
-
Разработан алгоритм построения полного асимптотического разложения решения типа квазифотона, сосредоточенного в окрестности пространственно-временного луча для волн на поверхности океана переменной медленно меняющейся глубины, проведено обоснование алгоритма с доказательством реализуемости всех его шагов;
-
Разработан алгоритм построения полного асимптотического разложения решения типа волнового вала, сосредоточенного в окрестности поверхности, сотканной из пространственно-временных лучей, для волн на поверхности океана переменной медленно меняющейся глубины, доказана разрешимость задач, возникающих на всех этапах схемы;
-
Построены в аналитическом виде с использованием метода характеристик эталонные решения уравнений Стокса и неразрывности с линейно и экспоненциально зависящими от координат вязкостью и плотностью в двумерном и трехмерном случаях,
-
Реализована методика тестирования вычислительного алгоритма решения краевой задачи для уравнений Стокса и неразрывности с помощью сравнения результатов расчета численными методами с построенными явно эталонными решениями.
Теоретическая и практическая значимость работы.
В работе получены новые теоретические результаты, важные для построения математических моделей и теоретического описания геофизических явлений. А именно, показано существование нового типа волн на поверхности океана – высокочастотных волновых паке-
тов, сосредоточенных в окрестности движущейся по этой поверхности точки (квазифотон) или линии (волновой вал) и экспоненциально убывающих при удалении от них. В отличие от солитонов, это линейные волны. Для данных решений построены полные асимптотические разложения по малому параметру - скорости изменения глубины океана. Проведено строгое доказательство корректности всей процедуры построения разложения. Предложенный метод может быть использован при исследовании других систем, описываемых уравнением Гамильтона-Якоби с малым параметром. Доказанная теорема о разрешимости неоднородной задачи Штурма-Лиувилля полезна при исследовании других проблем, в которых подобная задача возникает.
Получены новые аналитические решения двумерных и трехмерных уравнений Стокса и неразрывности с линейно и экспоненциально меняющимися вязкостью и плотностью. Данные решения представляют интерес для теоретической геофизики.
На практике построенные асимптотические решения типа квазифотона или волнового вала можно использовать для усовершенствования вычислительных схем. Во-первых, эти решения можно непосредственно встроить в структуру алгоритма. Это позволит сделать численный алгоритм быстрее. Или надежнее, что иногда даже более важно. Во-вторых, данные решения описывают характер поведения волны. А значит численная схема при работе должна давать такой же характер поведения. Это позволяет проверить корректность работы численного алгоритма.
В последней части работы не только найдены аналитически эталонные решения, но и продемонстрирован вариант их использования на примере конкретного численного алгоритма. Данные решения полезны при создании эффективных вычислительных алгоритмов решения геофизических задач, так как они дают возможность для тестирования вычислительных схем и выбора наиболее приемлемой для данного класса проблем. Отметим, что эталонные решения для экспоненциально меняющейся вязкости дают аппарат для исследования численных схем, предназначенных для расчета течений в случае резких контрастов вязкости, что представляет наибольшую трудность для вычислений.
Методы исследования. В работе использованы асимптотические методы теории волн, в частности, пространственно-временной лучевой метод. При построении и анализе асимптотического разложения решения использована теория уравнения Гамильтона-Якоби, методы теории операторов, теория аналитических функций, теория канонических систем, векторный анализ, теория формальных асимптотических рядов, теория уравнений в частных производных. При построении эталонных решений уравнений Стокса и неразрывности использованы методы теории уравнений в частных производных, в частности, метод характеристик, а также методы теории обыкновенных дифференциальных уравнений. Численное решение уравнений Стокса и неразрывности проводилось методами конечных разностей с использованием MATLAB.
Научная новизна. Все результаты, представленные в диссертации, являются новыми. Решения типа квазифотона ранее были построены только в задачах акустики и электродинамики, решение же типа волнового вала ранее вообще никем не строилось. Некоторые
эталонные решения уравнений Стокса и неразрывности в двумерном случае ранее были построены, но в совершенно другой ситуации. В трехмерном случае для переменной вязкости никаких эталонных решений ранее вообще не предлагалось. Построенные в работе точные решения для линейно и экспоненциально меняющейся вязкости получены впервые.
Степень достоверности и апробация результатов работы. Достоверность результатов обеспечивается использованием строгих математических методов и полными доказательствами корректности всех построений.
Результаты диссертационного исследования были представлены на 7 международных научных конференциях:
-
А.И. Попов, В.М. Бабич. Гравитационные волны квазифотонного типа на поверхности жидкости // Международная конференция "Дифференциальные уравнения и смежные вопросы"(23-я сессия), посвященная И.Г.Петровскому. Москва, МГУ. 29.05.2011 - 04.06.2011
-
A.I. Popov, V.M. Babich. Asymptotic solutions of equations describing waves spreading on the surface of heavy liquid concentrated near moving points // 9th International Conference of Numerical Analysis and Applied Mathematics. ICNAAM 2011. G-Hotels, Halkidiki, Greece.
19.09.2011 - 25.09.2011
-
A.I. Popov, V.M. Babich. Wave wall type solution for liquid surface waves // 10th International Conference of Numerical Analysis and Applied Mathematics. ICNAAM 2012. Kypriotis Hotels and Conference Center, Kos, Greece. 19.09.2012 - 25.09.2012
-
A.I. Popov, V.M. Babich. High Frequency Wave Packet Concentrated near the Moving Line on the Liquid Surface // International Conference on Advancement in Science and Technology 2012. "Contemporary mathematics, Mathematical Physics and their Applications". IIUM iCAST 2012. Kulliyyah of Science, International Islamic University Malaysia. 07.11.2012 - 10.11.2012
-
A.I. Popov. Benchmark solutions for Stokes flow algorithm testing // Mathematical Challenge of Quantum Transport in Nanosystems. NRU ITMO, Saint Petersburg, Russia. 12.03.2013 -15.03.2013
-
A.I. Popov. New approach for Stokes flow algorithm testing // International Conference on Computational Science. Barcelona, Spain. 05.06.2013 - 07.06.2013
-
A.I. Popov, I.S. Lobanov, I.Yu. Popov, T.V. Gerya. Stokes flow computing algorithm based on Woodbury formula // Mathematical Challenge of Quantum Transport in Nanosystems. ITMO University, Saint Petersburg, Russia. 23.09.2014 - 26.09.2014
Исследования частично проводились по темам, поддержанных грантами:
-
Методы математической физики в квантовой механике, квантовой теории поля и теории распространения волн (Санкт-Петербургский Государственный Университет, Физический факультет, 01.01.2010 - 31.12.2012)
-
Научный проект "Высокочастотный волновой пакет, сосредоточенный в окрестности движущейся линии на поверхности жидкости"для представления на научном мероприятии "International Conference on Advancement in Science and Technology 2012"(Грант РФФИ,
07.11.2012 - 31.12.2012)
3) Развитие математических методов для задач квантовой механики, нанофизики и тео-
рии распространения волн (Санкт-Петербургский Государственный Университет, 01.01.2012 - 31.12.2012)
4) Разработка математических методов исследования сложных физических систем (Министерство образования и науки РФ, 10.09.2013 - 15.12.2014)
Положения, выносимые на защиту:
-
Доказательство существования решения типа квазифотона, сосредоточенного в окрестности пространственно-временного луча для волн на поверхности океана переменной медленно меняющейся глубины и алгоритм построения его полного асимптотического разложения, обоснование алгоритма с доказательством реализуемости всех его шагов;
-
Доказательство существования решения типа волнового вала, сосредоточенного в окрестности поверхности, сотканной из пространственно-временных лучей, для волн на поверхности океана переменной медленно меняющейся глубины и алгоритм построения его полного асимптотического разложения, доказательство разрешимости задач, возникающих на всех этапах схемы;
-
Построенные в аналитическом виде с использованием метода характеристик эталонные решения уравнений Стокса и неразрывности с линейно и экспоненциально зависящими от координат вязкостью и плотностью в двумерном и трехмерном случаях, реализация методики проверки расчета стоксовых течений численными методами. Обоснование метода тестирования на базе построенный явных частных решений с доказательством разрешимости всех возникающих в процессе построения задач;
-
Методика тестирования вычислительного алгоритма решения краевой задачи для уравнений Стокса и неразрывности с помощью сравнения результатов расчета численными методами с построенными явно эталонными решениями. Реализация метода тестирования на основе построенных эталонных решений для конечно-разностного многосеточного алгоритма расчета течений Стокса с переменными вязкостью и плотностью.
Структура работы.
Пространственно-временной лучевой метод для волн на поверхности тяжелой
В статье [51] выведены аналитические интегральные соотношения между топографией границы и температурой как функцией волнового числа для слоев жидкости, вязкость которых меняется экспоненциально с глубиной. Аналогичные соотношения между гравитационными силами и температурой также выведены. Показано, что если вязкость меняется с глубиной не более чем на порядок, хорошо работает приближение постоянной вязкости. При вязкости, меняющейся на несколько порядков, сильно меняются топографические и гравитационные ядра выведенных интегральных уравнений. При очень сильном изменении вязкости можно использовать приближение твердой нижней границы.
В работе [22] проведено сравнение различных эталонных решений для устойчивой конвекции с постоянной вязкостью, конвекции с переменной вязкостью и зависящей от времени конвекции с внутренним разогревом. Авторы сравнивают числа Нуссельта, скорости, температуры, потоки тепла, топографию. Тестируются программы, использующие конечно-разностные, конечно-элементные и спектральные вычислительные схемы.
В статье [59] проведено сравнение различных методов исследования термохимической конвекции для течения Буссинеска с бесконечным числом Прандтля, в частности, сравниваются методы отслеживания химической неоднородности, цепочки маркеров и полевой метод. В случае изотермической неустойчивости Рэлея-Тейлора получено хорошее соответствие для начального роста неустойчивого нижнего слоя, однако, развитие во времени и локализация более поздних нестабильностей меньшего масштаба могут различаться в разных подходах. Для моделирования термической конвекции плотного слоя у дна мантии найдено хорошее соответствие для времени нескольких переворачиваний, но потом разница результатов может быть значительной. Эффекты хаотического перемешивания могут привести к отсутствию детального соответствия.
В работе [53] методом Мусхелишвили получено в замкнутой форме аналитическое решение для изолированного эллиптического включения в сдвиговом течении. Решения подходят для несжимаемых полностью упругих или полностью вязких систем. Ограничений на контрасты свойств материалов нет. Решение позволяет найти кинематические (функция тока, скорость) и динамические (давление, напряжение) параметры. Авторы подробно не выводят решение, а сосредотачиваются на его использовании. Авторы получают новые результаты для медленного вязкого течения с включением. Другое важное приложение - эталонные решения для тестирования кодов. Это особенно интересно, так как в модели допустимы высокие контрасты параметров. Все решения представлены также в MATLAB.
В статье [50] изучена неустойчивость Рэлея-Тейлора в геофизических задачах. В работе [41] исследовано влияние упругости на неустойчивость Рэлея-Тейлора в системе из вязко-упругого слоя, покрывающего вязкий слой меньшей плотности. Хотя части литосферы на определенных временных масштабах можно считать упругими, это обычно игнорируется в геофизических моделях динамики мантии. Недавно было показано, что упругие и вязко-упругие свойства могут сильно влиять на литосферные сдвиговые зоны. Задача в статье рассмотрена как численно, так и аналитически. Показано, что упругость при определенном соотношении параметров приводит к ускорению неустойчивости Рэлея-Тейлора. Это происходит из-за уменьшения доли вязкости в деформации вязко-упругого слоя, что наблю-дается при временных масштабах, меньших максвелловского времени релаксации їм = Q где г/ - вязкость, а G- упругий модуль сдвига. Для тектонических плит Земли, в большинстве случаев данный эффект не наблюдается, однако встречаются ситуации, когда результаты существенно отличаются от чисто вязкого случая. Это подтверждено численными расчетами.
В статье [60] разработаны эталонные решения для сравнения с результатами численного моделирования динамики и термической структуры зон субдукции. Эталонный случай основан на термомеханическом подходе, сочетающем кинематические и динамические методы. Предложены способы изучения влияния граничных условий и реологии на краевые течения и итоговую термическую структуру. Сравнение с численными подходами показывает, что в них необходим точный учет скачков скорости вдоль зон сейсмичности и высокое разрешение для термического пограничного слоя. Небольшая модификация граничных условий для краевых течений также необходима, чтобы избежать сингулярности давления, которая есть в угловых точках у аналитического решения.
В работе [52] предложены эталонные решения для проверки численных алгоритмов расчета динамики литосферы с субдукцией для различных вычислительных кодов (эйлеровых и лагранжевых, конечно-элементных и конечно-разностных с маркерами и без маркеров). В эталонных тестах для моделирования свободной поверхности используется слой с малой вязкостью и нулевой плотностью. Рассмотрена и чисто свободная граница. Показано, что сходимость в расчетах поведения зон субдукции при увеличении разрешения сильно зависит от схемы усреднения вязкости вблизи движущейся реологической границы. Выбор гармонического среднего приводит к наибыстрейшей субдукции, арифметического среднего - к медленной субдукции, а геометрическое среднее дает промежуточный результат. Проверен еще ряд норм средних и выбрана та, которая дает лучшее поведение при увеличении разрешения. Сравнение слабого верхнего слоя и чисто свободной поверхности показывает, что соответствие наблюдается только при высоком разрешении. Сравнение численной модели и лабораторного эксперимента показало, что численный лагранжев подход воспроизводит экспериментальный результат только при достаточно высоком разрешении. Коды, основанные на эйлеровом подходе со слабым верхним слоем воспроизводят экспериментальные результаты если используется гармоническое или геометрическое среднее. Таким образом, использование слабого верхнего слоя для моделирования свободной поверхности дает удовлетворительные результаты только при выборе высокого разрешения и подходящего способа усреднения.
В работе [24] предложены эталонные решения для сравнения с численными методами исследования конвекции с высоким числом Прандтля в трехмерной декартовой геометрии. Проведено сравнение численных результатов с результатами лабораторных экспериментов. Рассмотрены случаи бимодальной конвекции с постоянной вязкостью и конвекции в прямоугольной ячейке с вязкостью, зависящей от температуры. Сравниваются числа Нуссельта, скорости, температуры, потоки тепла. Рассмотрено зависящее от времени поведение и предел устойчивости для бимодальной конвекции. Протестированы восемь программ и проведено их сравнение. Также проведено сравнение с результатами лабораторных экспериментов, которое показало удовлетворительное соответствие в пределах погрешности эксперимента.
В статье [63] предложены эталонные решения для проверки численных методов расчета конвекции мантии в сферической геометрии и протестирован конечно-элементный код CitcomS для трехмерной сферической конвекции. Представлены два класса модельных вычислений: стоксовы течения, а также термическая и термохимическая конвекции. Для сток-сова течения вычислены скорость течения, топография на границе ядро-мантия при различных степенях сферических гармоник с помощью кода CitcomS и проведено сравнение с решениями, найденными методом матрицы пропагатора. Для термической и термохимической конвекции рассмотрены 24 случая, отвечающие различным параметрам модели: числа Рэлея меняются в пределах от 7103 до 105, контраст вязкости меняется от 1 до 107. Для каждого случая найдены усредненные по времени характеристики, включая числа Нуссельта для границы ядро-мантия, скорости, средние температуры, максимальные и минимальные скорости течения, температура на средней глубине мантии и среднее отклонение. Для термохимической конвекции, в дополнение к термической конвекции, была проведена оценка степени рассеяния первоначально плотных компонент, а также относительное изменение их объемов. Для девяти случаев с малым контрастом вязкости было проведено сравнение кода CitcomS с ранее опубликованными результатами. Для чисел Нуссельта и скоростей результаты отличаются не более, чем на 1 процент. Остальные пятнадцать случаев соответствуют либо большому контрасту вязкости в течение Стокса, либо анализу термохимической конвекции. Ранее подобные результаты получены не были. Эти результаты могут быть использованы для тестирования других кодов. Также проведена оценка эффективности параллельных вычислений для кода CitcomS на суперкомпьютере Ranger (3072 ядра) в Техасском компьютерном центре.
Комплексный вариант лучевого метода и квазифотоны
Пусть в ПВ-лучевом разложении 1т9 0, причем знак равенства имеет место лишь на некотором ПВ-луче /. Тогда при є С 1 частные суммы ПВ-лучевого ряда были бы существенно отличны от нуля лишь в малой окрестности /, т.к. е = Є Е1тв. “Волна” Фм представляла бы собой волновой пакет - квазифотон, движущийся вдоль по-верхности жидкости с групповой скоростью v = г, г . Строить комплексные решения уравнения эйконала сложно в силу того, что теория характеристик - вещественная теория. Однако, некоторый аналог техники ПВ-лучевых разложений здесь применить возможно, если 9 и Ф считать не функциями, а формальными степенными рядами. Этот подход уже применялся в работах ([3] - глава 3, [1] и [7]). Отличие похода, применяемого в настоящей статье от подхода в [1], [7] - отказ от использования так называемых квазисопряженных решений, что позволяет несколько упростить процедуру построения квазифотонов. Квазифотон у нас - это разложение вида (2.17):
Предполагается, что г = г(т), і = 1,2 - это некоторый фиксированный ПВ-луч. Именно в окрестности этого ПВ-луча мы построим квазифотон. Аналог переменных rf впервые был введен А.П. Качаловым [14] в другой задаче. Мы будем предполагать, что ПВ-луч - “опорный ПВ-луч квазифотона” проходит через начало координат т = 0, г = 0, г = 1, 2. Введем новые координаты s,ai,a2. Здесь а\, й2 - декартовы координаты пересечения ПВ-луча, проходящего через заданную точку, с плоскостью г = 0. s = т. В координатах s,ai,a2 точке г = 0, г = 0 соответствует точка s = 0, а\ = й2 = 0. При замене переменных т,г]1,г]2 на s,a\,a2 9 и Ф.,- превращаются в ряды по степеням а\,а2. Подставляя разложение (2.70) в уравнение и краевые условия, мы придем к тем же уравнениям (2.19)–(2.24). Прежде всего займемся уравнением для в = в -\-в 1 + , которому можно придать вид уравнения Гамильтона-Якоби:
Наши построения здесь практически совпадают с построением в монографии [3] (глава 3). Для того, чтобы разложение (2.70) представляло квазифотон, мы предположим, что в и 9 1 вещественны, а квадратичная форма Ітв 2 положительно определенная (т.к. квазифотонное решение должно экспоненциально убывать при удалении от ПВ-луча, а если бы 9 1 было комплексным, то егв в каком-то из направлений обязательно возрастало бы). Экспоненциальное убывание решения по всем направлениям обеспечивает в 2 .
Следовательно, в3 находится квадратурой: в3 = — Г ( т т ) дьт. Подставляем в3 в (2.77) и находим в( . в 1 также находится ибо в 1 = в3г]3. 2.4.3 Второй порядок Найдем теперь вгК Удовлетворение равенства (2.76) во втором порядке эквивалентно соотношению: (2.82) ат 2 овлоиргп і 2 at оил + о/! ол,- 9llrfrf + \ тттттттт ) rfrP = О 2 оилос г, 2 at atJ г, ї 0 О Введем Г = (Г,-,), Г,-,- = Г,-,-, матрицу квадратичной формы 9 2 = T rfrP. Гі7- = 2 + \ ( ) eildmj + \ (А) lj + h (тяРш) 8й + (шш) = 0 (2.83) ат 2 oU ioUtm 2 at Oudi 2 oUdiot,-. 2 at at,J п U U и и Тогда (2.82) можно записать как матричное уравнение Риккати: dT ппТ ҐУ — = ГАГ + ВГ + Г + 6. (2.84) ат где А = - ( Ш—) ,В = - (ЖёА ,С=- (Ж Л . овлоиргп atrou i rj at at- п U и и Положим Г = ZY . (2.85) Пусть Z и У - решения следующей системы: -г = BZ + СУ от -г- = —AZ — BTY ат Тогда Г удовлетворяет уравнению (2.84). Проверим это. Продифференцируем тождество Y_1Y = І по г: —і—Y + Y l-г- = 0, отсюда находим —,—: ат ат ат (2.86) = —У —У . (2.87) ат ат ;г = У-1 + = BZY l + С + У-М У-1 + У_1БТ = ат ат ат = ГАГ + ВГ + Г Т + С Потребуем, чтобы матрицы Z и У удовлетворяли двум соотношениям: (2.88) У Z — Z У = 0, (2.89) Z — Z Y = її (2.90) (/ - единичная матрица, - знак эрмитова сопряжения). Условия (2.89) и (2.90) обеспечивают симметричность Г и положительную определенность ІтпГ. Кроме того, равенство (2.90) обеспечивает существование У-1, то есть неравенство dety ф 0. Доказывается от противного. Действительно, если dety = 0 при каком-нибудь т, то найдется вектор D = - (У Z — Z Y) D:
Рассмотрим малую окрестность ПВ-луча: г = г(т) + rf, в І = ві(т) + pi, где rf и pi- малые величины. Раскладывая правые части (2.93) в ряд Тейлора по rf и pi, получим в первом порядке уравнения в вариациях: 2l = ( Q2H А пз і ( э2и \ dr двід r\ odiddjп J / 9 \ / 9\ (2.94) dpj / 9 Н \ 7 / 9 Н \ _ л = — atiati ) П [ а» ад Pi ат ос otJ r.1 of оді г\ J В окрестности ПВ-луча перейдем к координатам s, а\, а . Разложим if и рі в ряд по а\, аъ: дг]г дг]г г/ =—а\ +—ci2 + (2.95) да\ да2 dpi dpi Pi =—а\ +—ci2 + (2.96) да\ да2 Подставляем rf и рі в систему (2.94) и приравниваем коэффициенты при одинаковых степенях аі. Тогда в первом порядке получим: d I drf \ I д2И \ дг\і , / д2И \ 0Pj dr \ oak J oUioc,-1 п oak odiOUj п оа / ) /„ \ / о \ „ / о \ а 2.97 d I dpi \ I д2И \ on3 I д2И \ opj v = — г г 1 — г dr оаи огоэ дай огові дай Если ввести матрицы Кь = тр, Zn = -к- = т , то система (2.97) приобретает вид системы да bt да да (2.86). Итак, показали, что решения системы (2.86) можно записать в виде: дг]гдві = і\— ; % = т — (2.98) dak ак 2.4.4 Расчет высших приближений Уравнения для оказываются линейными. Приравнивая к 0 члены г-го порядка в равенстве (2.76), получаем: QQvr) / 92Н \ -двуг ( 92Н \ двут дву2 Ь Г]: Ь 77 7Г ТГ- 7Г ОТ Ос1о9 :з п дг}3 OVpiOVpj „ ОТ)1 drf где Б,г- известная функция, содержащая элементы, определенные на предыдущих этапах. С учетом обозначений (4.16) уравнение (2.99) записывается в виде: — V0 ,AV6 r — г/, B Ve = Sr (2.100) дт Здесь скобки - символ ск алярного произведения без комплексного сопряжения ( а, Ъ\ = Cljbj. Знаем, что 9у2 = (Гп,п) = (ZY ln,ri). Определим 7ву2 : \V 2 ) = Q— (ZY 1)-irfrf = (ZY l)-irf5ik + (ZY 1)-irf6ik = Vk (2.101) = 2 (ZY-1),. ril = 2 (ZY-lri), Принимая во внимание (2.101) и симметричность матрицы А, левую часть равенства (2.100) можно переписать в виде:
Перейдем к уравнениям переноса. В работах [1], [7] интегрирование уравнений переноса в классе формальных степенных рядов основывалось на понятии квазисопряженных решений. Здесь мы будем пользоваться более простой техникой, еще более близкой к технике стандартного лучевого метода. Положения анализа, использующиеся для решения рекуррентных уравнений справедливы, если функции заменить формальными степенными разложениями. все построения параграфов 3.2-3.3 проходят фактически без изменений. На этом главного члена получаем формулу:
Построение Фо
Здесь в Sra собраны все члены с б1-7, j п. При гг 2 это линейное уравнение. Решив уравнение (3.62), мы найдем вп. Таким образом, все коэффициенты формального степенного ряда для в определены. Ряд для в построен.
Оценим размер области, где волновой пакет существенно отличен от 0. Напомним, что волновой пакет Ф имеет вид (3.8): , rf у/є. Координаты (г, г]1,!]2) получаются из (т, 1, 2) линейным преобразованием, поэтому в координатах (т, 1, 2) размер области имеет тот же порядок. Итак, размер области, где волновой пакет существенно отличен от 0, имеет порядок О (у/є).
Все Ф - ищем в виде формальных степенных рядов по г]2. Рассмотрим нулевой порядок. Задача для Фо имеет вид: -к-20 — к2Фп = О, дФ0 dz = 0. z=—H (3.68) Здесь к2 - это формальный степенной ряд по rf (формула (3.12)). Естеcтвенно, что для Фо имеет место формула: Фо = А)(т, г/ ,г] ) cosh(k(z + Н)). (3.69) Формальный степенной ряд Ао будем искать из условия разрешимости задачи для Фі. Нам потребуется известная теорема о разрешимости неоднородной задачи Штурма-Лиувилля. Теорема сохраняет свой вид в случае, если решение ищется в случае формального степенного ряда.
Теорема 1. Рассмотрим однородную задачу Штурма-Лиувилля: . (3.71) Х/, f3j, j = 0,1 - вещественные числа. p(x), q(x) - формальные степенные ряды (ф.с.р.). Пусть существует решение в виде ф.с.р. уо ф 0. Тогда необходимое и достаточное условие существования решения в виде ф.с.р неоднородной задачи Штурма-Лиувилля
Введем вектор А: (є), г ЙП (є), г ЙГ2 (є) (3.80) Здесь vgr - групповая скорость. С учетом вышесказанного, уравнение (3.77) примет вид: div А = 0. (3.81) Введем некоторые понятия, необходимые для формулировки теоремы о неявной функции для формальных степенных рядов.
Предварительные сведения. Рассмотрим г + s формальных степенных рядов. оо У А3а{х)уа = а!-3 , j = 1,2, , г. (3.82) оо У А3а(х)уа = Ь 3 г , j = r + l,r + 2, , r + s. (3.83) H=i Q! = (Q!1,Q!2, ,Cis), (3.84) \a\ = a\ + «г + + as, cti Є Z+. (3.85) Заданные коэффициенты A3a є C c в некоторой области Q clr. Переменные x = (xl ) будем называть конечными переменными, у = (у1, у2, , ys) - малыми переменными. Пусть х и у - формальные степенные ряды: оо х = У В3а{а)Ъа, j = 1,2, , г, (3.86) оо у = У В3а{а)Ъа, j = r + l,r + 2, , r + s. (3.87) H=i ВЦа) Є С(Пі) (Пі С Mf) - функции, причем область значений B3nn п\(а) С Q. а = {а 1 ,о! 2 , ,а) - конечные переменные. Ъ= (\Р ,Ъ 2 , ,b s ) - малые переменные. Тогда х3 и у3, даваемые формулами (3.86), (3.87) можно подставлять в ряды (3.82), (3.83). При подстановке формальных степенных рядов для х3 и у3 функции А3а{х) заменяются рядами Тейлора-Маклорена. Сходимость не предполагается.
Определение. Обращение рядов (3.82), (3.83) - это нахождение таких рядов (3.86), (3.87), что при их подстановке в (3.82), (3.83) получится: оо У А3а(х)уа = a}1 , j = 1,2, , г, оо У А3а{х)уа = Ь 3 г , j = r + l,r + 2, , r + s. H=i Теорема 2. Рассмотрим ряды: У А3а(х)уа = a}1 , j = 1,2, , г, \а\=0 оо У А3а(х)уа = Ь 3 г , j = r + l,r + 2, ,r + s H=i (3.88) (3.89) Пусть определитель D( A1 (r) Ar Ar+l 7/1) /T+s 7/(s) ) J = D(xW,....,x(r\yW,....,yW) X(J =XQ , у )=0 Ф 0, (3.90) где дСг —— ОС (3.91) Тогда можно обратить ряды (3.88) и (3.89). От координат т, 1, 2 перейдем к координатам Sjd1, , которые вводятся следующим образом (см. Рис. 4): Рис. 4: Новые координаты а1 - длина дуги опорной кривой между точками Р и Р\. а2 - длина дуги между точками Р\ и М\ (длина проекции ПВ луча на плоскость г = 0).
Для того, чтобы перейти к новым координатам s a1 2 нужно определить ПВ луч. Его можно найти из решения канонической системы уравнений: d % 9Н d9T двы дЛ d9 дН т = s, — = ,— = 0, = — —т, — = — Н + всі as oO i as as at,1 as дв і Первые два из равенств, входящих в (3.92), приводят к соотношению: d % duj 9ri dr dki (3.92) (3.93) Вектор vgr = (vgri, v9r2) с компонентами т - - это по определению вектор групповой скорости. Начальные данные для системы (3.92) в виде ф.с.р. по а2 задаются на поверхности г = 0:
Линейно меняющаяся вязкость. Трехмерный случай
В предыдущих параграфах были найдены частные решения системы уравнений Стокса и неразрывности для линейно и экспоненциально меняющейся вязкости в двумерном и трехмерном случаях. Опишем схему проверки алгоритма. Рассмотрим прямоугольную область. Посчитаем значения, которые принимают эталонные решения (vx, vy и Р) на границе области. Теперь возьмем эти значения в качестве граничных условий. По теореме единственности, решение граничной задачи должно совпасть с нашим эталонным решением.
Итак, мы нашли точное аналитическое решение конкретной граничной задачи. Теперь решим ту же самую задачу с помощью численного алгоритма. Сравнение результатов позволит оценить качество работы алгоритма.
Воспользуемся данным методом для проверки алгоритма Variable-viscosity3Dests2 (http://www.solid-earth.net/5/issue1.html - Supplement к статье [33]). Рассматриваем прямоугольную область 0 х х Size, 0 у у size- Считаем, что f(rj) = А = const and г/ = ах + by + с. Точные аналитические (эталонные) решения будем обозначать за vXA,vyA, численные решения - за УХ}П,УУ}П.
Замечание. На всех рисунках "гг"означает "численное решение "а"означает "аналитическое решение"(эталонное).
Численный алгоритм основан на методе конечных разностей. В качестве граничных условий в программе устанавливаем следующие значения параметров (условия на границе прямоугольника): дП = {х = О, X = Xgize, У = 0, у = У size}
Ошибка в работе алгоритма определяется величиной отклонения численного решения от эталонного. Также мы оцениваем относительные ошибки трех типов: L , L\, L2 для различных контрастов вязкости. Результаты расчетов показывают, что алгоритм демонстрирует хорошую сходимость для малых контрастов вязкости, однако дает сбои при высоких контрастах. На Рис. 6 изображены численные и аналитические решения системы уравнений Стокса и неразрывности и их разность для распределения вязкости и плотности, изображенных на Рис. 7. На Рис. 8 показан график зависимости относительной ошибки от шага сетки в логарифмическом масштабе. Положительный наклон графика означает сходимость алгоритма (т.е. при уменьшении шага сетки точность решения увеличивается). Рис. 5: Прямоугольная область, для которой ведутся расчеты
Зависимость вязкости г/ от координат задаем, определяя вязкость в трех вершинах прямоугольника. 771 = 7/(0,0), 7/2 = 7/(0, ysize), V3 = V(xsize, 0)
Для случая линейно меняющейся вязкости были установлены следующие значения параметров: %size = У««е = 1) Сж = 0, Gy = 10, 7/i = 1, /Зі = 1, /?2 = 3 X 10 , С = 7/1, а = (т/з Vl)/Xsize, Ь = (7/2 — Ці) / У size) p = fi\ {ax + by + с) + /З2. Рисунки 6, 7, 8 соответствуют малым контрастам вязкости, рисунки 9, 10, 11 - большим контрастам вязкости. Рис. 7: Распределение вязкости и плотности ; двумерный случай, линейно меняющаяся вязкость, малый контраст вязкости (2 = 3 = 5). Рис. 8: Логарифм относительной ошибки vs Логарифм шага сетки; двумерный случай, линейно меняющаяся вязкость, малый контраст вязкости (г]2 = Щ = 5); синяя линяя - давление, зеленая - vx, черная - vy; обычная линия - ошибка типа L\, пунктирная линия - ошибка типа LQO, линия с точками - ошибка типа Z-2.
Распределение vx, vy и Р; двумерный случай, линейно меняющаяся вязкость, большой контраст вязкости (щ = Щ = 100). Рис. 10: Распределение вязкости и плотности ; двумерный случай, линейно меняющаяся вязкость, большой контраст вязкости (2 = 3 = 100).
Логарифм относительной ошибки vs Логарифм шага сетки; двумерный случай, линейно меняющаяся вязкость, большой контраст вязкости (щ = Щ = ЮО); синяя линяя -давление, зеленая - vx, черная - vy; обычная линия - ошибка типа L\, пунктирная линия -ошибка типа L00, линия с точками - ошибка типа L2. Рис. 12: Распределение vx, vy и Р; двумерный случай, экспоненциально меняющаяся вязкость, малый контраст вязкости (щ = Щ = 5).
Распределение вязкости г/ и плотности р; двумерный случай, экспоненциально меняющаяся вязкость, малый контраст вязкости (щ = Щ = 5).
Рисунки 12, 13, 14 соответствуют малым контрастам вязкости, рисунки 15, 16, 17- большим контрастам вязкости. Рис. 14: Логарифм относительной ошибки vs Логарифм шага сетки; двумерный случай, экспоненциально меняющаяся вязкость, малый контраст вязкости (щ = Щ = 5); синяя линяя - давление, зеленая - vx, черная - vy; обычная линия - ошибка типа L\, пунктирная линия -ошибка типа L , линия с точками - ошибка типа Z-2.
Распределение vx, vy и Р; двумерный случай, экспоненциально меняющаяся вязкость, большой контраст вязкости (щ = Щ = 100).
Распределение вязкости и плотности ; двумерный случай, экспоненциально меняющаяся вязкость, большой контраст вязкости (2 = 3 = 100). Рис. 17: Логарифм относительной ошибки vs Логарифм шага сетки; двумерный случай, экспоненциально меняющаяся вязкость, большой контраст вязкости (г]2 = Щ = 100); синяя линяя - давление, зеленая - vx, черная - vy; обычная линия - ошибка типа L\, пунктирная линия - ошибка типа L , линия с точками - ошибка типа L2.