Решения, сосредоточенные в окрестности пространственно-временных лучей, и эталонные решения в задачах о движении жидкости тема диссертации и автореферата по ВАК РФ 01.04.02, кандидат наук Попов, Антон Игоревич

• Попов, Антон Игоревич
• кандидат наук
• 2016

Оглавление диссертации кандидат наук Попов, Антон Игоревич

2.1 Энергетическое соотношение............................................. 17

2.1.1 Постановка задачи................................................ 17

2.1.2 Энергетическое соотношение....................................... 19

2.1.3 Вектор Умова У как “вектор работы”............................... 21

2.2 Пространственно-временной лучевой метод для волн на поверхности тяжелой

жидкости............................................................... 22

2.3 Уравнения переноса .................................................... 27

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.1 Введение .............................................................. 37

3.2 Исходные уравнения..................................................... 38

3.3 Пространственно-временной лучевой метод ............................... 38

3.4 Построение Ө........................................................... 40

3.5 Построение Ф^.......................................................... 45

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

1

4.2.3 Экспоненциально меняющаяся вязкость. Трехмерный случай........ 63

4.3 Схема метода и проверка корректности работы конкретного алгоритма .... 65

4.3.1 Двумерный случай. Пример для линейно меняющейся вязкости .... 66

4.3.2 Двумерный случай. Пример для экспоненциально меняющейся вязкости 69

4.3.3 Трехмерный случай. Пример для линейно меняющейся вязкости .... 71

4.3.4 Трехмерный случай. Пример для экспоненциально меняющейся вязкости 75

5 Заключение 79

Список литературы 81

6 Приложения 87

2

Введение

Актуальность. Проблема описания движений жидкости возникает во многих физических задачах. Большой интерес представляют геофизические задачи, актуальность которых связана с необходимостью прогнозирования природных явлений и обеспечения функционирования различных технических систем. Данные задачи чрезвычайно сложны и обычно допускают только численное решение. В то же время, прямые расчеты без предварительного математического анализа часто оказываются недостаточно эффективными. Кроме того, возникают вопросы верификации численных результатов, проверки корректности методов, качества работы численных алгоритмов. Для этой цели могут использоваться аналитические результаты, полученные строгими математическими методами. Разработке именно таких подходов и посвящена данная диссертационная работа. Рассматриваются две задачи о движениях жидкости, порожденные геофизическими явлениями: распространение волн на поверхности океана и течения в мантии Земли.

В задаче о поверхности океана показано существование волн, локализованных в окрестности движущейся по поверхности точки или линии. Использован асимптотический подход. В основе анализа - пространственно-временной лучевой метод, предложенный в теории дифракции В.М.Бабичем, В.С.Булдыревым, И.А.Молотковым. Рассмотрена задача об океане с медленно меняющейся глубиной. Скорость изменения глубины взята в качестве малого параметра задачи, по которому строится асимптотическое разложение. Анализируются решения двух типов: квазифотон (волна, сосредоточенная в окрестности движущейся точки) и волновой вал (волна, сосредоточенная в окрестности движущейся линии). В обоих случаях построено полное асимптотическое разложение соответствующего решения. Проведено доказательство корректности процедуры на всех ее этапах. Данное решение, полученное аналитически в рамках строгого математического подхода, можно использовать в качестве эталонного в данной задаче.

В задаче о течениях в мантии Земли предложен другой подход к построению эталонных решений. Данные течения описываются уравнениями Стокса и неразрывности с переменными вязкостью и плотностью. Подобного типа уравнения возникают и в других областях физики, в частности, в задачах о течениях в наноструктурах. Для выбранных конкретных типов зависимостей (линейной и экспоненциальной) вязкости и плотности от координат найдены аналитически явные решения уравнений в двумерном и трехмерном случаях. Данные эталонные решения очень удобно использовать для тестирования вычислительных алгоритмов. Схема действий такова. Рассмотрим построенное аналитическое решение в некоторой области. Поставим на границе области условия, задаваемые данным решением и рассмотрим краевую задачу для уравнений Стокса и неразрывности с данными граничными условиями. По теореме единственности решение данной краевой задачи должно совпадать с исходным аналитическим эталонным решением. Для проверки алгоритма решим данную краевую задачу с помощью тестируемой вычислительной схемы. Отклонение результата от аналитического показывает напрямую погрешность вычислений. При этом имеется возможность достаточно просто исследовать характер сходимости, скорость работы алгоритма и т.п. при

3

различных параметрах системы или вариациях вычислительной схемы, выбирая наилучший вариант для данной задачи.

Объектом исследования являются геофизические явления, связанные с динамикой жидкости. Предмет исследования - математические модели, описывающие волны на поверхности океана переменной глубины и течения в мантии Земли.

Цель диссертационного исследования - исследование волн на поверхности океана переменной глубины, сосредоточенных в окрестности движущихся точек и линий, асимптотическими методами и построение необходимых для тестирования вычислительных схем и алгоритмов эталонных решений уравнений Стокса и неразрывности с переменными вязкостью и плотностью.

Для достижения этой цели в диссертации:

1. Разработан алгоритм построения полного асимптотического разложения решения типа квазифотона, сосредоточенного в окрестности пространственно-временного луча для волн на поверхности океана переменной медленно меняющейся глубины, проведено обоснование алгоритма с доказательством реализуемости всех его шагов;

2. Разработан алгоритм построения полного асимптотического разложения решения типа волнового вала, сосредоточенного в окрестности поверхности, сотканной из пространственновременных лучей, для волн на поверхности океана переменной медленно меняющейся глубины, доказана разрешимость задач, возникающих на всех этапах схемы;

3. Построены в аналитическом виде с использованием метода характеристик эталонные решения уравнений Стокса и неразрывности с линейно и экспоненциально зависящими от координат вязкостью и плотностью в двумерном и трехмерном случаях,

4. Реализована методика тестирования вычислительного алгоритма решения краевой задачи для уравнений Стокса и неразрывности с помощью сравнения результатов расчета численными методами с построенными явно эталонными решениями.

Теоретическая и практическая значимость работы.

В работе получены новые теоретические результаты, важные для построения математических моделей и теоретического описания геофизических явлений. А именно, показано существование нового типа волн на поверхности океана - волн, сосредоточенных в окрестности движущейся по этой поверхности точки (квазифотон) или линии (волновой вал) и экспоненциально убывающие при удалении от них. Для данных решений построены полные асимптотические разложения по малому параметру - скорости изменения глубины океана. Проведено строгое доказательство корректности всей процедуры построения разложения. Предложенный метод может быть использован при исследовании других систем, описываемых уравнением Гамильтона-Якоби с малым параметром. Доказанная теорема о разрешимости неоднородной задачи Штурма-Лиувилля полезна при исследовании других проблем, в которых подобная задача возникает.

Получены новые аналитические решения двумерных и трехмерных уравнений Стокса и неразрывности с линейно и экспоненциально меняющимися вязкостью и плотностью. Данные решения представляют интерес для теоретической геофизики.

4

На практике построенные асимптотические решения типа квазифотона или волнового вала можно использовать для усовершенствования вычислительных схем. Во-первых, эти решения можно непосредственно встроить в структуру алгоритма. Это позволит сделать численный алгоритм быстрее. Или надежнее, что иногда даже более важно. Во-вторых, данные решения описывают характер поведения волны. А значит численная схема при работе должна давать такой же характер поведения. Это позволяет проверить корректность работы численного алгоритма.

В последней части работы не только найдены аналитически эталонные решения, но и продемонстрирован вариант их использования на примере конкретного численного алгоритма. Данные решения полезны при создании эффективных вычислительных алгоритмов решения геофизических задач, так как они дают возможность для тестирования вычислительных схем и выбора наиболее приемлемой для данного класса проблем. Отметим, что эталонные решения для экспоненциально меняющейся вязкости дают аппарат для исследования численных схем, предназначенных для расчета течений в случае резких контрастов вязкости, что представляет наибольшую трудность для вычислений.

Методы исследования. В работе использованы асимптотические методы теории волн, в частности, пространственно-временной лучевой метод. При построении и анализе асимптотического разложения решения использована теория уравнения Гамильтона-Якоби, методы теории операторов, теория аналитических функций, теория канонических систем, векторный анализ, теория формальных асимптотических рядов, теория уравнений в частных производных. При построении эталонных решений уравнений Стокса и неразрывности использованы методы теории уравнений в частных производных, в частности, метод характеристик, а также методы теории обыкновенных дифференциальных уравнений. Численное решение уравнений Стокса и неразрывности проводилось методами конечных разностей с использованием MATLAB.

Научная новизна. Все результаты, представленные в диссертации, являются новыми. Решения типа квазифотона ранее были построены только в задачах акустики и электродинамики, решение же типа волнового вала ранее вообще никем не строилось. Некоторые эталонные решения уравнений Стокса и неразрывности в двумерном случае ранее были построены, но в совершенно другой ситуации. В трехмерном случае для переменной вязкости никаких эталонный решений ранее вообще не предлагалось. Построенные в работе точные решения для линейно и экспоненциально меняющейся вязкости получены впервые.

Степень достоверности и апробация результатов работы. Достоверность результатов обеспечивается использованием строгих математических методов и полными доказательствами корректности всех построений.

Результаты диссертационного исследования были представлены на 7 международных научных конференциях:

1) А.И. Попов, В.М. Бабич. Гравитационные волны квазифотонного типа на поверхности жидкости // Международная конференция "Дифференциальные уравнения и смежные вопросы"(23-я сессия), посвященная И.Г.Петровскому. Москва, МГУ. 29.05.2011 - 04.06.2011

5

2) 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

3) 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

4) 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

5) 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

6) A.I. Popov. New approach for Stokes flow algorithm testing // International Conference on Computational Science. Barcelona, Spain. 05.06.2013 - 07.06.2013

7) 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

Исследования частично проводились по темам, поддержанных грантами:

1) Методы математической физики в квантовой механике, квантовой теории поля и теории распространения волн (Санкт-Петербургский Государственный Университет, Физический факультет, 01.01.2010 - 31.12.2012)

2) Научный проект "Высокочастотный волновой пакет, сосредоточенный в окрестности движущейся линии на поверхности жидкости"для представления на научном мероприятии "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)

Публикации.

По теме диссертации автором опубликовано 16 научных работ, из них 6 статей в журналах из Перечня ВАК, 4 - в журналах из списков Web of Science/Scopus.

Статьи, входящие в Перечень ВАК:

1. Popov, A.I. On the Stokes flow computation algorithm based on Woodbury formula [Text] / T. Gerya, I.S. Lobanov, A.I. Popov, I.Yu. Popov // Nanosystems: Physics, Chemistry, Mathematics.

- 2015. - Vol. 6. - P. 140-145. - DOI 10.17586/2220-8054-2015-6-1-140-145 (0,38 п.л. / 0,13 п.л.).

2. Попов, А.И. Исследование разрешимости задачи Штурма-Лиувилля при построении

6

асимптотического ряда [Текст] / А.И. Попов // Научно-технический вестник информационных технологий, механики и оптики. - 2015. - Т. 15. - No. 5. - С. 916-920. - ISSN 2226-1494. -DOI 10.17586/2226-1494-2015-15-5-916-920 (0,31 п.л. / 0,31 п.л.).

3. Popov, A.I. Benchmark solutions for nanoflows [Text] / T. Gerya, I.S. Lobanov, A.I. Popov,

I.Yu. Popov // Nanosystems: Physics, Chemistry, Mathematics. - 2014. - Vol. 5. - P. 391-399 (0,56 п.л. / 0,19 п.л.).

4. Попов, А.И. Волновые валы для волн на поверхности тяжелой жидкости [Текст] / А.И. Попов // Записки научных семинаров ПОМИ. - Санкт-Петебург, 2012. - Т. 409. - С. 151-175 (1,56 п.л. / 1,56 п.л.).

5. Попов, А.И. Асимптотическое решение уравнения Гамильтона-Якоби, сосредоточенное вблизи поверхности [Текст] / В.М. Бабич, А.И. Попов // Записки научных семинаров ПОМИ.

- Санкт-Петербург, 2011. - Т. 393. - С. 23-28 (0,38 п.л. / 0,19 п.л.).

6. Попов, А.И. Квазифотоны волн на поверхности тяжелой жидкости [Текст] / В.М. Бабич, А.И. Попов // Записки научных семинаров ПОМИ. - Санкт-Петербург, 2010. - Т. 379. -С. 5-23 (1,19 п.л. / 0,60 п.л.).

Статьи из списков Web of Science/Scopus:

1. Popov, A.I. Practical analytical solutions for benchmarking of 2-D and 3-D geodynamic Stokes problems with variable viscosity [Text] / T. Gerya, I.S. Lobanov, A.I. Popov, I.Yu. Popov,

S.I. Popov // Solid Earth. - 2014. - Vol. 5. - P. 461-476 (1,0 п.л. / 0,3 п.л.).

2. Popov, A.I. Numerical approach to the Stokes problem with high contrasts in viscosity [Text] / T. Gerya, I.S. Lobanov, A.I. Popov, I.Yu. Popov // Applied Mathematics and Computation. -2014. - Vol. 235. - P. 17-25 (0,56 п.л. / 0,19 п.л.).

3. Popov, A.I. High Frequency Wave Packet Concentrated near the moving Line on the Liquid Surface [Text] / A.I. Popov // Journal of Physics: Conference Series. - 2013. - Vol. 435. - DOI 10.1088/1742-6596/435/1/011001 (0,38 п.л. / 0,38 п.л.).

4. Popov, A.I. Wave wall type solution for liquid surface waves [Text] / A.I. Popov // AIP Conference Proceedings. - 2012. - Vol. 1479. - P. 728-731. - ISBN 978-0-7354-1089-3 (0,25 п.л. / 0,25 п.л.).

Прочие публикации:

1. Popov, A.I. Practical analytical solutions for benchmarking of 2-D and 3-D geodynamic Stokes problems with variable viscosity [Text] / T. Gerya, I.S. Lobanov, A.I. Popov, I.Yu. Popov,

S.I. Popov // Solid Earth Discussions. - 2013. - Vol. 5. - P. 2203-2281. - ISSN 1869-9537.

- DOI 10.5194/sed-5-2203-2013. - URL: http://www.solid-earth-discuss.net/5Z2203/2013/sed-5-2203-2013.html (4,94 п.л. / 1,65 п.л.).

2. Popov, A.I. New approach for Stokes flow algorithm testing [Text] / A.I. Popov // ICCS 2013, Abstract Book. - 2013 (0,06 п.л. / 0,06 п.л.).

3. Popov, A.I. Benchmark solutions for Stokes flow algorithm testing [Text] / A.I. Popov // MCQTN 2013, Abstract Book. - 2013. - P. 18-19. - 24 p. (0,13 п.л. / 0,13 п.л.).

4. Попов, А.И. Асимптотический анализ решений типа волновых пакетов, сосредоточенных в окрестности движущейся линии на поверхности жидкости [Текст] / А.И. Попов //

7

НИУ ИТМО. - Санкт-Петербург, 2012. - С. 262-263. - 270 с. (0,13 п.л. / 0,13 п.л.).

5. Popov, A.I. High Frequency Wave Packet Concentrated near the Moving Line on the Liquid Surface [Text] / A.I. Popov // iCAST 2012, Abstract Book. - 2012. - 200 p. (0,06 п.л. / 0,06 п.л.).

6. Попов, А.И. Гравитационные волны квазифотонного типа на поверхности жидкости [Текст] / В.М. Бабич, А.И. Попов // Международная конференция, посвященная 110-ой годовщине И. Г. Петровского: Тезисы докладов. - 2011. - С. 143-144. - 424 с. (0,13 п.л. / 0,13 п.л.).

Положения, выносимые на защиту:

1. Доказательство существования решения типа квазифотона, сосредоточенного в окрестности пространственно-временного луча для волн на поверхности океана переменной медленно меняющейся глубины и алгоритм построения его полного асимптотического разложения, обоснование алгоритма с доказательством реализуемости всех его шагов;

2. Доказательство существования решения типа волнового вала, сосредоточенного в окрестности поверхности, сотканной из пространственно-временных лучей, для волн на поверхности океана переменной медленно меняющейся глубины и алгоритм построения его полного асимптотического разложения, доказательство разрешимости задач, возникающих на всех этапах схемы;

3. Построенные в аналитическом виде с использованием метода характеристик эталонные решения уравнений Стокса и неразрывности с линейно и экспоненциально зависящими от координат вязкостью и плотностью в двумерном и трехмерном случаях, реализация методики проверки расчета стоксовых течений численными методами. Обоснование метода тестирования на базе построенный явных частных решений с доказательством разрешимости всех возникающих в процессе построения задач;

4. Методика тестирования вычислительного алгоритма решения краевой задачи для уравнений Стокса и неразрывности с помощью сравнения результатов расчета численными методами с построенными явно эталонными решениями. Реализация метода тестирования на основе построенных эталонных решений для конечно-разностного многосеточного алгоритма расчета течений Стокса с переменными вязкостью и плотностью.

1 Обзор имеющихся результатов

У большого количества гиперболических уравнений есть высокочастотные асимптотические решения, сосредоточенные в окрестности некоторой точки, движущейся вдоль луча (представляющего собой решение канонической системы) с групповой скоростью. В квантовомеханическом случае роль луча играет классическая траектория частицы. Широко известны сосредоточенные волны другого типа - солитоны. Они представляют собой решения нелинейных уравнений. В нашем же случае речь идет о высокочастотных линейных волновых пакетах. Исследование таких решений как для классического волнового уравнения, так и для квантово-механического уравнения Шредингера имеет довольно давнюю историю, беря начало с работ [4], [16]. Начиная с работы [8], подобные решения стали называть ква

8

зифотонами. Квазифотонные решения обычно строили очень сложными методами - либо методом В.П.Маслова, либо с помощью метода пограничного слоя в теории дифракции. В настоящей работе используется другой подход - пространственно-временной лучевой метод, родственный геометро-оптическому (или ВКБ) приближению. Построения такого типа проводились в [3], [9], [13], [14], [21]. Технически, процедура основана на использовании формальных степенных рядов (см. [10]). Схема построения аналогична получению классического пространственно-временного разложения. Однако, аналог эйконала в нашем случае комплексный, что закрывает возможность использования классического метода характеристик. Тем не менее, сходную процедуру реализовать удается.

В случае электромагнитного поля квазифотонные решения построены в работе [14]. Здесь их называют нестационарными (пространственно-временными) гауссовыми пучками. В указанной работе рассматривается данный тип решений системы уравнений Максвелла для среды, в которой электрическая и магнитная проницаемости представляют собой произвольные эрмитовы тензоры, зависящие только от координат. Подобные уравнения встречаются при описании распространения электромагнитных волн в кристаллах (с вещественными симметрическими тензорами электрической и магнитной проницаемости) или в ферромагнетиках. Заметим, что высокочастотный спектр таких решений сосредоточен вблизи фиксированной частоты.

В нашем случае рассматривается не волновое уравнение, а уравнение для волн на поверхности тяжелой жидкости, что вызывает необходимость серьезно модифицировать схему. Здесь исходное уравнение есть просто уравнение Лапласа для потенциала поля скоростей, а волновой характер движения появляется при рассмотрении граничного условия на свободной поверхности жидкости [3]. Малым параметром, по которому ведется разложение в асимптотическом ряде для решения, является скорость изменения глубины. Задачу удается свести к уравнению Гамильтона-Якоби, что открывает возможность использовать пространственновременной лучевой метод, который и позволяет решить задачу - построить полное асимптотическое разложение для решения типа квазифотона.

Следующая задача, рассмотренная в настоящей работе - это построение асимптотического решения типа волнового вала, т.е. решения, сосредоточенного в окрестности движущейся с групповой скоростью линии. Это естественное обобщение задачи о квазифотоне. Его построение требует существенно более сложной техники. Решений данного типа ранее построено не было ни для классического волнового уравнения, ни в квантово-механическом случае. В настоящей работе данное асимптотическое решение строится (как и квазифотон) для волн на поверхности жидкости.

В последней части работы найдены аналитически эталонные решения, а также продемонстрирован вариант их использования на примере конкретного численного алгоритма. Данные решения полезны при создании эффективных вычислительных алгоритмов решения геофизических задач, так как они дают возможность для тестирования вычислительных схем и выбора наиболее приемлемой для данного класса проблем.

В геофизике существует множество численных алгоритмов для расчета течений в мантии

9

Земли. В случае, когда мантия Земли неоднородна (существуют вкрапления высокой вязкости), некоторые алгоритмы начинают плохо работать. С точки зрения математики, задача сводится к решению уравнения Стокса с высокими контрастами вязкости.

В работе предложены методы проверки численных алгоритмов. Иногда одну и ту же задачу удается решить двумя разными способами. С одной стороны, предъявляется точное аналитическое решение задачи. С другой стороны, задача решается посредством применения численных алгоритмов. Сравнение результатов позволяет оценить качество работы численного метода. Нам удалось найти точные аналитические решения уравнений Стокса и неразрывности для линейно и экспоненциально меняющейся вязкости в двумерном и трехмерном случаях. Эти решения удобно использовать в качестве эталонных для проверки численных

методов.

Прежде чем говорить о построении эталонных решений, опишем современное состояние дел. Численное моделирование геофизических процессов - серьезная задача, для решения

которой используются мощные компьютеры в сочетании с продвинутыми вычислительными методами (см., например, [39] и ссылки в них). Одна из главных задач в настоящее время - решение системы уравнений Стокса и неразрывности в отсутствие трения в случае сильных контрастов вязкости. Жидкость считается несжимаемой. Проверка точности работы

численных алгоритмов с помощью эталонных решений широко используется в геофизике

(см., например, [22], [27], [28], [30], [31], [32], [44], [60]). В настоящее время аналитические и

численные решения написаны, в основном, для двумерного случая и включают в себя:

• Конвекция мантии с постоянной и переменной вязкостью (2D): [38] - предложена модель движения литосферных пластин и течения мантии, порожденных контрастами плотности, возникающими из-за охлаждения и утоньшения пластин в зонах субдукции, [51] - выведены аналитические интегральные соотношения между топографией границы и температурой как функцией волнового числа для для слоев жидкости, вязкость которых меняется экспоненциально с глубиной, [22] - проведено сравнение различных эталонных решений для устойчивой конвекции с постоянной вязкостью, конвекции с переменной вязкостью и зависящей от времени конвекции с внутренним разогревом;

• Термохимическая конвекция (2D): [59] - проведено сравнение различных методов исследования термохимической конвекции для течения Буссинеска с бесконечным числом Прандтля, в частности, сравниваются методы отслеживания химической неоднородности, цепочки маркеров и полевой метод;

• Механические и термомеханические течения в каналах и течения Куэтта для постоянной и переменной вязкости (2D): [58] - описаны основные модели и методы для описания указанных течений, [31] - предложен метод робастных характеристик для моделирования многофазных вязко-упруго-пластических термомеханических задач, сочетающий конечно-разностные методы и метод маркеров в ячейках, [29] - описаны основные численные методы расчета указанных геофизических течений;

10

• Течения вокруг деформируемых эллиптических включений (2D): [53] - методом Му-схелишвили получено в замкнутой форме аналитическое решение для изолированного эллиптического включения в сдвиговом течении;

• Нестабильность Рэлея-Тейлора (2D): [50] - изучена неустойчивость Рэлея-Тейлора в геофизических задачах, [41] - исследовано влияние упругости на неустойчивость Рэлея-Тейлора в системе из вязко-упругого слоя, покрывающего вязкий слой меньшей плотности;

• Термомеханические течения в зонах субдукции (2D): [60] - разработаны эталонные решения для сравнения с результатами численного моделирования динамики и термической структуры зон субдукции;

• Спонтанная субдукция на свободной поверхности (2D): [52] - предложены эталонные решения для проверки численных алгоритмов расчета динамики литосферы с суб-дукцией для различных вычислительных кодов (эйлеровых и лагранжевых, конечноэлементных и конечно-разностных, с маркерами и без маркеров;

• Конвекция мантии в евклидовой геометрии (3D): [24] - предложены эталонные решения для сравнения с численными методами исследования конвекции с высоким числом Прандтля в трехмерной декартовой геометрии;

• Конвекция мантии в сферической геометрии (3D): [63] - предложены эталонные решения для проверки численных методов расчета конвекция мантии в сферической геометрии и протестирован конечно-элементный код CitcomS для трехмерной сферической конвекции;

• Бесконечно малые и конечные по амплитуде нестабильности в процессе складкообразования (3D): [42] - показано, что линейная теория корректно описывает только малые амплитуды, для больших необходимы нелинейные уравнения, проведено моделирование литосферных деформаций.

Остановимся на результатах перечисленных работ более подробно. В работе [38] предложена модель движения литосферных пластин и течения мантии, порожденных контрастами плотности, возникающими из-за охлаждения и утоньшения пластин в зонах субдукции. Авторы оценивают количественно движущие силы, связанные с этими контрастами плотности, чтобы определить, могут ли они вызывать наблюдаемое движение пластин. Рассмотрена двумерная модель для оценки влияния реологии и граничных условий. В вязкой модели установившийся режим определяют по обнулению силы трения на границе. Если литосфера имеет ньютоновскую реологию, то силы на вычислительной сетке сильно зависят от ее шага, что приводит к трудностям в интерпретации. Улучшить ситуацию помогает введение твердо-упругой литосферы, разрушающейся при критическом напряжении, в вязкую модель. Модель распространена на трехмерную сферическую геометрию пластин. Скорость твердых

11

пластин согласована с решением вязкого уравнения Стокса на границе литосферы и астеносферы. Величины сил взяты из реальных литосферных измерений. Внутренние контрасты вязкости, связанные, например, с влиянием горячего слоя у нижней границы, считаются величинами меньшего масштаба и не учитываются. Движущие силы, связанные с субдук-цией, локализованы в соответствующих зонах. Силы, связанные с утоньшением пластин, распределены по вязкому слою. Баланс этих сил и определяет движение пластин. Модель представляет собой специфический случай конвекции пограничного слоя.

В статье [51] выведены аналитические интегральные соотношения между топографией границы и температурой как функцией волнового числа для слоев жидкости, вязкость которых меняется экспоненциально с глубиной. Аналогичные соотношения между гравитационными силами и температурой также выведены. Показано, что если вязкость меняется с глубиной не более чем на порядок, хорошо работает приближение постоянной вязкости. При вязкости, меняющейся на несколько порядков, сильно меняются топографические и гравитационные ядра выведенных интегральных уравнений. При очень сильном изменении вязкости можно использовать приближение твердой нижней границы.

В работе [22] проведено сравнение различных эталонных решений для устойчивой конвекции с постоянной вязкостью, конвекции с переменной вязкостью и зависящей от времени конвекции с внутренним разогревом. Авторы сравнивают числа Нуссельта, скорости, температуры, потоки тепла, топографию. Тестируются программы, использующие конечноразностные, конечно-элементные и спектральные вычислительные схемы.

В статье [59] проведено сравнение различных методов исследования термохимической конвекции для течения Буссинеска с бесконечным числом Прандтля, в частности, сравниваются методы отслеживания химической неоднородности, цепочки маркеров и полевой метод. В случае изотермической неустойчивости Рэлея-Тейлора получено хорошее соответствие для начального роста неустойчивого нижнего слоя, однако, развитие во времени и локализация более поздних нестабильностей меньшего масштаба могут различаться в разных подходах. Для моделирования термической конвекции плотного слоя у дна мантии найдено хорошее соответствие для времени нескольких переворачиваний, но потом разница результатов может быть значительной. Эффекты хаотического перемешивания могут привести к отсутствию детального соответствия.

В работе [53] методом Мусхелишвили получено в замкнутой форме аналитическое решение для изолированного эллиптического включения в сдвиговом течении. Решения подходят для несжимаемых полностью упругих или полностью вязких систем. Ограничений на контрасты свойств материалов нет. Решение позволяет найти кинематические (функция тока, скорость) и динамические (давление, напряжение) параметры. Авторы подробно не выводят решение, а сосредотачиваются на его использовании. Авторы получают новые результаты для медленного вязкого течения с включением. Другое важное приложение - эталонные решения для тестирования кодов. Это особенно интересно, так как в модели допустимы высокие контрасты параметров. Все решения представлены также в MATLAB.

В статье [50] изучена неустойчивость Рэлея-Тейлора в геофизических задачах.

12

В работе [41] исследовано влияние упругости на неустойчивость Рэлея-Тейлора в системе из вязко-упругого слоя, покрывающего вязкий слой меньшей плотности. Хотя части литосферы на определенных временных масштабах можно считать упругими, это обычно игнорируется в геофизических моделях динамики мантии. Недавно было показано, что упругие и вязко-упругие свойства могут сильно влиять на литосферные сдвиговые зоны. Задача в статье рассмотрена как численно, так и аналитически. Показано, что упругость при определенном соотношении параметров приводит к ускорению неустойчивости Рэлея-Тейлора. Это происходит из-за уменьшения доли вязкости в деформации вязко-упругого слоя, что наблюдается при временных масштабах, меньших максвелловского времени релаксации = GG где П - вязкость, а G- упругий модуль сдвига. Для тектонических плит Земли, в большинстве случаев данный эффект не наблюдается, однако встречаются ситуации, когда результаты существенно отличаются от чисто вязкого случая. Это подтверждено численными расчетами.

В статье [60] разработаны эталонные решения для сравнения с результатами численного моделирования динамики и термической структуры зон субдукции. Эталонный случай основан на термомеханическом подходе, сочетающем кинематические и динамические методы. Предложены способы изучения влияния граничных условий и реологии на краевые течения и итоговую термическую структуру. Сравнение с численными подходами показывает, что в них необходим точный учет скачков скорости вдоль зон сейсмичности и высокое разрешение для термического пограничного слоя. Небольшая модификация граничных условий для краевых течений также необходима, чтобы избежать сингулярности давления, которая есть в угловых точках у аналитического решения.

В работе [52] предложены эталонные решения для проверки численных алгоритмов расчета динамики литосферы с субдукцией для различных вычислительных кодов (эйлеровых и лагранжевых, конечно-элементных и конечно-разностных с маркерами и без маркеров). В эталонных тестах для моделирования свободной поверхности используется слой с малой вязкостью и нулевой плотностью. Рассмотрена и чисто свободная граница. Показано, что сходимость в расчетах поведения зон субдукции при увеличении разрешения сильно зависит от схемы усреднения вязкости вблизи движущейся реологической границы. Выбор гармонического среднего приводит к наибыстрейшей субдукции, арифметического среднего - к медленной субдукции, а геометрическое среднее дает промежуточный результат. Проверен еще ряд норм средних и выбрана та, которая дает лучшее поведение при увеличении разрешения. Сравнение слабого верхнего слоя и чисто свободной поверхности показывает, что соответствие наблюдается только при высоком разрешении. Сравнение численной модели и лабораторного эксперимента показало, что численный лагранжев подход воспроизводит экспериментальный результат только при достаточно высоком разрешении. Коды, основанные на эйлеровом подходе со слабым верхним слоем воспроизводят экспериментальные результаты если используется гармоническое или геометрическое среднее. Таким образом, использование слабого верхнего слоя для моделирования свободной поверхности дает удовлетворительные результаты только при выборе высокого разрешения и подходящего способа усреднения.

13

В работе [24] предложены эталонные решения для сравнения с численными методами исследования конвекции с высоким числом Прандтля в трехмерной декартовой геометрии. Проведено сравнение численных результатов с результатами лабораторных экспериментов. Рассмотрены случаи бимодальной конвекции с постоянной вязкостью и конвекции в прямоугольной ячейке с вязкостью, зависящей от температуры. Сравниваются числа Нуссельта, скорости, температуры, потоки тепла. Рассмотрено зависящее от времени поведение и предел устойчивости для бимодальной конвекции. Протестированы восемь программ и проведено их сравнение. Также проведено сравнение с результатами лабораторных экспериментов, которое показало удовлетворительное соответствие в пределах погрешности эксперимента.

В статье [63] предложены эталонные решения для проверки численных методов расчета конвекции мантии в сферической геометрии и протестирован конечно-элементный код CitcomS для трехмерной сферической конвекции. Представлены два класса модельных вычислений: стоксовы течения, а также термическая и термохимическая конвекции. Для стоксова течения вычислены скорость течения, топография на границе ядро-мантия при различных степенях сферических гармоник с помощью кода CitcomS и проведено сравнение с решениями, найденными методом матрицы пропагатора. Для термической и термохимической конвекции рассмотрены 24 случая, отвечающие различным параметрам модели: числа Рэлея меняются в пределах от 7 х 103 до 105, контраст вязкости меняется от 1 до 107. Для каждого случая найдены усредненные по времени характеристики, включая числа Нуссельта для границы ядро-мантия, скорости, средние температуры, максимальные и минимальные скорости течения, температура на средней глубине мантии и среднее отклонение. Для термохимической конвекции, в дополнение к термической конвекции, была проведена оценка степени рассеяния первоначально плотных компонент, а также относительное изменение их объемов. Для девяти случаев с малым контрастом вязкости было проведено сравнение кода CitcomS с ранее опубликованными результатами. Для чисел Нуссельта и скоростей результаты отличаются не более, чем на 1 процент. Остальные пятнадцать случаев соответствуют либо большому контрасту вязкости в течение Стокса, либо анализу термохимической конвекции. Ранее подобные результаты получены не были. Эти результаты могут быть использованы для тестирования других кодов. Также проведена оценка эффективности параллельных вычислений для кода CitcomS на суперкомпьютере Ranger (3072 ядра) в Техасском компьютерном центре.

В работе [42] показано, что линейная теория корректно описывает только малые амплитуды (для больших необходимы нелинейные уравнения), проведено моделирование литосферных деформаций. Для больших амплитуд представлены новые нелинейные уравнения. Численное моделирование складкообразования первоначально горизонтального слоя, возмущенного случайным шумом, показало, что оси складок в большинстве случаев перпендикулярны главному направлению укорачивания. Ветвление складок относительно нечувствительно к направлению укорачивания. Трехмерная неустойчивость складок приводит к уменьшению среднего дифференциального напряжения в пределах складчатого слоя. Это говорит о том, что подход "Christmas Tree"для представления напряжения в литосфере может быть непри

14

меним, если при деформации возникают складки.

Различные эталонные решения важны, так как описывают ситуации, с которыми численные алгоритмы могут столкнуться в процессе реального геофизического моделирования. Следовательно, расширение класса эталонных решений (2D и 3D случаи) играет важную роль в создании нового поколения вычислительных алгоритмов, целью которых является выработка общего подхода к деформации твердого тела и течению жидкости, использование адаптивного шага сетки в локальном и глобальном масштабах, лучший учет инерции и т.д. Например, в [45] представлены методы, приемлемые для моделирования больших деформаций геоматериала для некоторой модификации метода конечных элементов. В представленном коде ELIPSIS сочетаются эйлеровы и лагранжевы точки. Неизвестные переменные вычисляются в эйлеровых точках, а лагранжевы точки (частицы) несут в себе информацию об истории переменных. Этот метод идеально подходит для моделирования жидкоподобного поведения твердых тел, которое регулярно используется в геофизике. Представлены эталонные решения из геофизической области. В [26] предложен основанный на MATLAB конечноэлементный метод MILAMIN для расчета двумерных геофизических течений. MILAMIN позволяет достичь хорошей точности в задачах с неоднородным геоматериалом. Проведена оптимизация глобальных матричных вычислений. Эта техника оптимизации применима ко многим ситуациям, в частности, в работе продемонстрировано ее применение к термическим и несжимаемым течениям (Стокса). Проведено сравнение с другими открытыми кодами в данной области. В [55] предложена модель конвекции сжимаемой мантии с большим контрастом вязкости в трехмерной сферической оболочке на базе сеточного метода. Используется yin-yang код. Решение для скорости и давления строится мультисеточным методом, а температурное поле - методом конечных объемов. Сходимость мультисеточного солвера в условиях больших изменений вязкости составляет проблему. Предложена новая интерполяция для давления, которая серьезно улучшает ситуацию. Тесты показали, что он работает при величине изменения вязкости до 19 порядков. Соответствие с результатами уже имеющихся кодов проверено. Показана возможность распараллеливания для работы на кластерах. В [54] предложена модель высокой точности для описания течений мантии разных масштабов и тектоники плит с использованием параллельных вычислений. В глобальное описание течения включено движение плит. Холодные термические аномалии нижней мантии соединены с океаническими плитами узкими каналами, что меняет скорость океанических плит. Потери за счет вязкости в пределах изогнутой литосферы составляют 5-20 % от общих потерь в литосфере и мантии. В статье [30] предложена новая конечно-разностная схема с адаптивной сеткой для расчета несжимаемого стоксова течения. Новая схема дискретизации позволяет минимизировать вычислительные затраты. Численным способом продемонстрировано, что адаптивная сетка не приводит к осцилляциям давления и обладает тем же порядком точности, что и соответствующая неадаптивная схема дискретизации. На ряде примеров показано преимущество адаптивной сетки перед обычной. Приведены геофизические примеры как литосферного, так и планетарного масштабов.

Цель настоящей работы - существенное расширение класса эталонных решений для те

15

чений Стокса в двумерном и трехмерном случаях. В отличие от предыдущих исследований, мы предпочитаем не брать за основу какую-то конкретную модель, а получить аналитические решения задачи в общем случае, что поможет создать значительное число тестовых задач. Планируется решить задачу Стокса в случае несжимаемой жидкости для (а) линейно и (б) экспоненциально меняющейся вязкости. Показано как создать двумерные и трехмерные тестовые задачи на основе общих аналитических решений. В конце работы приводятся результаты проверки некоторых широко используемых численных алгоритмов (см. [31], [32]) на основе полученных эталонных решений.

Введение диссертации (часть автореферата) на тему «Решения, сосредоточенные в окрестности пространственно-временных лучей, и эталонные решения в задачах о движении жидкости»

Структура работы.

Работа изложена на 87 стр., состоит из 6 глав, 33 рисунков, 18 таблиц, списка литературы, 4 приложений (19 стр.).

16

2 Квазифотоны для волн на поверхности тяжелой жид-

кости

В данной главе строятся квазифотоны (особого класса волновые пакеты) для волн на поверхности тяжелой жидкости. Эта задача актуальна в связи с исследованием распространения и развития волн в океане. В частности, квазифотонное решение можно рассматривать как начальную линейную стадию развития волны цунами. Квазифотоны известны для волн многих типов (см. [1], [3], [4], [9], [14], [16]). Строить формулы для квазифотонов можно было бы на основе методики дифракционного пограничного слоя или техники комплексных лагранжевых многообразий В.П. Маслова [16].

В настоящей работе используется иной подход (см. [1]): квазифотон представляется асимптотическим пространственно-временным лучевым (ПВ-лучевым) рядом, члены которого являются формальными степенными рядами. Целью данной работы является построение указанного асимптотического ПВ-лучевого ряда.

Для достижения поставленной цели решались следующие задачи:

1. Физически осмысленный и математически корректный вывод энергетического соотношения для данной системы.

2. Решение уравнений переноса в рамках вещественного варианта пространственно-временного лучевого метода (ПВЛМ).

3. Построение комплексного варианта ПВЛМ.

4. Построение квазифотонного решения.

2.1 Энергетическое соотношение

2.1.1 Постановка задачи

Для полноты изложения напомним некоторые сведения из теории волн на поверхности тяжелой жидкости. Эти сведения известные или близкие к известным. Но нам не удалось найти доказательства этих фактов в литературе. Поэтому выведем их сами.

Будем исходить из известных уравнений для идеальной жидкости, заполняющей область:

—H(x1,x2) < z < 0, < x1,x2 <

= —grad(p — po) — gradpgz, (2.1)

где р = const - плотность жидкости, g - ускорение свободного падения, p0- давление на поверхности жидкости.

Течение предполагается безвихревым:

"V = —ргп^Ф, Ф = Ф(х1,х2,z,t). (2.2)

17

Жидкость предполагается несжимаемой: divdX = 0, поэтому потенциал скоростей Ф будет гармонической функцией:

ДФ =

д2 д2 д2 \

dx2 dx2

Условие непроницаемости дна (при z = — H(x1,x2), H > 0):

vn]z=-H = 0, где n - нормаль к поверхности дна.

Нормаль к поверхности дна: N = , ддН, l) ,

Vn = — (gradФ, ) = — g- cos (n,Xij — ддХФ- cos ан dx2

1)

дн

дФ__

дх' /(dH )2+( Й )2+1

vn]z^H = 0, что эквивалентно

J

(ф Xi) 7(W+()2+i

n = ]N ].

^n, — ддф cos ^n, V)

_ дФ i

dz

n, Х2

/дФ dH дФ \ dz dx1 dx1

dH дФ\

+ дХ2дхх)

= 0.

z=-H

(2.3)

dx2

7(^W)T

Из уравнения (2.1) следует, что:

рдФ = (р — Ро) + pgz + x(t), где x(t) - произвольная функция, не зависящая от x1,x2,z. Положим x(t) = 0.

Тогда руф = (р — ро) + pgz.

В частности, при р = р0 на поверхности жидкости:

дФ

Ж

=

z=h

(2.4)

где h- значение z при р = р0, т.е. возвышение волны над поверхностью z = 0

(уровень z = 0 отвечает поверхности жидкости без волн). Считаем, что h-мало, а функция Ф достаточно гладкая (т.е. мало меняется с изменением z), поэтому равенство (2.4) можно записать в виде: уф ]z=0 = gh, что эквивалентно формуле:

д2Ф _ dh at2 ,=о = .

(2.5)

На поверхности вертикальная компонента скорости vz равна

дһ dt ,

то есть:

д^ _ _ дФ

д? = Vz = — д^.

(2.6)

Из (2.5) и (2.6) получаем краевое условие для Ф при z = 0:

д2Ф дФ

dt2 + дд7у)

= 0.

z=0

(2.7)

Итак, решаем уравнение Лапласа при следующих краевых условиях:

ДФ = 0

д2Ф + g дФ at2 + g az

=0

z=0

дФ + H дФ + H дФ дz + дх1 + Hx2 дх2

=0

z=-H

18

2.1.2 Энергетическое соотношение

Рассмотрим выражение для энергии жидкости, находящейся в бесконечно тонком столбике xi < Xi < xi + dxi, x2 < < x2 + dx2, —H < z < h

Эта энергия складывается из потенциальной энергии Ep и кинетической энергии Ek. E = Ek + Ep = Edxidx2, где E - плотность энергии, рассчитанная на единицу площади

поверхности.

Найдем выражение для потенциальной энергии. Считаем, что глубина H медленно меняется. Выберем "нулевой уровень” так, что потенциальная энергия столба жидкости равна 0,

если его верхняя граница находится на уровне z = 0.

h h2

dEp = g = р—gdxidx2

Согласно формуле (2.4): h = фу

2

dxidx2

dEp =р g( д)

Найдем кинетическую энергию . Поскольку "V = —gradФ, для элементарного объема

кинетическая энергия такова: р(дга^ф) dzdxidx2. Тогда для столбика сечением dxidx2:

/ h

dE^ = / у

(

Т.к. h мало, то dEk J р(дга^ф) d^ dxidx2. Полная энергия данного столбика: dE =/ р $2

Һ g

) ......

7 (grad$)2, \ ,

/ р--------dz dxidx2.

-H 2 /

Полная энергия цилиндрического столба жидкости такова:

2 dxidx2 +

z=0

(2.8)

E g/2)Ф

z=0

0

Плотность энергии: E^ J Р(дго^Ф)^. +

-н

Нам потребуется выражение для фЕ.

р (дга^Ф)2^х1^х2^л

Q

(2.9)

-Р Ф2

2g Ф

z=0

dE dt

2 t

рФ.Ф

z=0 g

tt

z=0

z=0

dxidx2 р(grad$)2dxidx2dz

Q

dxidx2 р ($xi$txi + $X2$tX2

Q

+ $tz) dxidx2dz

о

В первом слагаемом воспользуемся краевым условием (2.7), а во втором - проинтегрируем

по частям. Получим:

уд = — JJ рФtФz ]z=0dxidx2+

z=0

JJ рФt(Фxl cos(d^, Xi) + cos(d^, X^) + cos(d^, ))dE—

dQ

рФt (ФХ1Х1 + ФХ2X2 + ) dxidx2dz

Q

19

В силу гармоничности Ф (ЛФ = ФЖ1Ж1 + ФХ2Х2 + Ф^ = 0) последнее слагаемое в вышеприведенной формуле равно 0.

Учитывая, что

Фж1 cos(xt, Xi) + Фж2 cos(dt, X2) + Ф^ cos(dt, , dn

получаем: dE Г Г , /У дФ дҒ = РФ*Ф^tz=odxidx2 РФУdndE z=0 dQ

Второй интеграл разобьем на сумму трех: по поверхности жидкости, по боковой поверхности цилиндра и по дну.

= * JJ рФуФ ]z=o + JJ рФ^Фdxidx2+

+ J J рФ; (Фх1 cos(6^, Xi) + Фх2 cos(dt, Ж2)) dE + J J рФ^дПdE

lateral surface z=-H

Первые два слагаемых сокращаются, последнее слагаемое равно 0 в силу непроницаемости дна: dn L-н = 0.

J рФу (Фх1 cos(6t, X1 ) + Фх2 cos(xt, x2)) dE =

lateral surface

0

= J d/ J (рФ^Фж! cos(dt, Xi)+ рФ^Ф^2 cos(dt, X2))dz,

dD -H

(2.10)

где dD - граница основания цилиндрического столбика жидкости. Введем новый вектор t - вектор Умова:

S =(S1,S2);

дФ дФ dtdxj

dz.

(2.11)

Тогда формула (2.10) записывается в виде:

=—/ (*^' )d/.

dD

(2.12)

Применяя теорему Остроградского-Гаусса, получаем:

J (t, dt)d/ = JJ div^Sdxidx2 =

dD D

S*D — площадь D, M- некоторая точка внутри D.

Пусть область D сжимается к некоторой точке N. Тогда:

- SD (по теореме о среднем) M CD

= lim / (^t, X^)d/.

N(xi,x2) D^N SD У

dD

Плотность энергии с точке равна: 1 E E = li^ ——. D^N SD

— + dzvS = 0. dt

С учетом последних преобразований формулу (2.12) можно переписать в том виде, в котором

обычно записываются различные “законы сохранения” математической физики:

(2.13)

20

2.1.3 Вектор Умова S как “вектор работы”

Равенство (2.12) имеет ту же форму, что и знаменитая теорема Пойнтинга в электродинамике (в отсутствие тока). Приблизительно в то же время, что и Пойнтинг, формулу вида (2.13) вывел в случае уравнений теории упругости Н.А. Умов (1874)(см. [12])

Вектор у несколько условно интерпретируют как вектор потока энергии. Умову же принадлежит и другая интерпретация S*. В случае теории упругости он показал, что “поток” вектора J Sndl (n - нормаль к боковой поверхности цилиндра) можно интерпретировать dD

как работу в единицу времени упругих сил, приложенных к S со стороны, противоположной той, куда направлена нормаль n. В нашем случае вектор "S можно трактовать аналогичным образом. Проделаем соответствующие выкладки.

— / (^, )dl = / (-Sn)dl

dD dD

0

dD -H

\ cos(n,xi) +

dz.

d Ф дФ . . дФ . .

л— Л— cos(n, -— cos(n, Ж2).

Подставляя формулу pффф = (P — p0) + pgz в интеграл и учитывая, что = —дг^Ф, = —дго^Ф - , то есть = — дф,

получаем:

j ғ )di = p дФдф ds = jj (p—po) дф ds+Д pgz дф ds =

dD S S S

= jj (p — po)(—Vn)ds + jj pgzдфds = Wi + W2.

S S

Здесь S- боковая поверхность. S = {x.y.z : — H(x1.x2) < z < 0.x1,x2 G dD} dS = dzdl — элемент поверхности S.

W1 интерпретируется как работа сил давления над жидкостью в объеме Q в единицу времени. Интеграл W2 преобразуем:

ГГдФ ГГ д Ф /У .

W2 / pg^—dS = / / pg^—dS = p^ / / z(gradФ) - WdS.

S до до

(2.15)

dQ— вся поверхность цилиндра. Интеграл по верхнему основанию цилиндра равен 0 в силу того, что z = 0, интеграл по нижнему основанию цилиндра равен 0 в силу условия непро-

ницаемости дна (2.3). Применяя к (2.15) формулу Остроградского-Гаусса jj "F - T^dS =

jjj div^dx1dx2dz, получаем: до

о

W2 = pg^^/* div (z(gradФ))dж1dж2dz.

о д д д

div(z(gradФ)) = — (zФ^1) ^Фхз) (zФz) =

dxi dx2 dz

= z У + z Йф + z ЙФ + Ф, = ^АФ + Ф,

(АФ = 0 в силу гармоничности Ф).

21

Тогда W2 = pg dx1dx2dz.

Q

dD S

(p - Po)(-Vn)dE + pg

///

(-vz )dx1dx2dz.

(2.16)

Последний интеграл интерпретируется как работа силы тяжести над жидкостью в объеме Q в единицу времени. Итак, увеличение в единицу времени энергии в объеме Q складывается из работы в единицу времени сил давления (“внутренних” сил жидкости) и сил тяжести (т.е. силы, внешней для жидкости). Построения здесь согласуются с общими положениями о работе внутренних и внешних сил в деформируемых средах (см. [9]).

2.2 Пространственно-временной лучевой метод для волн на поверх-

ности тяжелой жидкости

Пусть глубина H - медленно меняющаяся функция: H H(бх1,бх2), где 6 - малый

параметр задачи. Анзатц ПВЛ метода - это разложение

Ф - Л (т,ф ,^2,z)(i6)j .

j=0

Т — 6t, Ф — 6X1, ^2 — 6X2.

Роль частоты и волнового вектора здесь играют

дӨ д 7 Ө ) д 7 дӨ дӨ ) ^7 Ө \ 7 д д )

= — дТ = — , <%) = Д б)' \ дъ 'аД)'

(2.17)

(2.18)

Подставляя анзатц (2.17) в уравнение Лапласа АФ = 0 и приравнивая коэффициенты при

одинаковых степенях, приходим к последовательности равенств:

ТоФо = 0, ТоФ1 + Т1Фо = 0,

' LoФ, + Т1Ф;-1 + L2 Ф,-2 = 0'

АоФ = - к2Ф,

Р1Ф = 2V'ӨV'Ф + А'Ө - Ф,

д2

д72,

Краевые условия (2.3), (2.7) приводят к соотношениям:

V =Г -д-А) А' = + -

V д^1 'д7^ ' д72 + д72

Р2Ф = -А^Ф,

- 76)2 7Ж)2

7 дФ, . д2Ө^ дӨдФ,

- ^Ф'+ дТ2 Ф'-1 + 2

j-i

д2 Ф,-2^ д2т /

= 0,

^=о

дФ

+ V'H^Ө • Ф^-1 - VH^Ф^-21

= 0,

z=-H

j = 0,1,2,..., Ф-1 = Ф-2 = 0.

(2.19)

(2.20)

(2.21)

(2.22)

(2.23)

(2.24)

д2Фо

д-2

- к2Фо = 0,

дФо

^2Фо^)

= 0,

^=о

д Фо д-

При j = 0 мы получим необычную спектральную задачу, где роль собственного числа играет

^2 и входит это собственное число в краевое условие:

= 0.

z=-H

(2.25)

22

Из (2.25) легко следует, что

Фо = Ao (т, , ^2) ' ch [k (z + H)] , z2 = gk - th (kH)

(2.26)

Таким образом, дисперсионное соотношение имеет вид:

z = z (k) = ^g^th (kH), где = (ki,k2) = , .

dz dz \ dz h

dk1, dk^ dk k

Выведем соотношение ( S) = —Г (E). Напомним выражения для основных величин.

(2.27)

о

Е= Р / (gradФ)2 dz + Ф2

</ g z=0

-H

Ф^Фж^. dz.

дӨ

дт

дӨ d^i

1 дӨ

- дt

1 дӨ

- дж^ дz дk^.

-Өо.

Ө^.

z = —

В нулевом приближении имеем:

Ф = Фо^1 ө

Ө

— Фо cos -

-

(2.28)

Поясним, что означает знак “—” в (2.28). С комплексными экспонентами работать удобнее, чем с тригонометрией. Соглашение: работаем с Ф = ФоЛө, а физический смысл приписываем ЯеФ.

Тогда

Ө Ө 1 дӨ

(^Ф)^ - ^)j cos - - Фо sin - - (2.29)

z - у- д^ j у

дФо Ө дӨ

Ф/ - cos---Фо Si^ - —

дt - - дт

Функция Лагранжа L есть разность кинетической и потенциальной энергий, то есть в нашем

(2.30)

случае имеет вид:

о

L = Р / ОгаЛФ)2 dz - Ф2

2g т=о

-H

Подставляем в (2.31) выражения (2.29) и (2.30):

(2.31)

_р

2д

L

(дфф0 cos 2 - Фо sin 2

\ 2 о

Фо sin (2 - dz + Р J ФZdz-

-н

. дә^2

дт/

(2.32)

т=о

23

где ЭФ о _ дФр дт дФр _ дФр dt d-rdt'dxj dgj '

Преобразуем формулу (2.32):

L - Р у Е ((ЭФ?)

- 2g ((ЭФ° )2 cos'

+Р J ф^"^-

-H

,2 Ө ' е

2 cos2 Ө - 2Фоcos Ө sin Ө + ф0 (sin2 z=0 - 2тФ?фоЭӨ cos е sin е 1^=о + ф2 (g)2 sin2 Ө

elz=0

е) "Ъ-) +

z=0

cos е и sin е - быстро меняющиеся функции (ибо 6 - мало и малое изменение Ө вызовет значительные изменения е). Остальные функции в формуле (2.32) меняются мало. Считая, что эти функции на периоде изменения cos Ө, sin Ө постоянны, усредняем лагранжиан L по

данному периоду:

Принимая во внимание, что:

2

cos - d

6

2п

6

_ П

ө ' ө, co^ si^ - d

с с

0

_0,

получаем:

+1 Фо(dj) ^dz + Р

(L) - Р y E C ()

-Hj=i\

- (1 (ЭФР0 )2 + 2 Ф2 (ЭӨ )2)

/ фр,Л -

z=0

Слагаемыми дФ° _ - 6 и дФо _ - 6 можно пренебречь в силу того, что 6 << 1.

dt дт dxj dxj '

Итак, среднее значение лагранжиана представляется в виде:

dz+Р

то есть

0 / 0 \

(L) - 2 у 2 (ө2 + ө2^ dz + Р / у Ф2 dz\

-H '-ДГ /

0

Ф2 dz

-н

- 4) Ф2ө0

(2.33)

z=0

Второе слагаемое в (2.33) интегрируем по частям:

/ Ф2dz _ ФФ^]Z=-H - У ФФ^^dz-

-н -н

Учитывая анзатц ПВЛ метода,

Ф — e Ф^' (i6)j,

j=0

для Ф^ получаем:

Ф

z

УӨ _ e

d Фр + dz +

<9Ф1. ,

i6 + dz

дФ2

dz

(i6)2 +

УӨ e

дФ0

dz

24

Здесь мы пренебрегли слагаемыми более высокого порядка малости. В силу (2.25) имеем:

д Фр dz

= 0,

,=-н

откуда следует, что:

Ф, tz=-H = 0-

Преобразуем интеграл:

/ Ф,dz = ФФ, lZ=-H

-н

= -1ФФ^ g tt

z=0

0 0

- J ФФ,,dz = ФФ,],=0 - / -н -н

0

-J ФФ,, dz-

-н

ФФ,, dz =

Последнее равенство выполнено в силу граничных условий (2.7):

Ф,],=0 = - -Ф^

g

z=0

Найдем Фtt. В силу (2.28) имеем: Ф — Ф0 cos Ө. Следовательно,

Фй '

1 д Ф0

6 д?

дФэ д?

Ө

cos----

6

Ө sin - -

6

дӨ dt

- —-Ф0 co^ -

62 6

1^ - Ө

-Ф0 sin -

66

дӨ

д? '

. Ө

----Ф0 sin -

6

дӨ д?-

/д^2 1^ .

Ы - 6Ф0 sin 6' = -Ф0 cos 6 . (дӨ)2 = -Ф0 cos 6 . Ө2.

Ө —Ф0 cos -62 6

Ө /д^2

W

6

Ө /д^2

W

Здесь мы пренебрегли величинами меньшего порядка.

Кроме того,

- ^(z_ (ф0 co, 6^

д2Ф0 Ө Ө

-7——- co^ = k Ф0 co^ -дz2 6 6

0

С учетом вышесказанного, получаем следующее выражение для J ФZdz:

-н

0

Ф2dz = -Ф0Ө0 cos2 -g

к2Ф^ cos2 -dz.

-н

6

Тогда

0

ФЖ + d2)dz + Ф2ө0

4g

-н

z=0

0

4 у k4gdz - 44 Ф^^^

-н

= 0-

z=0

(2.34)

Итак, показали, что (L) = плотности полной энергии:

0, то есть: (E^) = (Ep). Аналогично находится выражение для

(E)- Р

+ Ф2Ө0

J Ф0(Ө2 + Ө2№ + Ф2Ө0

-н g

= Ф2 2g Ф°

z=0

- 4 J к2Ф2dz+

-H

(2.35)

0 =

z=0

z=0

z=0

25

Нам потребуется соотношение для компонент вектора Умова S (см. (2.11)):

(S.) = -P 1 (dt ЛеФdX- Леф) dz -

^-P / Фо sin Ө . дӨ • Ф. sin Ө . (-)d^ .

то есть

/ ° % \ ° ) ^-^ (-^)sin2 -Ф0 - d^ -

-H -H

(2.36)

о

Теперь выведем связь между J Ф(Дл и Ф0^=о. Функция Фо удовлетворяет условиям:

-H

Дифференцируем (2.37) по к:

(4Ф - к2Фо = 0.

(у дФ0 - ^2фо)1,=о=о. дФоI _ о

dz lz=-H '

д2 дФ^_ ,2 дФо

дл2 дк дк

2кФо _ 0,

д дФо 2 дФо

ддл дк дк

_ 0.

z=0

д дФо дл дк

_ 0.

z=-H

Домножим (2.38) на Фо и проинтегрируем по л от -H до 0:

о о о

/Фо' дд dz - к2 / 'Ф"^2 _2к /Ф2^".

-H -H -H

(2.37)

(2.38)

(2.39)

(2.40)

(2.41)

2^Ж Ф")

Проинтегрируем по частям в первом слагаемом (2.41):

о

Г Ф _дГ дФо dz _ Ф д дФо 1о

J Фо дz2 дк dZ Фо дz дк ]-H -H

_ Фо

одz дк I-H

1Ф д дФо дФо дФо *\

\Фо дz дк дк дz 7

о

Г д дФо дФо dz ___

J 3z дк ' 3z dZ

-H

о

дФо дФо 1о + Г дФо д2Фо dz _

дк д^ —H d дк дz2

-H

о

+ 1

-H

z=0

dz

дФ° - k2Фоdz _

(Ф4 g

дФо

дк

дФо дк

Фо

z=0

+ к2

1 дфк0 Ф,dz _

-H

2^ 2

g дк Фо

z=0

+ к2 1 дфО Ф"dz

-H

I 2^ д^ /R \ + Ддк ФД<

^4

g

Подставляем в (2.41):

2^ д^ g дк

Ф2о

о

+к2 /

-H

дФо дк

ФсДл

дФо дк

о

- ФсДл _ 2к

-H

Ф2dz,

z^

26

то есть

2^д^ g dk

Ф2

z=0

о

(2.42)

Вернемся к выражению (2.36) для среднего значения Sj:

(S) = р^k^ Ф2^л = рфkj - gdkф2 -H g

= № - (E) = (E) = Vgrj (E).

z=0

= P^2 ф2

dk k % 0

Мы пришли к важному равенству:

(Ф = - (E),

(2.43)

где (E) = (Ek) + (Ep).

2.3 Уравнения переноса

2.3.1 О решении уравнений переноса

Решение уравнений переноса мы приведем, следуя монографии [3]. Рассмотрения здесь достаточно стандартны. Обратимся к уравнениям и краевым условиям (2.25) - (2.27) при

д2ф

—1 - к2Ф1 = -2VӨ^Фо - Д'Ө - Фо, dz2

д Ф1

dz

- Ф1)

g /

= -g 1

z=0

2дӨдФ0

2 дт дт

д2Ө^ \

+ аГ2 Ф0)

z=0

д Ф1

dz

-VH^Ө. Ф01,=-н

(2.44)

(2.45)

(2.46)

z=-H

j = Е

Задача (2.44) - (2.46) нахождения Ф1 - это неоднородная задача Штурма-Лиувилля, причем соответствующая однородная задача имеет нетривиальное решение Ф0 = А0 (т, ^1,^2) * ch [k (z + H)]. Соответствующее условие разрешимости и приводит к уравнению переноса для Ф0. Умножим (2.44) на Ф0 и проинтегрируем полученное выражение от -H до 0:

0

-н

0 k2Ф^ dz = -

-н

Ф0 (2^Ө^Ф0 + Д'ө • Фо)dz

(2.47)

Интегрируя по частям левую часть равенства, получаем:

0

J ф0

-н

(дУ -k2 Ф^ dz=ф0

0 dz

iz=0

)z=-H

/ dy дф0dz k2Ф,Фоdz =

-н -н

(ф0дУ - ф1 дФ°)

iz=0

)z=-H

0

+ J

-н

Ф1

д2Фо dz2

- k2Ф^ dz = (Ф0дф!

дФоА

1 6^7

iz=0

)z=-H

(2.48)

Ф

Таким образом:

д Ф1

дz

дФ0

- Ф-щ)

z=0

z=-H

0

Ф0 (2VӨ^Ф0 + Д'ө • Фо)dz

-н

(2.49)

27

Принимая во внимание формулы (2.25), (2.45) и (2.46), левую часть равенства можно запи-

сать в виде:

z=0

(Ф0ЗФ!- Ф. ЗФс)]z=0 =Ф^ДФ. -

\ 0 dz 1 dz/lz=-H 0у g 1

— Ф1 дФ;0l,=° — Ф'0 (—V'H V'^. Ф0)1,=-н = Ф^ Ф. — g (2 дӨдФ?+й У))

1 (2 дӨдФ0 + й ФР)

-Ф. ЗФс] =

1 dzlz=-н

- Ф. Ф0 + V'HV# • Ф01,=-н.

z=0

(2.50)

z=0

Итак, (2.49) переписывается в виде:

-(2дӨдФ0 + фе)

дтдт

0

= Ф0 (2V'^VФ0 + А'ө • Ф0)^л.

-H

+ Ф^'И ^Й1,__н =

z=0

(2.51)

/ 0 \ 0 0

V'^ Ф0^л - V'Ф0^л - А'ө + ^ө • Ф0^л =

у -н у -н -н

= J Ф0А'ө^л + ^ө / 2Ф0^Ф0 — Ф0]z=-н - V'(—И)^

0

У Ф0 (Ф0А'ө + 2^ө^Ф0) dz + V'H • ^ө • Ф0]z=-H.

-н

Слагаемое в (2.51) при z = 0 - это производная

(2.52)

(-д^) Ф0)

z=0

Ф0 / дө g у дт

дФ0 д2ө

"дҒ + дТ2 °)

z=0

что эквивалентно уравнению:

дТ(д (-э Ф')

z=0

+ ^V'^ Ф0^^ - ^'ө^ = 0.

(2.53)

Домножим обе части равенства на р:

рЕД Ф0

2дДу 0

z=0

+ ^v', / Ф2^л • ^^ө • 1 р = о.

(2.54)

Учитывая, что: (Е) = Фд

z=0

- р и J Ф(Дл - S - р можно переписать (2.54) в виде:

\ / -н

-1 (^) С'Ч) о

д^^ у у у

(2.55)

Воспользуемся тем, что:

(S) = —Г (Е) и = — дӨ = —ө0 и запишем (2.55) в виде:

Д (Е) И(Е) у,Д =° дт + Д , J ,

Ут

Ө0

что эквивалентно:

у _д_ Ө()\дт

(Е)+ (v', (Е) vgr + (Е^ дт^ + vgr

V^)=°.

28

(dr + (vgr.= 0, где s - координата вдоль ПВ-луча (расстояние от заданной точки на ПВ-луче до начальной поверхности).

Итак,

д

— (E) + (V, (S)) = 0. (2.56)

Уравнение (2.11) можно записать в виде классического уравнения Гамильтона-Якоби:

+ H&к,) = o, н = ө, = ө, = g-.

(2.57)

Канонические уравнения, с помощью которых решается уравнение (2.49), в частности

дают:

Решение (s) канонической

дН ddr =0 ddi дН dd ^н ӨдН

дө. d7 = . d7 = - . dS = - + дө!.

системы (2.58) задает пространственно временной луч.

-dsөт = dsн(^1,^2,ө1,ө2)= E (g+ dH

i=1,2

7дИ dH + ди/ =0

i=1"E + 0

Равенство dS = 0 показывает, что уравнение (2.55) можно записать в виде:

^SE^ + д^(;'=дН) -

Уравнение (2.59) - это есть равенство нулю дивергенции вектора A с компонентами:

(E). (E) ^. j = 1.2.

(2.58)

(2.59)

(2.60)

Введем в пространстве т. ^,.^2 систему координат a,. a2.s, связанную с решениями канонической системы.

Мы решаем каноническую систему (2.58) , взяв в качестве начальных данных

т )s=0 0 7ts=0 ai. d)s=0 d0(a1.a2).

Ө1 ts=0 =

дӨо да2' .

=

т 7 ds

Значения s.a,.^ мы и возьмем в качестве координатной системы. Вблизи плоскости z = 0

эта система регулярна. Вектор А, имеющий в этой системе координат компоненты (z) . 0. 0

в декартовой системе т. . ^2 будет иметь компоненты (2.60). Напомним, что дивергенция вектора в декартовой системе координат ж1. ж2. . xm, имеющая классическое выражение div A = Y2 ддХд, имеет выражение

i=1

divA = 7 (A"d).

(2.61)

. жт', а J- якобиан: J = '.,x

1д

J iiS((E) J ) = 0.

J = D^Xi.Xg)

D(s. a,. a2)

где A!' - компоненты вектора А в системе координат ж1' x , а J -якобиан: j = ^(^1'-.

Запишем равенство (2.59) в системе координат, принимая во внимание равенство (2.61):

(2.62)

29

Поэтому на пространственно временном луче т = s, Xj = Xj(s, a1, a2), j = 1,2:

(E) J = (E) Jl„0 .

(2.63)

В силу того, что (E) = 2 (E) = P

2^2 ф2

g 0

и пользуясь формулой Ф0 = A0 (т, ^1, ^2)-ch [k (z + H)] z=0

(см. (2.26)) (2.63) можно записать в виде:

J - A0ch2(kH) = j. A2ch2(kH))s=0. (2.64)

Обозначим Д0(а1,а2) = ^/( J - A0ch2(kH))s=0. Получаем, что: A0 = ^/jCy'kH), где Ф0 = Ach(k(z + H)).

2.3.2 Расчет высших приближений

Высшие приближения ПВ-лучевого метода находятся так же, как и в монографии [3]. Опишем кратко соответствующие построения.

Выпишем уравнения для определения Ф^'+1:

д 2ф

j2+1 - k2Фj+1 + 2VW% + Д'Ө • Ф^ - Д%_1 = 0,

(2.65)

^дӨдФ, д2Ө^ д2 Ф

Р j+1 + 2+ дт2 j - *7Й2

j-1

= 0,

z=0

д Ф'+1 дz

+ V'HV'Ө - Ф^' - V'HV'Фj-1

= 0.

z=-H

(2.66)

(2.67)

Пусть Фj-1 найдено на предыдущем этапе рассмотрения, а Ф^- определено с точностью до

решения однородной задачи:

Ф2' (т,^1 ,^2,z) = Ф2'(т,^1,^2,z) + Aj (т,^1,^2)Ф0(т,^1,^2,z).

(2.68)

Задача Штурма-Лиувилля (2.65)- (2.67) для нахождения Фj+1 такова, что соответствующая однородная задача имеет нетривиальное решение Ф0. Такая задача имеет решение тогда и только тогда, когда выполнено условие разрешимости, которое выводится на том пути, на котором выводилось условие (2.53). Соответствующие вычисления приводят (см. [3], глава 3) к соотношениям вида:

Р (E) Д + ^j(s) = 0.

где ^j (s) определено на предыдущих этапах рассмотрения рекуррентной системы уравнений (2.65)- (2.67).

2.4 Комплексный вариант лучевого метода и квазифотоны

2.4.1 Общая схема построения

В этой главе мы непосредственно переходим к построению квазифотонов. До этого мо-

мента для удобства записи индексы у координат ставились снизу. Здесь же будем ставить их

наверху.

30

Пусть в ПВ-лучевом разложении ТтӨ > 0, причем знак равенства имеет место лишь на некотором ПВ-луче l. Тогда при 6 1 частные суммы ПВ-лучевого ряда

м

Фм - еФ'.

j=0

(2.69)

были бы существенно отличны от нуля лишь в малой окрестности l, т.к.

ie е

е-1

“Волна” Фм представляла бы собой волновой пакет - квазифотон, движущийся вдоль поверхности жидкости с групповой скоростью v = ((dr, ddr). Строить комплексные решения уравнения эйконала сложно в силу того, что теория характеристик - вещественная теория.

Однако, некоторый аналог техники ПВ-лучевых разложений здесь применить возможно,

если Ө и Ф' считать не функциями, а формальными степенными рядами. Этот подход уже применялся в работах ([3] - глава 3, [1] и [7]). Отличие похода, применяемого в настоящей статье от подхода в [1], [7] - отказ от использования так называемых квазисопряженных решений, что позволяет несколько упростить процедуру построения квазифотонов. Квазифотон у нас - это разложение вида (2.17):

Ф е

Ф'(W1,^2,zW, j=0

где в свою очередь Ө и Ф' - формальные степенные ряды (см. [10]):

Ө = ө(0) + ө(1) +..............

Ф' = Ф(0) + Ф(1) +

(2.70)

(2.71)

(2.72)

Здесь

Ө(0) = Ө0(т), Ө(1) = Ө'(т)п' + Ө2(т)п2, Ө(2) = Ө11(т)(п1)2 + 2Ө12(т)п'п2 + Ө22(т)(п2)2, ........

или в общем виде:

Ө(г) = Ө'1 (т )п'' - - , Ф^ = Ф'1 (nW1 - - , (2.73)

п' = - Пт), i = 1,2. (2.74)

Предполагается, что '(т), i = 1, 2 - это некоторый фиксированный ПВ-луч. Именно в

окрестности этого ПВ-луча мы построим квазифотон. Аналог переменных ф впервые был введен А.П. Качаловым [14] в другой задаче. Мы будем предполагать, что ПВ-луч - “опорный ПВ-луч квазифотона” проходит через начало координат т = 0, ' = 0, i =1, 2. Введем новые координаты s, щ, a2. Здесь щ, а2 - декартовы координаты пересечения ПВ-луча, проходящего через заданную точку, с плоскостью т = 0. s = т .В координатах s, щ, а2 точке т = 0, ' = 0 соответствует точка s = 0, щ = а2 = 0. При замене переменных т, ф,П2 на s,ai,a2 Ө и Ф' превращаются в ряды по степеням щ,а2. Подставляя разложение (2.70) в уравнение и краевые условия, мы придем к тем же уравнениям (2.19)-(2.24).

31

Прежде всего займемся уравнением для Ө = Ө(0) + Ө(1) + , которому можно придать вид

уравнения Гамильтона-Якоби:

дӨ + H (< '.<2.Ө,.Ө,) =0, Н= k <2.75)

Наши построения здесь практически совпадают с построением в монографии [3] (глава 3).

Для того, чтобы разложение (2.70) представляло квазифотон, мы предположим, что Ө(0) и Ө(1) вещественны, а квадратичная форма 1тӨ(,) положительно определенная (т.к. квази-фотонное решение должно экспоненциально убывать при удалении от ПВ-луча, а если бы Ө(1) было комплексным, то е'^^ в каком-то из направлений обязательно возрастало бы). Экспоненциальное убывание решения по всем направлениям обеспечивает Ө(,).

В окрестности ПВ-луча перейдем к переменным т, п' = <' — <' (т):

дӨ дТ

dd d^' дп' dr

(s' + < '(т).

дп'у

= 0.

(2.76)

+ Н

2.4.2 Нулевой и первый порядки

(2.77)

(2.78)

Подставляем (2.71) в (2.76) и выделяем слагаемые с одинаковыми степенями п'. Для нулевого и первого порядка получим:

— Ө'(т)+ Н (<'(Т)..,.) = °. щп'— Ө'^^п'+ ()оп'+ Ө'' ()0п' = 0.

Здесь и в дальнейшем символ (F)0 означает значение функции F на луче <' = <'(т).

Ө'^ имеет мнимую часть. ТтӨ'^- положительно определенная квадратичная форма. Приравниваем мнимую часть (2.78) к 0:

d<' / дН\ *

dr \дӨу/0

Это равенство выполнено для любого nj. Полагая

п' = d ) _ df

\ дӨ^/0 dr

и учитывая положительную определенность 1тӨ(,) = 7тӨ'^'n'nj, получаем, что

7тӨ'^n'nj = 0 эквивалентно соотношению:

d<' = / дН)

dr дӨп^ 0.

Из (2.80) следует, что в (2.78) 2 слагаемых сокращаются и остается:

ddi = _ d дН)

dT \ д^7 0.

7тӨ'^'

nj = 0.

(2.79)

(2.80)

(2.81)

(-1

В формуле (2.81) знаем значение (dj (ибо значения всех функций на ПВ-луче известны).

Следовательно, Ө^' находится квадратурой: Ө^' = — J (ddj) dr. Подставляем Ө^' в (2.77) и

находим Ө(0). Ө(1) также находится ибо Ө(1) = Ө^'nj.

32

2.4.3 Второй порядок

Найдем теперь dij. Удовлетворение равенства (2.76) во втором порядке эквивалентно со-

отношению:

У"П + 2 (аДҢт)0 ө"'Cj + 1 (аҖ?)0 #'jn'n"+ +1 "ҖЩ)0ө"+ 2")0=0

(2.82)

Введем Г = (Г^), r^j = Г^^, матрицу квадратичной формы Ө(2) = 2.

17 d2H 4 ^'j + a2^\ ^i' + a2^ = 0

2 ^ag^ad^J0 ө + 2 ^аө^, ag^0 ө + 2 ^ag^ag^0 0

матричное уравнение Риккати:

) Ө'' dmj +

Тогда (2.82) можно записать как

dГ

— = ГАГ + ВГ + ГВT + C. dT

где A = -^ а2н 4 B = _7_^и_4 C =_7

где A аө^^^, В ag^ad^^^, C ag^ag^0.

Положим

(2.83)

(2.84)

Г _____ 2 dd j + 1

Г ij = 2 dT + 2

Г = ZY-1

(2.85)

Пусть Z и Y - решения следующей системы:

J dT = BZ + CY [ dT = -AZ - BTY

Тогда Г удовлетворяет уравнению (2.84). Проверим это.

Продифференцируем тождество Y-1Y = I по т:

dY-1 Y + Y-1 dY О <-.Дгчт.ттс trcv<-.TTT;rnT dY-1 .

d^ dT

(2.86)

= 0, отсюда находим dT :

dY-1 = -Y-1 dYY-1 dr dr

(2.87)

dr =

dT

= ГАГ + В Г + ГВТ + C

Потребуем, чтобы матрицы Z и Y удовлетворяли двум соотношениям:

dz Y-i + z dY__1 dT dT

= BZY-1 + C + ZY-1AZY-1 + ZY-1BT =

(2.88)

YTZ - ZTY = 0,

YZ*Y =

(2.89)

(2.90)

(I - единичная матрица, - знак эрмитова сопряжения).

Условия (2.89) и (2.90) обеспечивают симметричность Г и положительную определенность 1тГ. Кроме того, равенство (2.90) обеспечивает существование Y-1, то есть неравенство

det Y = 0. Доказывается от противного. Действительно, если det Y = 0 при каком-нибудь т, d1

то найдется вектор D =

1 (Y*Z - Z*Y) D:

= 0 такой, что YD = 0. Скалярно умножим вектор D на

у dm у

(в, 1 (Y*Z - Z*Y) = (D, D) > 0

(2.91)

33

С другой стороны:

D,1 (Y*Z - Z*Y) D ) = 1 (D, Y*ZD) - 1 (D, Z*YD) = 1 (YD, ZD) = 0 (2.92)

i / i i i

(так как YD = 0).

Мы пришли к противоречию, доказывающему, что равенство det Y = 0 невозможно.

Равенство (2.89) обеспечивает симметричность Г:

Г - Гт = ZY-1 - (YT)-1ZT = (YT)-1(YTZ - ZTY)Y-1 = 0,

равенство (2.90) - положительную определенность 1тГ:

1тГ = 1 (Г - Г*) = ^ (ZY-1 - (Y*)-1Z*) = ^(Y*)-1(Y*Z - Z*Y)Y-1 = 1 (YY*)-1. Таким образом, квадратичная форма Ө(2) = 1 (Гр,р) удовлетворяет всем нужным условиям.

Решая систему (2.86), находим Z и Y, а через них выражаем Ө(2).

Матрицы Z и Y удовлетворяют вдоль луча l линейной системе, являющейся системой уравнений в вариациях для канонической системы (2.58), при этом оказывается:

У Z = ().

ОЩ/ уощ/

Действительно, канонические уравнения для ПВ-лучей имеют вид:

0H ddi 0H

= ау, -аӯ (2'93)

Рассмотрим малую окрестность ПВ-луча: Д = ^:(т) + р:, = Ө:(т) + p:, где ф и p:- малые

величины. Раскладывая правые части (2.93) в ряд Тейлора по ф и p:, получим в первом

порядке уравнения в вариациях:

Di = ^_Р2^^ р.

dr 0 п + 0

dpi ^7_d^ p.

dT у д^д^/ 0 n уд$гдө^/ 0 p.

(2.94)

В окрестности ПВ-луча перейдем к координатам s, a1, a2. Разложим и p: в ряд по a1, a2:

Pi

= ^- ai oai

dp: ai oai

дп:

О— а2 + да2

, dp:

Д— a2 + 0^2

(2.95)

(2.96)

Подставляем и p: в систему (2.94) и приравниваем коэффициенты при одинаковых степе-

нях a:. Тогда в первом порядке получим:

7 ^7дпД ] dT у dak / ] ^7 дрг \ dT у дак у

д2^ 4 dnj + дӨ^д^Д/ о dak

7 д2^ 4 dnj ^д$гд$Д< 0 dak

д2Н 4 др^ дӨ^дӨ^/ о дак

7 д2^ 4 др^ У д^дӨ^Д 0 дак

(2.97)

Если ввести матрицы 4Д = , Z:k = = дОг, то система (2.97) приобретает вид системы

(2.86).

Итак, показали, что решения системы (2.86) можно записать в виде:

ДУ Z =(ад.

oak/ у oak/

(2.98)

34

Уравнения для Ө(г), r > 2 оказываются линейными. Приравнивая к 0 члены r-го порядка

в равенстве (2.76), получаем:

2.4.4 Расчет высших приближений

Э2Н \

Э2Н \

ЭӨ^^ ЭӨ^^ 0

ЭӨ(г)

"дГ +

n-

ЭӨ(г) dnj

ЭӨ(г) ЭӨ(2) dn- dn-

(2.99)

где Sr- известная функция, содержащая элементы, определенные на предыдущих этапах.

С учетом обозначений (4.16) уравнение (2.99) записывается в виде:

ЭӨ(г)

(W(2), AW(r)) - (n,BW(r)) = Sr

(2.100)

Здесь скобки - символ скалярного произведения без комплексного сопряжения bj = aj bj.

Знаем, что Ө(2) = 1 (Гщп) = 1 (ZY-1n,n). Определим VӨ(2):

(^Ө(2^) k = A (ZY-')j.n-nj = (ZY-'),..n-^*-k + (ZYnjZ-k = = 2(ZY-1 )k. n- = 2(ZY-1n)k

(2.101)

Принимая во внимание (2.101) и симметричность матрицы A, левую часть равенства (2.100) можно переписать в виде:

ддТ! - ((AZY-1 + BT) n, W(r)) = - ((AZ + BTY) Y-1n, ^Ө(г)) =

= + (-YY-'n. v^)

(2.102)

Здесь мы использовали выражение для -—г.

Итак,

ЭӨ(г) /dY

Ҷ- n

'Г

(2.103)

В выражении для Ө(г) от переменных (т, n) переходим к переменным (тщ), где = Y ^n. Тогда Ө(г) = Ө(г)(т,n) = Ө(г)(т, YY-1 n) = Ө(г)(т, Y^).

ЭӨ(г) (т, Y (т )^) эт

ЭӨ(г)(т, n) ЭӨ(г) dn-

dn- Эт

Эт

дӨ(r)(т,n) ЭӨ(гЬ^ dY 3

(n)+-

Эт

dn-

(2.104)

Таким образом, уравнение (2.103) теперь запишется в виде:

ЭӨ(г)(т, Y (т )^)

дт

'Т

(2.105)

Производная по т в равенстве (2.105) берется при фиксированном ^. Ө(г) находится квадратурой по т от известного выражения Sr. Таким образом, мы нашли асимптотический ряд для Ө.

2.4.5 Заключение

Перейдем к уравнениям переноса. В работах [1], [7] интегрирование уравнений переноса

в классе формальных степенных рядов основывалось на понятии квазисопряженных реше-

ний. Здесь мы будем пользоваться более простой техникой, еще более близкой к технике

стандартного лучевого метода.

35

Положения анализа, использующиеся для решения рекуррентных уравнений переноса, справедливы, если функции заменить формальными степенными разложениями. При этом все построения параграфов 3.2-3.3 проходят фактически без изменений. На этом пути для главного члена получаем формулу:

Список литературы диссертационного исследования кандидат наук Попов, Антон Игоревич, 2016 год

Список литературы

1. Бабич, В.М. Квазифотоны и пространственно-временной лучевой метод [Текст] / В.М. Бабич // Записки научных семинаров ПОМИ. - Санкт-Петербург, 2008. - том 342. - С. 5-13.

2. Бабич, В.М. Асимптотические методы в задачах дифракции коротких волн [Текст] / В.М. Бабич, В.С. Булдырев // Москва: Наука. - 1972. - 456 с.

3. Бабич, В.М. Пространственно-временной лучевой метод: линейные и нелинейные волны [Текст] / В.М. Бабич, В.С. Булдырев, И.А. Молотков // Ленинград: Издательство ЛГУ.

- 1985. - 272 с.

4. Бабич, В.М. Построение асимптотики решения уравнения Шредингера, сосредоточенной в окрестности классической траектории [Текст] / В.М. Бабич, Ю.П. Данилов // Записки научных семинаров ЛОМИ. - Ленинград, 1969. - том 15. - С. 47-65.

5. Бабич, В.М. Асимптотическое решение уравнения Гамильтона-Якоби, сосредоточенное вблизи поверхности [Текст] / В.М. Бабич, А.И. Попов // Записки научных семинаров ПОМИ. - Санкт-Петербург, 2011. - том 393 - С. 23-28.

6. Бабич, В.М. Квазифотоны волн на поверхности тяжелой жидкости [Текст] / В.М. Бабич, А.И. Попов // Записки научных семинаров ПОМИ. - Санкт-Петербург, 2010. - том 379

- С. 5-23.

7. Бабич, В.М. О квазифотонах волн Релея [Текст] / В.М. Бабич, А.В. Поскряков // Записки научных семинаров ПОМИ. - Санкт-Петербург, 2006. - том 332 - С. 7-18.

8. Бабич, В.М. Комплексный пространственно-временной лучевой метод и квазифотоны [Текст] / В.М. Бабич, В.В. Улин // Записки научных семинаров ЛОМИ. - Ленинград, 1981. - том 117 - С. 5-11.

9. Белов, В.В. Квазиклассическое траекторно-когерентное приближение для уравнений типа Хартри [Текст] / В.В. Белов, А.Ю. Трифонов, А.В. Шаповалов // Теоретическая и математическая физика, 2002. - том 130, Выпуск 8. - С. 460-492.

10. Бохнер, С. Функции многих комплексных переменных [Текст] / С. Бохнер, У. Т. Мартин ; пер. с англ. Б. А. Фукса. - Москва: Изд-во иностр. лит., 1951. - 300 с.

11. Виноградов, И.М. Математическая энциклопедия [Текст] / И.М. Виноградов // Советская энциклопедия. Москва: Наука. - том 4. - 1984. - 608 с.

12. Гуло, М.М. Н.А. Умов [Текст] / М.М. Гуло // Москва: Просвещение. - 1977.

13. Качалов, А.П. Нестационарные электромагнитные гауссовы пучки в неоднородной анизотропной среде [Текст] / А.П. Качалов // Записки научных семинаров ПОМИ. - Санкт-Петербург, 2000. - том 264. - С. 83-100.

81

14. Качалов, А.П. Система координат при описании квазифотонов [Текст] / А.П. Качалов // Записки научных семинаров ЛОМИ. - Ленинград, 1984. - том 140. - С. 73-76.

15. Кирпичникова, Н.Я. Волны на поверхности тяжелой жидкости с точки зрения пространственно-временного лучевого метода и его модификаций (линейная теория) [Текст] / Н.Я. Кирпичникова // Записки научных семинаров ЛОМИ. - Ленинград, 1987.

- том 165. - С. 91-101.

16. Маслов, В.П. Комплексный метод ВКБ в нелинейных уравнениях [Текст] / В.П. Маслов // Москва: Наука. - 1977. - 584 с.

17. Попов, А.И. Волновые валы для волн на поверхности тяжелой жидкости [Текст] / А.И. Попов // Записки научных семинаров ПОМИ. - Санкт-Петербург, 2012. - том 409. - С. 151-175.

18. Седов, Л.И. Механика сплошной среды [Текст] / Л.И. Седов // Физматлит. - Москва: Наука, 1983. - том 1. - 528 с.

19. Albers, M. A local mesh refinement multigrid method for 3D convection problems with strongly variable viscosity [Text] / M. Albers // Journal of Computational Physics. - 2000.

- Vol. 160. - P. 126-150.

20. Babich, V.M. Quasiphotons of waves on the surface of the heavy liquid [Text] / V.M. Babich,

A. I. Popov // Journal of Mathematical Sciences. - 2011. - Vol. 173 (3). - P. 243-253.

21. Bagrov, V.G. The semiclassical trajectory-coherent approximation in quantum mechanics [Text] / V.G. Bagrov, V.V. Belov, A.Yu. Trifonov // New York: Ann. Phys. - 1996. - Vol. 246. - No. 2. - P. 231-290.

22. Blankenbach, B. A benchmark comparison for mantle convection codes [Text] / B. Blankenbach, F. Busse, U. Christensen, L. Cserepes, D. Gunkel, U. Hansen, H. Harder, G. Jarvis, M. Koch, G. Marquart, D. Moore, P. Olson, H. Schmeling, T. Schnaubelt // Geophysical Journal International. - 1989. - Vol. 98. - P. 23-38.

23. Buiter, S.J.H. The numerical sandbox: comparison of model results for a shortening and an extension experiment in Analogue and Numerical Modelling of Crustal-Scale Processes [Text] / S.J.H. Buiter, A.Y. Babeyko, S. Ellis, T.V. Gerya, B.J.P. Kaus, A. Kellner, G. Schreurs, Y. Yamada // London: Geological Society, Special Publications. - 2006. - Vol. 253. - P. 29-64.

24. Busse, F.H. 3 D convection at infinite Prandtl number in Cartesian geometry - a benchmark comparison [Text] / F.H. Busse, U. Christensen, R. Clever, L. Cserepes, C. Gable, E. Giannandrea, L. Guillon, G. Houseman, H.-C. Nataf, M. Ogawa, M. Parmentier, C. Sotin,

B. Travis // Geophys. Astrophys. Fluid Dyn. - 1993 - Vol. 75. - P. 39-59.

25. Crameri, F. A comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the 'sticky air' method [Text] / F. Crameri, H. Schmeling, G.J.

82

Golabek, T. Duretz, R. Orendt, S.J.H. Buiter, D.A. May, B.J.P. Kaus, T.V. Gerya, P.J. Tackley // Geophysical Journal International. - 2012. - Vol. 189. - P. 38-54.

26. Dabrowski, M. MILAMIN: MATLAB-based finite-element method solver for large problems [Text] / M. Dabrowski, M. Krotkiewski, D.W. Schmid // Geochemistry, Geophysics, Geosystems. - 2008. - Vol. 9. - Q04030.

27. Deubelbeiss, Y. Comparison of Eulerian and Lagrangian numerical techniques for the Stokes equations in the presence of strongly varying viscosity [Text] / B.J.P. Caus // Phys. Earth Planet. Inter. - 2008. - Vol. 171. - P. 92-111.

28. Duretz, T. Discretization errors and free surface stabilization in the finite difference and marker-in-cell method for applied geodynamics: A numerical study [Text] / T. Duretz, D.A. May, T.V. Gerya, P.J. Tackley // Geochemistry, Geophysics, Geosystems. - 2011. - Vol. 12(7).

- P. 1-26. - Q07004.

29. Gerya, T.V. Introduction to Numerical Geodynamic Modelling [Text] / T.V. Gerya // Cambridge: Cambridge University Press. - 2010. - 358 p.

30. Gerya, T.V. An adaptive staggered grid finite difference method for modeling geodynamic Stokes flows with strongly variable viscosity [Text] / T.V. Gerya, D.A. May, T. Duretz // Geochemistry, Geophysics, Geosystems. - 2013. - Vol. 14(4). - P. 1200-1225.

31. Gerya, T.V. Characteristic-based marker-in-cell method with conservative finite-difference schemes for modeling geological flows with strongly variable transport properties [Text] / T.V. Gerya, D.A. Yuen // Phys. Earth Planetary Interiors. - 2003. - Vol. 140. - P. 293-318.

32. Gerya, T.V. Robust characteristics method for modeling multiphase visco-elasto-plastic thermo-mechanical problems [Text] / T.V. Gerya, D.A. Yuen. // Phys. Earth Planetary Interiors. - 2007. - Vol. 163. - P. 83-105.

33. Gerya, T. Practical analytical solutions for benchmarking of 2-D and 3-D geodynamic Stokes problems with variable viscosity [Text] / T. Gerya, I.S. Lobanov, A.I. Popov, I.Yu. Popov,

S. I. Popov // Solid Earth. - 2014. - Vol. 5. - P. 461-476.

34. Gerya, T. Numerical approach to the Stokes problem with high contrasts in viscosity [Text] /

T. Gerya, I.S. Lobanov, A.I. Popov, I.Yu. Popov // Applied Mathematics and Computation.

- 2014. - Vol. 235. - P. 17-25.

35. Gerya, T. Benchmark solutions for nanoflows [Text] / T. Gerya, I.S. Lobanov, A.I. Popov, I.Yu. Popov // Nanosystems: Physics, Chemistry, Mathematics. - 2014. - Vol. 5(3). - P. 391-399.

36. Gerya, T. On the Stokes flow computation algorithm based on Woodbury formula [Text] / T. Gerya, I.S. Lobanov, A.I. Popov, I.Yu. Popov // Nanosystems: Physics, Chemistry, Mathematics. - 2015. - Vol. 6(1). - P. 140-145.

83

37. Practical analytical solutions for benchmarking of 2-D and 3-D geodynamic Stokes problems with variable viscosity [Text] / T. Gerya, I.S. Lobanov, A.I. Popov, I.Yu. Popov, S.I. Popov // Solid Earth Discussions. - 2013. - Vol. 5. - P. 2203-2281.

38. Hager, B.H. A simple global model of plane dynamics and mantle convection [Text] / B.H. Hager, R.J. O'Connell //J. Geophys. Res. - 1981. - Vol. 86. - P. 4843-4867.

39. Ismail-Zadeh, A. Computational methods in Geodynamics [Text] / A. Ismail-Zadeh, P.J. Tackley // Cambridge: Cambridge University Press. - 2010. - 332 p.

40. Kathalov, A.P. Inverse Boundary Spectral Problems [Text] / A.P. Kathalov, Ya.V. Kurylev, M. Lassas // Chapman and Hall: CRC Monogr. and Surv. in Pure and Appl. Math. - Boca Raton, 2001.- Vol. 123. - 290 p.

41. Kaus, B.J.P. Effects of elasticity on the Rayleigh-Taylor instability: implications for large-scale geodynamics [Text] / B.J.P. Kaus, T.W. Becker // Geophysical Journal International. - 2007. - Vol. 168. - P. 843-862.

42. Kaus, B.J.P. 3D Finite amplitude folding: implications for stress evolution during crustal and lithospheric deformation [Text] / B.J.P. Kaus, S.M. Schmalholz // Geophys. Research Lett. - 2006. - Vol. 33(14). - L14309.

43. J.B. Keller Surface waves on water of non-uniform depth [Text] / J.B. Keller // J. Fluid Mech. - 1958. - Vol. 6. - P. 607-614.

44. Moresi, L. The accuracy of finite element solutions of Stokes' flow with strongly varying viscosity [Text] / L. Moresi, S. Zhong, M. Gurnis // Phys. Earth Planet. Inter. - 1996. - Vol. 97 (1-4). - P. 83-94.

45. Moresi, L. A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials [Text] / L. Moresi, F. Dufour, H.-B. Mulhaus // J. Comput. Phys. - 2003. - Vol. 184. - P. 476-497.

46. Popov, A. SLIM3D: A tool for three-dimensional thermomechanical modeling of lithospheric deformation with elasto-visco-plastic rheology [Text] / A. Popov, S. Sobolev // Phys. Earth Planet. Inter. - 2008. - Vol. 171 (1-4). - P. 55-75.

47. Popov, A.I. High Frequency Wave Packet Concentrated near the moving Line on the Liquid Surface [Text] / A.I. Popov // Journal of Physics: Conference Series. - 2013. - Vol. 435.

48. Popov, A.I. Wave wall type solution for liquid surface waves [Text] / A.I. Popov // AIP Conference Proceedings. - 2012. - Vol. 1479. - P. 728-731.

49. Ralston, J. Gaussian beams and the propagation of singularities [Text] / J. Ralston // Englewood Cliffs: MAA Studies Math. - New York, 1982. - Vol. 21. - P. 206-248.

84

50. Ramberg, H. Instability of layered system in the field of gravity [Text] / H. Ramberg // Phys. Earth Planet. Inter. - 1968. - Vol. 1. - P. 448-474.

51. Revenaugh, J. Dynamic topography and gravity anomalies for fluid layers whose viscosity varies exponentially with depth [Text] / J. Revenaugh, B. Parsons // Geophysical Journal International. - 1987. - Vol. 90(2). - P. 349-368.

52. Schmeling, H. A benchmark comparison of spontaneous subduction models - Towards a free surface [Text] / H. Schmeling, A.Y. Babeyko, A. Enns, C. Faccenna, F. Funiciello, T.V. Gerya, G.J. Golabek, S. Grigull, B.J.P. Kaus, G. Morra, S.M. Schmalholz, J. van Hunen // Phys. Earth Planet. Inter. - 2008. - Vol. 171 (1-4). - P. 198-223.

53. Schmid, D.V. Analytical solutions for deformable elliptical inclusions in general shear [Text] / D.V. Schmid, Yu.Yu. Podladchikov // Geophysical Journal International. - 2003. - Vol. 155. - P. 269-288.

54. Stadler, G. The dynamics of plate tectonics and mantle flow: from local to global scales [Text] / G. Stadler, M. Gurnis, C. Burstedde, L.C. Wilcox, L. Alisic, O. Ghattas // Science. - 2010.

- Vol. 329. - P. 1033-1038.

55. Tackley, P.J. Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid [Text] / P.J. Tackley // Phys. Earth Planet. Inter. - 2008. - Vol. 171. - P. 7-18.

56. Tackley, P.J. Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations [Text] / P.J. Tackley, S.D. King // Geochemistry, Geophysics, Geosystems. - 2003. - Vol. 4(4).

57. Torrance, K. E. Thermal convection with large viscosity variations [Text] / K.E. Torrance, D.L. Turcotte // Journal of Fluid Mechanics. - 1971. - Vol. 47. - P. 113-125.

58. Turcotte, D. L. Geodynamics [Text] / D.L. Turcotte, G. Schubert // Cambridge: Cambridge University Press. - 2002.

59. van Keken, P.E. A comparison of methods for the modeling of thermochemical convection [Text] / P.E. van Keken, S.D. King, H. Schmeling, U.R. Christensen, D. Neumeister, M.-P. Doin // Journal Geophys. Research: Solid Earth. - 1997.- Vol. 102 (B10). - P. 22477-22495.

60. van Keken, P.E. A community benchmark for subduction zone modeling [Text] / P.E. van Keken, C. Currie, S.D. King, M.D. Behn, A. Cangleoncle, J. He, R.F. Katz, S.-C. Lin, E.M. Parmentier, M. Spiegelman, K. Wang // Phys. Earth Planet. Inter. - 2008. - Vol. 171 (1-4).

- P. 187-197.

61. Zhong, S. Analytic solutions for Stokes' flow with lateral variations in viscosity [Text] / S. Zhong // Geophysical Journal International. - 1996. - Vol. 124. - P. 18-28.

85

62. Zhong, S. The role of plates and temperature-dependent viscosity in phase change dynamics [Text] / S. Zhong, M. Gurnis // Journal of Geophysical Research. - 1994. - Vol. 99. - P. 15903-15917.

63. Zhong, S. A benchmark study on mantle convection in a 3-D spherical shell using CitcomS [Text] / S. Zhong, A. McNamara, E. Tan, L. Moresi, M. Gurnis // Geochemistry, Geophysics, Geosystems. - 2008. - Vol. 9(10).

86

6 Приложения

Приложение 1

Построение Ф0

Все ищем в виде формальных степенных рядов по п2. Рассмотрим нулевой порядок.

Задача для Фо имеет вид:

дФо dz

к2Фо = 0,

= 0.

z=—H

(6.1)

Естественно, что для Фо имеет место формула:

Фо = Ао(т,п*,п2) ' cosh(k(z + H)).

(6.2)

Здесь к2 - это формальный степенной ряд по п2 (формула (3.12)). С формальными степенными рядами мы будем производить только действия сложения и умножения. Поэтому с ними можно работать как с обычными функциями (см. [8]). Фо ищем в виде:

Фо = X Ф0(^ п1, т )(п2)' = ЕЕ Ф0, (z + H )"(п2)' = Е ^"(z + H)", (6.3)

i=0 ^=о "=о "=о

где ^" Y2^=о ФО, (п2)^ - формальный степенной ряд по (п2)3 Мы предположили, что ФО раскладывается в степенной ряд по (z + H). Докажем это предположение. Будем искать решение Фо в виде формального степенного ряда по п2 и покажем, что коэффициенты этого разложения ФО выражаются только через функции, представимые степенными рядами по (z + H).

k = Е а (п2 )j, (6.4)

^=о

к2 = Е(п2)' Е а а^-1, (6.5)

^=О ^'=О

Фо = Е ФО(п2)', (6.6)

^=О

m m—i

к2Фо = Е(п2)т Е ФО Е "j"m—i—j . (6.7)

т=О 7=О ^'=О

Подставляем Фо и к2ФО в систему (6.1) и приравниваем коэффициенты при одинаковых степенях п2. Коэффициент при (п2)т:

d 2ФОт

dz2

m—1 m—i

аО фm = Е ФО Е aj am—i—j

i=0 ^'=О

дФ(^

dz

= 0,

z=—H

(6.8)

(6.9)

87

фт - функция, а не формальный степенной ряд, поэтому решение задачи (6.8)-(6.9) трудно

сти не представляет.

Замечание.

Известно, что уравнение вида y" — k2y = f (x) имеет решение:

1

y = C1cosh(kx) + C2sinh(kx) + — / f(^)sinhk(x — )d^. k Х0

(6.10)

В нашем случае, с учетом краевого условия (6.9), уравнение (6.8) имеет решение:

Ф^ = Amcosh(ao(z + H))

aj am-i-j sinh(ao(z — ())d(.

(6.11)

Последовательно решаем уравнения для Ф^ для различных m. Пользуемся формулой (6.11).

Ф° = A0cosh(a0(z + H)).

(6.12)

При дальнейшем применении формулы (6.11) в выражении для Ф^* может получиться либо sinh(...), либо cosh(...), либо (z + H)nsinh(...), либо (z + H)ncosh(...), либо их линейная комбинация. Все эти функции - это сходящиеся на всей оси степенные ряды по (z + H). Итак, доказали, что Ф° раскладывается в степенной ряд по (z + H). Обосновали формулу (6.3). Продолжим построения.

д 2Фо

dz2

= n(n — 1)(z + H)n-2 =

n=2

^n+2(n + 2)(n + 1)(z + H)n.

n=0

(6.13)

Подставляя (6.3) и (6.13) в (6.1), получаем:

д2ф

д 2---^2фо (^n+2(n + 2)(n + 1) — ^2^^(z + H)n = О,

(6.14)

откуда находим рекуррентное соотношение для коэффициентов:

^n+2(n + 2)(n +1) — k2^n = 0 ^n+2 =

k2

(n + 2)(n + 1)

(6.15)

Краевое условие (3.16) приводит к соотношению:

д Фо dz

z=-H

0 = ^nn(z + H)n-1

n=1

= ^1-

z=-H

Учитывая (6.15), получаем:

^2n+1 = 0,

k2

^2n = 2n(2n — 1) ^2n-2 = -

k2n

(2П)! ^0-

Таким образом,

ф0 ^0

n=0

k2n (2n)!

(z + H)2n = ^0 - cosh(k(z + H)).

(6.16)

(6.17)

(6.18)

(6.19)

Итак, решение Ф0 в классе формальных степенных рядов имеет тот же вид, что и в классе

функций:

Ф0 = А0(т, у1, у2) - cosh(k(z + H)).

88

Приложение 2

Теорема 1. Рассмотрим оДнороДную заДачу Штурма-Лиувилля:

Ту = (р(Ду')' - q(x)у = е

^оУ = (аоУ + ^УЭ = 0, (6.20)

^1У = (воУ + Ду') = 0.

p(x) > 0, Imq = 0. (6.21)

aj, Д, j = 0,1 - вещественные числа. p(x), q(x) - формальные степенные ряДы (ф.с.р.).

Пусть существует решение в виДе ф.с.р. уо = 0. ТогДа необходимое и Достаточное условие существования решения в виДе ф.с.р неоднородной заДачи Штурма-Лиувилля

Ту = -F,

< Ду = A,

^iy = в.

(6.22)

таково:

P(xi) ВУо

Pi

p(xo) —уо

ai

F (x)y0(x)dx.

(6.23)

ЗДесь F(x), A, B - ф.с.р.

Доказательство. Необходимость.

(p(x)y')'- q(x)y = -F-

(6.24)

Домножим уравнение (6.24) на уо и проинтегрируем по x от хо до xi:

/ ((РУ')' - ДУ)Уоdx / (ру')'Уodx Д^ УУУо^х =

J Х0 И Х0 И X0

ЛЖ1 ЛЖ1

- / ру'уоdx Д qууodx = У Х0 У X0

Д у(руо)'dx Д qууodx.

Х0 х0 х0

Х0

Х1 Л Х1

= РУ'Уо

= РУ'Уо

X1

- руо у

Х0

Х0

Х1

Х0

Подставим в (6.25) граничные условия:

уо (хо) = - ао Уo(xo),

ai

У'(хо) = — - аоУ(xo), ai ai

Уо (xi) = - в° Уо(хХ

Pi

у'(x1) = -B - У(x1).

pi pi

(6.25)

(6.26)

(6.27)

(6.28)

(6.29)

89

Тогда

/ ((py')'— qy)yodx = p(xi)B-yo(xi) — p(xi)yyo(xi) — Jx) ei ei

—p(xo)—yo(xo) + p(xo)ao yyo(xo) + p(xi) yoy(xi) — ai ai ei

ao /^i

—p(xo) yoy(xo)^ y((pyo)' — qyo)dx = ai </xo

С другой стороны,

= p(xi) Byo

Pi

— p(xo) -Ayo ai X=X1 1

(6.30)

'Ж1

X=Xo

FX1 PX1

/ ((py')'— qy)yodx = — / Ғуо^ж.

Jxo </Xo

(6.31)

</Ж0 </Ж0

Из сравнения (6.30) и (6.31) получаем формулу (6.23), то есть необходимость доказана. Достаточность.

Пусть ф - решение задачи Коши:

(pO' — q^ = — ғ,

X=Xo

(6.32)

ф'

= B.

= A

X=Xo

Задача Коши всегда разрешима. Следовательно, ф существует. Пусть y = ф — yo.

(py')'— qy = (p^')'— q^ — (pyo)' + qyo = -С

(6.33)

то есть y удовлетворяет нужному уравнению. Проверим краевые условия. Домножим первое уравнение в (6.32) на yo и проинтегрируем по x от xo до xi:

((p^')' — q^)yo dx =

(pO'yodx

q^yodx =

Xi

= p^'yo

Xo

p^'yo dx

q^yodx =

= p(^'yo — yoф)

xi /<X1

^((pyo)'— qyo)dx =

Xo AXo

= p((y' + yo )yo — yo (y + yo))

Xi

Xo

= p(y' yo — yo y)

X1

Xo

(6.34)

то есть

((p^')' — q^)yodx = p(xi)y'(xi)yo(xi) — p(xi)yo (xi)y(xi) —

—p(xo)y'(xo)yo(xo) + p(xo)yo (xo)y(xo).

(6.35)

Подставляя граничные условия (6.28)-(6.29) в (6.35), приходим к формуле:

90

/*Х1 во

/ ((p^')' - q^)y0dx = p(xi)y'(xi)y0(xi) + p(xi)в0y0(xi)y(xi) -

Ухо Pi

-p(x0)y'(x0)y0(x0) - p(x0)a0y0(x0)y(x0). ai

(6.36)

С другой стороны,

^X1 ҒЖ1

/ ((py')'- qy)y0dx = - / Fy0dx.

Jxo ^xo

Х1

(6.37)

хо

Пусть

p(xi) By0

Pi

Тогда из (6.36)-(6.38) следует, что:

- p(x0) ^Ay0

ai

X=X1 i

Х=Хо

M1

= - / F(x)y0(x)dx.

Ухо

(6.38)

p(xi) By0

Pi

- p(x0) ^Ay0

X=X1 ai

X=Xo

= p(xi)y'(xi)y0(xi) + p(xi)p0y0(xi)y(xi) -

Pi

-p(x0)y'(x0)y0(x0) - p(x0)a0y0(x0)y(x0).

ai

(6.39)

Условие (6.39) должно быть выполнено для любых точек x0 и xi. Зафиксируем x0 и будем менять xi. Так как соотношение (6.39) должно быть выполнено всегда, то части уравнения, отвечающие x0 и xi должны быть независимы друг от друга:

-p(x0)—y0(x0) = -p(x0)y'(x0)y0(x0) - p(x0)a0y0(x0)y(x0). ai ai

(6.40)

p(xi) By0 (xi) = p(xi)y'(xi)y0(xi) + p(xi) p0 У0 (xi)y(xi). Pi Pi

y0(xi) = 0. Доказывается от противного. Если y0(xi) = 0, то из краевого условия:

Ду = (в0У + Piy')

= 0 следует, что: y0

Х = Х1

(6.41)

= 0. Значит y0 - решение задачи Коши:

X=X1

(p(x)y0)' - q(x)yo = Ф

Уо

= 0,

Х = Х1

у0

= 0.

Х = Х1

Следовательно, y0 = 0, что противоречит предположению теоремы. По условиюp > 0, Pi = °. Поделим обе части (6.41) на p(xi)y0(xi) и домножим на Pi, получим:

B = Piy'(xi) + в0У(х1),

(6.42)

то есть получили нужное условие при x = xi. Вернемся к соотношению (6.40):

-p(x0) —У0(Х0) = -p(x0)y'(x0)y0(x0) - p(x0) — У0(Х0)y(x0). ai ai

(6.43)

91

Учитывая, что p > О,

условие при x = x0:

y0(x0) = 0 (доказывается аналогично y0(x1)), a1 = 0, получаем нужное

a1y'(x0) + a0y(x0) = A-

(6.44)

92

Приложение 3

Предварительные сведения. Рассмотрим r + s формальных степенных ряДов.

X Ai(x)y" = a(a', lal=o j = 1, ^, (6.45)

X Aa(x)ya = t(a-r', j= = r + 1,r + 2, ,r + s . (6.46)

a = (ai, a2, ' ,as) , (6.47)

]a] = ai + a2 + + as, ai G Z+ . (6.48)

Заданные коэффициенты Aa G в некоторой области Q G Rr.

Переменные x = (xi, x2, , xr) буДем называть конечными переменными, y = (yi, y2, , ys) - малыми переменными.

Пусть x и y - формальные степенные ряДы:

x(a' = X Bi Mb", j =1,2.........r. (6.49)

lal=o

y(a-r' = X Ba (a)&", j = r + 1,r + 2, , r + s. (6.50)

H=i

Ba (a) G C^(Qi) (Qi G Rr) - функции, причем область значений B"oo o'(a) G Q. a =

(a(i',a(2', , a(r') - конечные переменные. b = (b(i',b(2', , b(s') - малые переменные.

ТогДа xa и ya, Даваемые формулами (6.49), (6.50) можно поДставлять в ряДы (6.45), (6.46). При подстановке формальных степенных ряДов Для xa и ya функции Aa(x) заменяются ряДами Тейлора-Маклорена.

СхоДимость не предполагается.

Определение. Обращение ряДов (6.45), (6.46) - это нахождение таких ряДов (6.49), (6.50), что при их подстановке в (6.45), (6.46) получится:

X^Aa(x)ya = a(a', j = 1,2, ,r, lal=o (6.51)

X Aa (x)ya = b(a-r', j = r + 1, r + 2, , r + s . l"l=i (6.52)

Теорема 2. Рассмотрим ряДы:

X Aa(x)ya = a(a\ j = 1,2, lal=o (6.53)

X Aa (x)ya = b(a-r', j = r + 1, r + 2, , r + s . H = i (6.54)

93

J =

Aa G СД(П), Q c Rr .

Пусть определитель

D(A1„, ,)(x)......A[,, ,)^Aa+1y".....^Aa+sy")

lai=i

H=i

гДе