Аналитические решения уравнений Стокса в криволинейных координатах для тестирования численных алгоритмов тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Макеев Илья Владимирович

• Макеев Илья Владимирович
• кандидат наук
• 2016

Оглавление диссертации кандидат наук Макеев Илья Владимирович

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

3.1.2 Построение решений для течений над полостью

3.1.3 Картина линии тока для течений над полостью

3.2 Течение в полости, имеющей форму эллиптической подковы

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

3.2.2 Построение решений

3.2.3 Геометрическая картина течений

3.3 Проверка конечно-разностной схемы на эталонных решениях

3.4 Выводы

Заключение

Литература

Приложение

Приложение

Введение диссертации (часть автореферата) на тему «Аналитические решения уравнений Стокса в криволинейных координатах для тестирования численных алгоритмов»

Введение

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

Нахождение решений уравнений Стокса является актуальной задачей для различных областей науки. В задачах геодинамики уравнениями Стокса описывается течение Земной мантии [3,4]. В ряде исследований система уравнений Стокса дополняется уравнением переноса тепла и используется для описания процесса конвекции [5,6]. Уравнения Стокса используются при моделировании течения крови [7], течений в наноструктурах [8,9], при описании течений жидкости, возникающих в процессе химического травления [10]. Стоксовы течения исследуются в седиментационном анализе при определении диаметра частиц по скорости их оседания, а также в гидродинамической теории смазки [11].

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

соответствующей криволинейной системы координат. В работе [15] найдены решения для круговой полости в полярной системе координат, в [1 6] решения уравнений Стокса получены в эллиптической системе координат для области, ограниченной двумя эллипсами. Аналогичный подход применяется для трехмерных течений в тех случаях, когда граница имеет сферическую или цилиндрическую форму. В [17,18] находятся аналитические решения уравнений Стокса в сферической системе координат, в [19-21] уравнения решаются в цилиндрических координатах.

Наиболее распространенным подходом, используемым при построении аналитических решений уравнений Стокса, является метод разделения переменных. В ряде исследований применяется метод функций Грина [2,22]. В некоторых работах для нахождения аналитического решения используют метод декомпозиции [23,24]. Область разбивается на регионы, в каждом из которых задается общий вид функции решения. Для решений в регионах ставятся условия склейки, после удовлетворения которых находится решение для всей области.

Решения уравнений Стокса для систем, возникающих при моделировании физических процессов, как правило, не могут быть получены в явном виде. В таких случаях решения находят численно. Множество исследований было посвящено разработке эффективных алгоритмов численного решения уравнений Стокса. Наибольшее распространение получили алгоритмы, основанные на конечно-разностных схемах [25-27], алгоритмы, использующие метод конечных элементов [28,29] и метод конечных объемов [30-32].

Одним из подходов, позволяющих сократить время сходимости численных алгоритмов, является применение многосеточных методов [33-35]. Многосеточные алгоритмы для решения уравнений Стокса, основанные на методе конечных элементов, рассматриваются в статьях [36-39]. В [3,40] представлены многосеточные алгоритмы, использующие конечно-разностные схемы. В [41] описываются многосеточный метод для решения уравнений Стокса, основанный на алгебраическом подходе.

Для некоторых задач геодинамики возникает необходимость находить решения уравнений Стокса с учетом неоднородной вязкости среды. Вязкость горных пород, как правило, зависит от их температуры. Если температура верхних и нижних слоев среды значительно отличается, то моделирование течений для таких систем приводит к численному решению уравнений при высоких контрастах вязкости [42-45]. Разработке эффективных алгоритмов для подобных задач посвящен ряд исследований [46-52]. Значительное влияние на точность численного решения оказывает конфигурация сетки. Если скачок в значениях вязкости среды происходит не на границе ячейки сетки, а внутри нее, это может приводить к росту ошибки при построении численного решения. Одним из подходов к решению проблемы является предварительный учет распределения вязкости среды при генерации ячеек сетки и подбор разного разрешения сетки для участков области с различными перепадами вязкости [53,54]. Интересно, что такого же типа стоксовы течения с высокими контрастами вязкости наблюдаются и в совершенно других пространственных масштабах -при описании течений в наноструктурах [55-57].

Ряд задач геодинамики, связанных с моделированием течений во внутренних областях планет, приводит к необходимости решать уравнения Стокса в сферической полости [4] (чаще всего рассматривается область, ограниченная двумя концентрическими сферами).

Естественным шагом к построению решений уравнений Стокса для областей со сферически симметричной границей становится переход из декартовой в сферическую систему координат. Для численного решения уравнений в таких областях было разработано множество сеток. В статьях [58,59] строятся численные решения на ортогональных сетках, ребрами ячеек которых становятся параллели и меридианы сферы.

Решения уравнений на ортогональных сетках типа "YmYang" рассматриваются в [60-62].

Ряд сеток был получен путем проекции правильных многогранников на поверхность сферы. В [63] такие сетки строятся для икосаэдра, в [64-66] для куба,

в [67, 68] для октаэдра, в [69] для тетраэдра. Сетки, полученные таким путем, не являются ортогональными, но хорошо подходят, например, для использования метода конечных элементов.

В некоторых работах при моделировании течения мантии в сферической полости применяется спектральный метод построения решения [70-72].

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

Помимо сферической системы координат для ряда задач оказывается удобным находить численное решение уравнений Стокса, используя цилиндрическую систему. Течения в некоторых задачах геофизики рассматриваются в цилиндрической системе координат [6,75]. При моделировании течения крови в сосудах рассматриваются осесимметричные течения, описываемые уравнениями Стокса [7]. Некоторые программы, предназначенные для численного решения уравнений Стокса, имеют модификации для цилиндрической формы границы [76]. Уравнения Стокса, в том числе с переменной вязкостью, активно используются для моделирования течений в наноструктурах [21], особенно в моделях, основанных на кристаллитной теории жидкости [77].

Важным этапом разработки программ численного решения уравнений является процесс тестирования. Результатом тестирования может стать оценка корректности работы алгоритма. Одним из способов тестирования является сравнение результата работы численного алгоритма с аналитическим решением, полученным для некоторой тестовой задачи. Аналитическое решение в таком случае выступает в качестве эталона. Тестовая задача может быть специализирована для проверки корректной работы программы при варьировании

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

Ряд исследований в области геодинамики был посвящен нахождению аналитических решений, ставших основой для создания тестовых задач, используемых для проверки численных алгоритмов. В статьях [22,53] рассматриваются аналитические решения уравнений Стокса для течений с переменной вязкостью в двумерных прямоугольных областях. Эталонные аналитические решения для трехмерных течений с переменной вязкостью выводятся в [78]. В [79] представлено аналитическое решение уравнений Стокса для течения в бесконечном канале. Аналитическое решение, задающее стоксово течение с переменной плотностью, рассматривается в [80].

Аналитическое решение уравнений Стокса с постоянной вязкостью, описывающее двумерное течение над полостью, выводится в [24]. Такие решения находят применение при тестировании численных методов решения уравнений Стокса для двумерных областей [10,40,81,82].

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

Другим способом тестирования является сравнение результатов работы нескольких численных алгоритмов между собой при решении некоторой общей для них задачи. В работах [84,85] рассматривается тестовая задача для процесса субдукции. В статьях [86,87] производится сравнение численных алгоритмов при

моделировании процесса термальной конвекции.

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

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

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

Научная новизна результатов, полученных в работе:

1. Найдены частные аналитические решения уравнений Стокса с переменными вязкостью и плотностью в цилиндрической и сферической системах координат для таких параметров системы, при которых вязкость, плотность и напряженность

гравитационного поля являются функциями координаты г: ц = ?](г), р = р(г),

О = О(г).

2. Предложены и протестированы численные схемы для многосеточного метода решения уравнений Стокса с переменной вязкостью на ортогональных смещенных сетках в цилиндрической и сферической системах координат. Написан комплекс программ в среде Ма1ЬаЬ.

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

4. Выведено явное решение краевых задач для уравнений Стокса и неразрывности для течений в области, имеющей форму эллиптической подковы.

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

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

Методы исследования, используемые в данной работе, можно разделить на две группы - аналитические и численные. К аналитическим методам относятся такие методы, как метод Фурье решения дифференциальных уравнений в частных производных, методы теории функций комплексной переменной, методы решения обыкновенных дифференциальных уравнений. К численным методам относятся метод Гаусса-Зейделя, используемый для решения систем линейных уравнений, многосеточный метод решения дифференциальных уравнений, а также методы и приемы написания программ в пакете MatLab.

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

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

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

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

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

Апробация работы. Основные результаты работы докладывались на международных конференциях «Mathematical Challenge of Quantum Transport in Nanosystems» (Санкт-Петербург, 2014 г. и 2015 г.), а также III и IV Всероссийском конгрессе молодых ученых (Санкт-Петербург, 2014 г. и 2015 г.).

Публикации. Основные результаты по теме диссертации изложены в 8 печатных работах, из них 3 в журналах из Перечня ВАК, 1 в журнале, индексируемом в Web of Science и Scopus.

1. Makeev I.V. A benchmark solution for 2D Stokes flow over cavity / I.Y.Popov, I.V. Makeev // Zeitschrift fur Angewandte Mathematik und Physik - 2014. - V. 65. - P. 339-348 (входит в Web of Science и Scopus).

2. Makeev I.V. Analytical benchmark solutions for nanotube flows with variable viscosity / I.V. Makeev, I.V. Blinova, I.Y. Popov // Nanosystems: Physics, Chemistry, Mathematics - 2015. - V. 6. - № 5. - P. 672-679 (входит в Перечень ВАК).

3. Макеев И.В. Эталонные решения уравнений Стокса с переменной вязкостью в цилиндрических и сферических координатах для тестирования алгоритмов / И.В. Макеев, И.Ю. Попов, И.В. Блинова // Научно-технический вестник информационных технологий, механики и оптики - 2016. - Т. 16. - №1.

- С. 161-167 (входит в Перечень ВАК).

4. Makeev I.V. Steady Stokes flow between confocal semi-ellipses / I.V. Makeev, I.Y. Popov // Nanosystems: Physics, Chemistry, Mathematics - 2016. - V. 7 - № 2.

- P. 324-331 (входит в Перечень ВАК).

5. Макеев И.В. Эталонные решения уравнений Стокса для осесимметричных течений / И.В. Макеев // Сборник трудов IV Всероссийского конгресса молодых ученых - СПб: Университет ИТМО - 2015. - С. 271-273.

6. Макеев И.В. Решения гидродинамических уравнений в сферической системе координат / И.В. Макеев // Труды студенческого центра прикладных математических исследований - СПб: Университет ИТМО. -2014. - С. 44-47.

7. Макеев И.В. Эталонные решения гидродинамических уравнений при малых числах Рейнольдса / И.В. Макеев // Сборник тезисов докладов конгресса молодых ученых, Выпуск 2. - СПб: Университет ИТМО - 2014. - С. 266.

8. Makeev I.V. Benchmark solutions for Stokes flow through nanostructures / I.V.Makeev // Mathematical Challenge of Quantum Transport in Nanosystems: Book of Abstracts of the International Conference (Saint Petersburg, September 9-11, 2015). -Saint Petersburg: ITMO University, 2015. - P. 15.

Построенные в диссертационной работе аналитические решения уравнений и разработанный комплекс программ внедрены и используются для тестирования алгоритмов расчета стоксовых течений в компании ООО «Корнинг СНГ».

Объем и структура работы. Диссертация состоит из оглавления, введения, трех глав, заключения, списка литературы и приложений. Полный объем диссертации составляет 103 страницы с 48 рисунками. Список литературы содержит 88 наименований.

Глава 1

Аналитические решения для уравнений Стокса в цилиндрических и сферических координатах

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

1. 1 Уравнения Стокса в цилиндрических координатах

^стема уравнений Навье-Стокса, используемая в гидродинамике для описания течения несжимаемой жидкости с переменной вязкостью и плотностью имеет вид [3]:

— ду- ду-

— а1к +рОг = р( + у, ), (1.1)

—хк 'к ' -I к дхк ( )

где с - тензор напряжений; р - плотность жидкости; О - напряженность гравитационного поля, у - вектор скорости.

Уравнения (1.1) дополняются уравнением неразрывности:

V(рv) = 0, (1.2)

При описании стационарного течения система (1.1) примет вид:

- -у.

Ск +РО =р(ук-г-) (1.3)

/-V I к I I I V к ^

дхк дхк

При малых числах Рейнольдса можно пренебречь членом в правой части (1.3). В этом случае система (1.3) сведется к системе уравнений Стокса [3]:

д

^ а = р (14)

В цилиндрических координатах (г ,ф, х) система уравнений Стокса примет вид:

=-рОг,

(У-*\=-рО,, (1.5)

(V-*) 2 = -рв,,

Уравнение неразрывности (1.2) можно переписать в форме:

ру^ + д(ру) | 1 д(рф | а(руг) = 0 г дг г дф дх

В координатах (г,ф, х) дивергенция тензора напряжений имеет следующие компоненты:

^.а)г = а ^Оп^ОФ^ 1 а + а,

г дг г г дф дх

до 2а 1 да да

(V.*) =_ф +_гФ +1_фф +_ф,

ф дг г г дф дх

^ л дагх а 1 дафх дахх

(V - а)г = —— + —- н---— н--—.

дг г г дф дх

Тогда уравнения системы (1.5) могут быть записаны в форме:

Т + +1 Т + Т-дР = -рОг, (1.6)

дг г г дф дх дг

т + т + 1 -1 дР = -Р09, (1.7)

дг г г дф дх г дф

н 1 Т + Т-дР = -рОх, (1.8)

дг г г дф дх дх

где т - девиаторный тензор напряжений, Р - давление.

В цилиндрических координатах т имеет следующее компоненты

~ ду„

Тг = ,

дг

г„ = + У-),

г -ф г

~ -у,

т„ = 2л —-

22 / /-ч

-2

Г -ф -Г г

-Г -2

1 -V

ф + - -Г2 х

-2 Г -ф'

где л - вязкость жидкости. В развернутой форме система (1.5) имеет вид:

„ -л -V „ -Ч „1 -V ~ 1 -уф „ у 1 -л -У„

2 — —^ + 2л—Г + 2л--^ - 3^ — —--2л+ ^ — — +

-Г -Г -Г Г -Г Г -ф Г Г -ф -ф

-л уф

1 - 2уГ 1

2 +л —^ + -л

Г -ф

- Ч

1 -л-^

Г -ф -Г -ф Г2 ' Г2 -ф2 Г ' -ф-Г

-л -у -л -у -у -у -Р

+—- —- + —- —- + л—- + л

-2 -Г -2 -2 -2-Г -2

2 ^ ^

(1.9)

- Ч

1 -л ду -л -Уф 1 -л 1 ^2у ---—L н---—ф---- у +л--— + л- 2

Г -Г -ф -Г -Г Г -Г Г -Г-ф -Г

1 -у 1 -У

3 л —т —L + л----л-- + 2

Г -ф Г -Г Г

ф 2± л-Уф , _

2 Г2 -ф -ф Г2 -ф

„1 -л

+ 2 — — У +

(1.10)

2л

1 - Уф -л -Уф 1 -л -Уг - Уф 1 -Ч 1 -Р

■ +

- + —

Г2 -ф2 -2 -2 Г -2 -ф

л—ф+л——--~— = -рОф

-2 Г -2-ф Г -ф

-л -уг -л -уг -2у -2У 1 -у2 1 -У 1 -л -у<

■ +

■+л

-Г -Г -Г -2 -Г

2

+ л

-Г-2

+ л

+л

Г -Г Г -2 Г -ф -2

+ -

1 -л -У 1 - Уф 1 -2У _ -2У _ -л -Уг -Р

■ + л

+л

Г -ф -ф Г -ф-2 Г -ф

+ 2л

-2 2

-2 -2 -2

= -рО2

(111)

2

1.2 Построение решения в осесимметричном случае

Решим задачу построения частных аналитических решений системы уравнений (1.2), (1.9)—(1.11) для таких параметров системы, при которых вязкость, плотность и напряженность гравитационного поля являются функциями координаты г : л = л(г), р = р(г), О = 0(г).

Будем искать компоненты скорости и давление в следующей форме: у = у (г), \ = \(г), ^ = ^(г), Р = Р(г). Тогда уравнения (1.9)—(1.11) преобразуются к виду:

2П-Ч + 2лЧ + 2л<- 2л\у - Р' = -рОг, (1.12)

г г

Л Ч;-~Л 'Ч; + ЛЧ"; + V;=-PG;, (113)

г г г

л-< + лЧ +лч, = -рОг. (1.14)

г

Уравнение (1.2) примет вид:

р1 у +р\ + ру'г = 0. (1.15)

г

Разделив переменные в (1.15), имеем:

г А р \л —L = -(— + — )йг.

Уг г р

Проинтегрировав, находим первую компоненту скорости:

с

Vг =—, (1.16)

гр

где с - произвольная константа.

Зная функцию для первой компоненты скорости (1.16), из уравнения (1.12) находим давление:

Р(г) = \ (рОг +2л - V + 2лЧ + 2лу"г - 2л1уг )&. (1.17)

¿г г

Путем элементарных преобразований уравнение (1.13) может быть приведено

к виду:

^чЛ1)-^^)= -р; (1.18)

Л г г л г л

Будем искать решения уравнения (1.18) в виде функций у (г) = (г) + ^ (г), где (г) является частным решением уравнения (1.18), а ^ (г) является общим решением однородного уравнения:

ч;+Ч;(Л+1) - Ч/ Л+-1)=0. (1.19)

Л г г л г

Выполним следующую подстановку в уравнении (1.19):

2 = ■

(1.20)

(1.21)

В результате получим уравнение первого порядка относительно 2 :

' л' \ 2 1 л' 1 Л 2 + 2(— + -) + 2-------- = 0.

л г г л г

Произведем следующую подстановку: 2 = — у, где у = у(г).

лг

Уравнение (1.21) преобразуется к уравнению Риккати:

у + у2 — -л'-- = 0. (1.22)

лг г

Заметим, что функция у = л является частным решением уравнения (1.22). Будем искать решение уравнения (1.22) в виде функции у = и+л, где и = и(г). Тогда функция и должна удовлетворять уравнению Бернулли:

, 01 1 2

и = -2 — и--и .

г лг

Далее следует учесть, что и ( г ) = 0 является решением уравнения (1.23).

(1.23)

Выполним следующую подстановку: и =—q, где а = а( г). Уравнение (1.23)

Г

преобразуется к виду:

а

лг

а2 = 0.

Проинтегрировав, находим:

а =

/-V

лг

Возвращаясь к функциям у , получим уравнение:

у лг г J лг

у

у

ф

1

1

Разделив переменные и выполнив интегрирование, находим:

Г 11 1

1п(| уф I) =1 (- + — --ж + С21.

ф 1 Г2 лГ2 2 1 ,

1-^1 + с„ 1 лг

Решение и(г) = 0 уравнения (1.23) приводит к следующему решению уравнения (1.19):

уф = СГ.

В результате, имеем общее решение однородного уравнения (1.19):

Уф1(г) = Сх/(Г) + С2Г ,

Г

Г 1 1 1

где /(г) = ехр( I (— + —--)ф2), с1, с2 - произвольные константы.

^ Г лг Г 1 12 /2 Г-Ц^ + С

лг

Частное решение уравнения (1.18) находим методом Лагранжа в виде

Уф0(г) = С,(г)/(г) + С2(г)г .

Здесь производные функций С1(г), С2( г) удовлетворяют следующей системе:

С! ' / (г) + С2 г = 0, С! '/'(Г) + С2 ' =-рОф .

л

Решением системы будут функции вида

С>(г) = ! *,

1 л(/ - г)

С2 (г) =

* л(/ - г)

В результате получаем следующее решение уравнения (1.13), задающее вторую компоненту скорости:

Уф (Г ) = СХ/(Г) + С2Г + С (г )/(г ) + С2 (г )г. (1.24)

Решим уравнение для третьей компоненты скорости. Пусть у = w(г), тогда уравнение (1.14) запишется в форме:

w' + (Л + = -рО.. (1.25)

Л г л

Выполним следующую подстановку: w = ш—, где ш = ш(г). Уравнение (1.25)

Лг

преобразуется к виду

1 р г ш — =--

Лг л

Проинтегрировав и выполнив обратную подстановку, находим функцию w :

1 г

-д ^ро^г+с1).

w =--(I г рО„аг + с

Лг *

В результате, находим решение уравнения (1.14) в виде:

г 2 г2

у = —(| г1р0^йг1 + сх)йг2 + с 2. (1.26)

1 Лг2 1

Формулы (1.16), (1.17), (1.24), (1.26) задают частные аналитические решения уравнений Стокса (1.5) и неразрывности (1.2).

Рассмотрим модификацию предыдущей задачи для случая, когда плотность среды зависит от первой и третьей координаты, а напряженность гравитационного поля равна нулю: р = р(г,,), л = л( г), о = 0.

Будем искать решения следующего типа: у = у (г), у, = у, (г), у = у (г, ^),

Р = Р(г).

Уравнение неразрывности (1.2) примет вид:

1 г ' ду, др _ {л

р-у +р'чг + ру'г +р+ у2 = 0. (1.27)

г д, д,

Предположим, что функция плотности представима в форме р = р (г)р (г). В этом случае уравнение (1.27) может быть преобразовано к виду:

Л 1 др(г) ,ду 1 дрЛ£) .

(- у +--^^ у + уг) = -(—- +--у).

г р(г) дг д, р(,) д,

Положим, что V = VI (гКг (г). Тогда, завершая процесс разделения переменных, имеем:

1 Л 1 др(г) .д^, (г) 1 да (г) . чч

-(- V +--^^ V + V ) = -( г24 ' +--^^ V, (г)) = Л.

/ч\ г / \ ^ г г/ \ ^ / \ ^ г 2 V //

VI (г) г а (г) дг дг р(г) дг

Для нахождения частного решения, рассмотрим случай Л = о:

^ 2(г) ^ ^ (г) = 0 (1.28)

дг р (г) дг

1 1 др(г) , (л

-V + —- Уг + уг = 0 (1.29)

г р (г) дг Проинтегрировав (1.28), имеем

V 2(г) = , (1.30)

Р2( г )

где й - произвольная константа.

Уравнение (1.29) совпадает с уравнением (1.15), решение которого было получено ранее:

с

гр1

(1.31)

Уравнение (1.11) преобразуется к виду:

' г —т + п--- + 2п—Т = 0. (1.32)

дг дг ' дг2 г дг дг2

Разделение переменных в уравнении (1.32) дает

1 .дп ду,(г) д2V,(г) 1 дг,(гУ „ 1 д2V,(г)

(^ > +п—+ п--2^2) = -2--= Л.

71ул (г) дг дг дг2 г дг гг2 (г) дг

Для частного случая Л = 0, имеем

д= 0, (1.33)

дг

дп <д$ Лт) д2г ,(г) 1 дv ,(г) ,л 0,ч

/ .и 7 + п—^^+п—= 0. (1.34)

дг дг дг г дг

Проинтегрировав (1.33), находим

V =

Уг 2 = аг + b,

(1.35)

где а, Ь - произвольные константы.

Уравнение (1.34) является частным случаем уравнения (1.14) при о = 0. Таким образом, решением уравнения (1.34) будут функции вида (1.26) при о = 0:

у,1 = с

Лг

В результате, получаем частное решение уравнения (1.32)

йг. (1.36)

= (а, + Ь) Г—йг. (1.37)

Лг

Уравнение (1.9) примет вид:

1 1 д ду 2л - у; + 2 Л Ч + 2лу"г - 2Л — Уг +Л — (—^) - Р' = 0. г г дг д,,

Список литературы диссертационного исследования кандидат наук Макеев Илья Владимирович, 2016 год

Литература

1. Слёзкин, Н.А. Динамика вязкой несжимаемой жидкости / Н.А. Слёзкин. -М.: Государственное издательство технико-теоретической литературы, 1955. - 521с.

2. Pozrikidis, C. Introduction to Theoretical and Computational Fluid Dynamics / C. Pozrikidis. - Second edition. - Oxford University Press, 2011. - 1296 p.

3. Gerya, T. Introduction to Numerical Geodynamic Modeling / T. Gerya. - Cambridge University Press, 2010. - 358 p.

4. Ismail-Zadeh, A. Computational Methods for Geodynamics / A. Ismail-Zadeh, P. Tackley. - Cambridge University Press, 2010. - 332 p.

5. Davies, D. R. A hierarchical mesh refinement technique for global 3-D spherical mantle convection modelling / D. R. Davies, J. H. Davies, P. C. Bollada, O. Hassan, K. Morgan, P. Nithiarasu // Geosci. Model Dev. - 2013. - V. 6. - P. 1095-1107.

6. Jarvis, G.T. Effects of curvature on two-dimensional models of mantle convection: Cylindrical polar coordinates/ G.T. Jarvis // J. Geophys. Res. - 1993. - V. 98.

- p. 4477-4485.

7. Keener, J. Mathematical Physiology II: Systems Physiology / J.Keener, J.Sneyd. -Second edition. - New York: Springer, 2009. - 580 p.

8. Popov, A.I. Benchmark solutions for nanoflows / A.I. Popov, I.S. Lobanov, I.Yu. Popov, T.V. Gerya // Nanosystems: Phys., Chem., Math. - 2014. - V. 5. - P. 391-399.

9. Kirby, B.J. Micro- and Nanoscale Fluid Mechanics: Transport in Microfluidic Devices / B.J. Kirby. - Cambridge University Press, 2010. - 512 p.

10. Driesen, C. H. An accurate boundary-element method for Stokes flow in partially covered cavities / C. H. Driesen, J. G. M. Kuerten // Computational Mechanics. - 2000.

- V. 25. - P. 501-513.

11. Петров, Н.П. Гидродинамическая теория смазки. Избранные работы / Н.П. Петров. - АН СССР, 1948. -556с.

12. Moffat, H.K. Viscous and resistive eddies near a sharp corner/ H.K. Moffat// J. Fluid Mech. - 1964. - V.18. - P. 1-18.

13. van der Woude, D. Stokes flow in a rectangular cavity by rotlet forcing / D. van der Woude, H. J. H. Clercx, G. J. F. van Heijst, V. V. Meleshko // Physics of Fluids. -2007. - V. 19. - P. 083602-19.

14. Gürcan, F. Eddy genesis and transformation of Stokes flow in a double-lid driven cavity / F.Gürcan, P.H. Gaskell, M.D. Savage, M.C.T.Wilson // J. Mec. Eng. Sci. -2003. - V. 217. - P. 353-363.

15. Krasnopolskaya, T.S. Steady Stokes flow in an annular cavity / T.S. Krasnopolskaya, V.V. Meleshko, G.W.M. Peters, H.E.H. Meijer // Q. J. Mech. Appl. Math. - 1996. - V. 49. - P. 593-619.

16. Saatdjian, E. On the solution of Stokes equations between confocal ellipses /

E. Saatdjian, N. Midoux, J. C. Andre // Phys. Fluids. - 1994. - V. 6. - P. 3833-3846.

17. Ramkissoon, H. Stokes flow past a slightly deformed fluid sphere / H. Ramkissoon // Zeitschrift für angewandte Mathematik und Physik ZAMP. - 1986. - V. 37.

- P. 859-866.

18. Happel, J. Low Reynolds Number Hydrodynamics / J.Happel, H.Brenner. -Prentice-Hall, 1965. -553 p.

19. Hackborn, W.W. An analysis of a Stokes flow in an annular region / W.W. Hackborn // Canadian Applied Mathematics Quarterly. - 2000. - V. 8. - № 2.

- P. 171-183.

20. Murad, A. Stokes Flow between Two Cylinders / A.Murad, S.K. Sen // Journal of Bangladesh Academy of Sciences. - 2012. - V. 36. - № 1. - P. 123-135.

21. Belonenko, M.B. Soliton-induced flow in carbon nanotube / M.B. Belonenko, S.A.Chivilikhin, V.V.Gusarov, I.Yu. Popov, O.A.Rodygina // Europhys. Lett. - 2013.

- V. 101. - P. 66001-3.

22. Zhong, S. Analytic solutions for Stokes flow with lateral variations in viscosity / S. Zhong // Geophysical Journal International. - 1996. - V. 124. - P. 18-28.

23. Wang, C.Y. Flow over a surface with parallel grooves / C.Y. Wang // Phys. Fluids.

- 2003. - V. 15. - P. 1114-1121.

24. Driesen, C.H. Low-Reynolds number flow over partially covered cavities / C.H. Driesen, J.G.M. Kuerten, M. Streng // Journal of Engineering Mathematics. -

1998. - V. 34. - P. 3-21.

25. Chen, G. A fast finite difference method for biharmonic equations on irregular domains and its application to an incompressible Stokes flow / G. Chen, Z. Li, P. Lin // Advances in Comput. Math. - 2008. - V. 29. - P. 113-133.

26. Gerya, T.V. An adaptive staggered grid finite difference method for modeling geodynamic Stokes flows with strongly variable viscosity / T.V. Gerya, D.A. May, T. Duretz // Geochemistry, Geophysics, Geosystems. - 2013. - V. 14. - P. 1200-1225.

27. Gerya, T.V. Characteristic-based marker-in-cell method with conservative finite-difference schemes for modeling geological flows with strongly variable transport properties / T.V. Gerya, D.A. Yuen // Phys. Earth Planet. In. - 2003. - V. 140.

- P. 293-318.

28. Arnold, D.N. A stable finite element for the Stokes equations / D. N. Arnold, F. Brezzi, M. Fortin // Calcolo. - 1984. - V. 21. - P. 337-344.

29. Dohrmann, С. A stabilized finite element method for the Stokes problem based on polynomial pressure projections / C. Dohrmann, P. Bochev // International Journal of Numerical Methods in Fluids. - 2004. - V. 46. - P. 183-201.

30. Harder, H. A finite-volume solution method for thermal convection and dynamo problems in spherical shells / H. Harder, U. Hansen. // Geophysical Journal International. - 2005. - V. 161. - P. 522-532.

31. Ye, X. A discontinuous finite volume method for the Stokes problems / X. Ye // SIAM J. Numerical Analysis. - 2006. - V. 44. - P. 183-198.

32. Wang, J. A new finite volume method for the Stokes problems./ J.Wang; Y.Wang; X.Ye // International Journal of Numerical Analysis and Modeling. - 2010. - V. 7.

- № 2. - P. 281-302.

33. Wesseling, P. An Introduction to Multigrid Methods / P. Wesseling. -R.T. Edwards, Inc, 2004. -296 p.

34. Федоренко, Р.П. Релаксационный метод решения разностных эллиптических уравнений / Р.П. Федоренко // Ж. вычисл. матем. и матем. Физ. - 1961. - T.1.

- № 5. -С. 922-927.

35. Федоренко, Р.П. О скорости сходимости одного итерационного процесса

/ P.n. OegopeHKO // bmhhch. MaTeM. h MaTeM. oh3. - 1964. - T.4. - № 3.

- C. 559-564.

36. Weber, H. A multigrid method for a Petrov-Galerkin discretization of the Stokes equations / H. Weber // Numerische Mathematik. - 1993. - V. 66. - P. 525-541.

37. Wang, M. Multigrid methods for the Stokes equations using distributive Gauss-Seidel relaxations based on the least squares commutator / M. Wang, L. Chen // Journal of Scientific Computing. - 2013. - V. 56. - P. 409-431.

38. Brenner, S.C. A nonconforming multigrid method for the stationary stokes problem / S.C. Brenner // Math. Comp. - 1990. - V. 55. - P. 411-437.

39. John, V. A coupled multigrid method for nonconforming finite element discretizations of the 2D-Stokes equation / V. John , L. Tobiska // Computing. - 2000.

- V. 64. - P. 307-321.

40. Altas, I. Multigrid Solution of Automatically Generated High-Order discretizations for the biharmonic Equation / I.Altas, J.Dym, M. Gupta, P. Manohar // SIAM J. Scientific Computing. - 1998. - V. 19. - P. 1575-1585.

41. Notay, Y. A new algebraic multigrid approach for Stokes problems / Y. Notay // Numerische Mathematik. - 2016. - V. 132. - P. 51-84.

42. Tackley, P.J. Effects of strongly temperature-dependent viscosity on time-dependent, three-dimensional models of mantle convection / P.J. Tackley // Geophys.Res.Lett. - 1993. - V. 20. - P. 2187-2190.

43. Kageyama, A. Geodynamo and mantle convection simulations on the Earth Simulator using the Yin-Yang grid / A.Kageyama, M.Yoshida // Journal of Physics: Conference Series. - 2005. - V. 16. - P. 325-338.

44. Moresi, L. Numerical investigation of 2D convection with extremely large viscosity variations / L. Moresi, V. Solomatov // Physics of Fluids. - 1995. - V. 7.

- P. 2154-2162.

45. Trompert, R.A. The application of a finite volume multigrid method to three-dimensional flow problems in a highly viscous fluid with a variable viscosity / R.A.Trompert, U.Hansen // Geophys.Astrophys.Fluid.Dyn. - 1996. - V. 83.

- P. 261-291.

46. Furuichi, M. Development of a Stokes flow solver robust to large viscosity jumps using a Schur complement approach with mixed precision arithmetic / M. Furuichi, D.A.May, P.J.Tackley // Journal of Computational Physics. - 2011. - V. 230.

- P. 8835-8851.

47. Popov, A.I. On the Stokes flow computation algorithm based on Woodbury formula / A.I. Popov, I.S. Lobanov, I.Yu. Popov, T.V. Gerya // Nanosystems: Physics, Chemistry, Mathematics. - 2015. - V.6. - № 1. - P. 140-145.

48. Lobanov, I.S. Numerical approach to the Stokes problem with high contrasts in viscosity / I.S. Lobanov, I.Yu. Popov , A.I. Popov, T.V. Gerya // Applied Mathematics and Computation. - 2014. - V. 235. - P. 17-25.

49. Auth, C. Multigrid solution of convection problems with strongly variable viscosity / C. Auth, H. Harder // Geophysical Journal International. - 1999. - V. 137.

- P. 793-804.

50. Kameyama, M. Multigrid iterative algorithm using pseudo-compressibility for three-dimensional mantle convection with strongly variable viscosity / M. Kameyama,

A. Kageyama, T. Sato // Journal of Computational Physics. - 2005. - V. 206.

- P. 162-181.

51. Stemmer, K. A new method to simulate convection with strongly temperature- and pressure-dependent viscosity in a spherical shell: Applications to the Earth's mantle / K. Stemmer, H. Harder, U. Hansen // Physics of the Earth and Planetary Interiors. -2006. - V. 157. - P. 223-249.

52. Deubelbeiss, Y. Comparison of Eulerian and Lagrangian numerical techniques for the Stokes equations in the presence of strongly varying viscosity / Y. Deubelbeiss,

B.J.P. Kaus // Phys. Earth Planet. In. - 2008. - V. 171. - P. 92-111.

53. Moresi, L. The accuracy of finite element solutions of Stokes's flow with strongly varying viscosity / L. Moresi, S. Zhong, M. Gurnis // Physics of the Earth and Planetary Interiors. - 1996. - V. 97. - P. 83-94.

54. Albers, M. A local mesh refinement multigrid method for 3D convection problems with strongly variable viscosity / M. Albers // Journal of Compuational Physics. - 2000.

- V. 160. - P. 126-150.

55. Edited by Li D. (Ed.) Encyclopedia of microfuidics and nanofuidics. - New York: Springer. - 2008.

56. Chivilikhin S. A. Model of fluid flow in a nano-channel / S.A. Chivilikhin, I.Yu. Popov, V.V. Gusarov, A.I. Svitenkov // Russian J. Math. Phys. - 2008. - V.15.

- P.409-411.

57. Li, T. Structured and viscous water in subnanometer gaps / T. Li, J. Gao, R. Szoszkeiwicz, U. Landnan, E. Riedo // Phys. Rev. B. - 2007. - V. 75. - P. 115415-6.

58. Iwase, Y. An interpretation of the Nusselt-Rayleigh number relationship for convection in a spherical shell / Y. Iwase, S. Honda // Geophysical Journal In ternational. - 1997. - V. 130. - P. 801-804.

59. Zebib, A. Infinite Prandtl number thermal convection in a spherical shell / A. Zebib, G. Schubert, J. M. Straus // Journal of Fluid Mechanics. - 1980. - V. 97. - P. 257-277.

60. Tackley, P. Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid / P. Tackley // Physics of the Earth and Planetary Interiors. - 2008. - V. 171. - P. 7-18.

61. Kameyama, M. Multigrid-based simulation code for mantle convection in spherical shell using yin-yang grid / M. Kameyama, A. Kageyama, T. Sato // Physics of the Earth and Planetary Interiors. - 2008. - V. 171. - P. 19-32.

62. Baumgardner, J.R. Icosahedral Discretization of the Two-Sphere / J.R. Baumgardner, P.O. Frederickson // SIAM Journal on Numerical Analysis. - 1985.

- V. 22. - P. 1107-1115.

63. Baumgardner, J.R. Three dimensional treatment of convection flow in the Earth's mantle / J.R. Baumgardner // Journal of Statistical Physics. - 1985. - V. 39.

- P. 501-511.

64. Choblet, G. Oedipus: a new tool to study the dynamics of planetary interiors / G. Choblet, O. Cadek, F. Couturier, C. Dumoulin // Geophys. J. Int. - 2007. - V. 170.

- P. 9-30.

65. Choblet, G. Modelling thermal convection with large viscosity gradients in one block of the 'cubed sphere' / G. Choblet // Journal of Computational Physics. - 2005.

- V. 205. - P. 269-291.

66. Ronchi, C. The"cubed sphere": a new method for the solution of partial differential equations in spherical geometry / C. Ronchi, R. Iacono, P.S. Paolucci // Journal of Computational Physics. - 1996. - V. 124. - P. 93-114.

67. Tabata, M. Finite element approximation to infinite Prandtl number Boussinesq equations with temperature-dependent coefficients - Thermal convection problems in a spherical shell / M. Tabata // Future Generation Computer Systems. - 2006. - V. 22.

- P. 521-531.

68. Tabata, M. A stabilized finite element method for the Rayleigh-Benard equations with infinite Prandtl number in a spherical shell / M. Tabata, A. Suzuki // Computer Methods in Applied Mechanics and Engineering. - 2000. - V. 190. - P. 387-402.

69. Zhong, S. Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection / S. Zhong, M.T. Zuber, L. Moresi , M. Gurnis // Journal of Geophysical Research. - 2000. - V. 105. - P. 11063-11082.

70. Glatzmaier, G.A. Numerical simulations of mantle convection: Time-dependent, three-dimensional, compressible, spherical shell / G.A. Glatzmaier // Geophys. Astrophys. Fluid Dyn. - 1988. - V. 43. - P. 223-264.

71. Young, R.E. Finite-amplitude thermal convection in a spherical shell / R.E. Young // Journal of Fluid Mechanics. - 1974. - V. 63. - P. 695-721.

72. Zhang, S. Some effects of lateral viscosity variations on geoid and surface velocities induced by density anomalies in the mantle / S. Zhang, U. Christensen // Geophys. J. Int. - 1993. - V. 114. - P. 531-547.

73. Evonuk, M. A 2D study of the effects of the size of a solid core on the equatorial flow in giant planets / M. Evonuk , G.A. Glatzmaier // Icarus. - 2006. - V. 181.

- P. 458-464.

74. Gerya, T.V. Robust characteristics method for modeling multiphase visco-elasto-plastic thermo-mechanical problems / T.V. Gerya, D.A. Yuen // Physics of the Earth and Planetary Interiors. - 2007. - V. 163. - P. 83-105.

75. Gumis, M. Generation of long wavelength heterogeneity in the mantle by the dynamic interaction between plates and convection / M. Gumis, S. Zhong // Geophys. Res. Lett. - 1991. - V. 18. - P. 581-584.

76. Thieulot, C. ELEFANT: a user-friendly multipurpose geodynamics code / C. Thieulot // Solid Earth Discuss. - 2014. - V. 6. - P. 1949-2096.

77. Rodygina, O. A. Crystallite model for flow in nanotube caused by wall soliton / O. A. Rodygina, S. A. Chivilikhin, I.Yu. Popov, V. V. Gusarov // Nanosystems: Phys., Chem., Math. - 2014. - V. 5. - № 3. - P. 400-404.

78. Popov, I.Yu. Practical analytical solutions for benchmarking of 2-D and 3-D geodynamic Stokes problems with variable viscosity / I. Yu. Popov, I. S. Lobanov, S. I. Popov, A. I. Popov, T. V. Gerya // Solid Earth. - 2014. - V. 5. - P. 461-476.

79. Turcotte, D.L. Geodynamics / D.L. Turcotte, G. Schubert. - Third edition. -Cambridge University Press, 2014. - 636 p.

80. Kramer, S.C. An implicit free surface algorithm for geodynamical simulations / S.C. Kramer, C.R. Wilson, D.R. Davies // Physics of the Earth and Planetary Interiors. -2012. - V. 194. - P. 25-37.

81. Higdon, J.J.L. Stokes flow in arbitrary two-dimensional domains: shear flow over ridges and cavities / J.J.L. Higdon // Journal of Fluid Mechanics. - 1985. - V. 159.

- P. 195-226.

82. Arad, M. A highly accurate numerical solution of a biharmonic equation / M. Arad, A. Yakhot, G. Ben-Dor // Numerical Methods for Partial Differential Equations. - 1997.

- V. 13. - P. 375-391.

83. Tosi, N. Semi-analytical solution for viscous Stokes flow in two eccentrically nested spheres / N. Tosi, Z. Martinec // Geophys. J. Int. - 2007. - V. 170. - P. 1015-1030.

84. van Keken, P. A community benchmark for subduction zone modeling / P. van Keken, C. Currie, S.D. King, M.D. Behn, C. Amandine, J.He, R.F. Katz, S.Lin,

E.M. Parmentier, M. Spiegelman, K. Wang // Physics of the Earth and Planetary Interiors. - 2008. - V. 171. - P. 187-197.

85. Schmeling, H. A benchmark comparison of spontaneous subduction models-Towards a free surface/ H. Schmeling, A.Y. Babeyko, A. Enns, C. Faccenna,

F. Funiciello, T. Gerya, G.J. Golabek, S. Grigull, B.J.P. Kaus, G. Morra, S.M. Schmalholz, J. van Hunen // Physics of the Earth and Planetary Interiors. - 2008.

- V. 171. - P. 198-223.

86. Travis, B. J. A benchmark comparison of numerical methods for infinite Prandtl number thermal convection in two-dimensional Cartesian geometry / B. J. Travis, C. Anderson, J. Baumgardner, C. W. Gable, B. H. Hager, R. J. O'Connell, P. Olson, A. Raefsky, G. Schubert // Geophysical and Astrophysical Fluid Dynamics. - 1990.

- V. 55. - P. 137-160.

87. Busse, F.H. 3D convection at infinite Prandtl number in Cartesian geometry - a benchmark comparison / F. H. Busse, U. Christensen, R. Clever, L. Cserepes, C. Gable, E. Giannandrea, L. Guillou, G. Houseman, H. C. Nataf, M. Ogawa, M. Parmentier, C. Sotin, B. Travis // Geophysical and Astrophysical Fluid Dynamics. - 1994. - V. 75.

- P. 39-59.

88. Makeev, I.V. Analytical benchmark solutions for nanotube flows with variable viscosity / I.V. Makeev, I.V. Blinova, I.Y. Popov // Nanosystems: Physics, Chemistry, Mathematics - 2015. - V. 6. - № 5. - P. 672-679.

ПРИЛОЖЕНИЕ 1. Реализация операции smoothing в программе Stokes3D-cylinder.

% return new approximation for velocity and pressure (vx,vy,vz,pr) % and distribution of residuals (resx,resy,resz,resc)

% (RX,RY,RZ,RC) is a distribution of right parts for equations % (etaxy,etaxz,etayz,etan) is a viscosity distribution

% notation (x,y,z) corresponds to (r,phi,z) in cylindrical % coordinates

function[vx,resx,vy,resy,vz,resz,pr,resc]=

Stokes_cylindrical_smoother (prnorm,etaxy,etaxz,etayz,etan, iternum,krelaxs,krelaxc,xnum,ynum,znum,xstp,ystp,zstp,RX,vx,RY,vy, RZ,vz,RC,pr)

xkf1=1/xstp; ykf1=1/ystp; zkf1=1/zstp;

xkf=1/xstpA2;

ykf=1/ystpA2;

zkf=1/zstpA2;

xkf2=2/xstpA2;

ykf2=2/ystpA2;

zkf2=2/zstpA2;

xykf=1/xstp/ystp;

xzkf=1/xstp/zstp;

yzkf=1/ystp/zstp;

%offset from zero for r component shift=1;

% Gauss-Seidel vx, xy, xz, P iteration circle for niter=1:1:iternum; for i=1:1:ynum+1;

for j=1:1:xnum+1;

for k=1:1:znum+1; if (j<xnum+1)

if (i==1 || i==ynum+1 || j==1 ||

j==xnum || k==1 || k==znum+1 );

vx(i,j,k)=RX(i,j,k);

if(i==1)

vx(i,j,k)=2*RX(i,j,k) -vx(i + 1,j,k);

end

if(i==ynum+1)

vx(i,j,k)=2*RX(i,j,k) -vx(i-1,j,k);

if(k==1)

vx(i,j,k)=2*RX(i,j,k) -vx(i,j,k+1);

end

if(k==znum+1)

vx(i,j,k)=2*RX(i,j,k) -vx(i,j,k-1);

end

else

r_j=(j-1)*xstp+shift;

resxcur=RX(i,j,k)+(pr(i-1,j,k-1)-pr(i-1,j-1,k-1))/xstp; %dtau_r_r/dr

resxcur=resxcur-(xkf2*(etan(i-1,j,k-1)*(vx(i,j + 1,k)-vx(i,j,k))-etan(i-1,j-1,k-1)*(vx(i,j,k) -vx(i,j-1,k))));

%dtau_r_phi/dphi

tau_r_phi=ykf1*(etaxy(i,j,k-1)*((1/r_j)*ykf1*(vx(i+1,j,k)-vx(i,j,k)) +xkf1*(vy(i,j + 1,k) -vy(i,j,k))-(1/r_j)*(vy(i,j + 1,k)+vy(i,j,k))*0.5));

tau_r_phi=tau_r_phi-ykf1*(etaxy(i-1,j,k-1)*((1/r_j)*ykf1*(vx(i,j,k)-vx(i-1,j,k))+xkf1*(vy(i-1,j + 1,k)-vy(i-1,j,k))-(1/r_j)* (vy(i-1,j+1,k)+vy(i-1,j,k))*0.5));

resxcur=resxcur-(1/r_j)*tau_r_phi;

% dtau_r_z/dz

resxcur=resxcur-zkf1*(etaxz(i-1,j,k)*(xkf1*(vz(i,j+1,k)-vz(i,j,k))+ zkf1*(vx(i,j,k+1) -vx(i,j,k))) -etaxz(i-1,j,k-1)*(xkf1*(vz(i,j+1,k-1)-vz(i,j,k-1))+zkf1*(vx(i,j,k)-vx(i,j,k-1))));

%tau_r_r

tau_r_r=xkf1*(etan(i-1,j,k-1)+etan(i-1,j-1,k-1))* (vx(i,j + 1,k) -vx(i,j -1,k))*0.5;

%tau_phi_phi

tau_phi_phi=(1/r_j)*(etan(i-1,j,k-1)+

etan(i-1,j-1,k-1))*(ykf1*((vy(i,j+1,k)+vy(i,j,k))*0.5-(vy(i-1,j+1,k)+vy(i-1,j,k))*0.5) + vx(i,j,k));

resxcur=resxcur-(1/r_j)*(tau_r_r-tau_phi_phi);

% Koefficient at vx(i,j,k)

kfxcur=-xkf2*(etan(i-1,j,k-1)+etan(i-1,j-1,k-1))-

(1/(r_jA2))*(etan(i-1,j,k-1)+etan(i-1,j-1,k-1))-

(1/(r_jA2))*ykf*(etaxy(i,j,k-1)+etaxy(i-1,j,k-1))-

zkf*(etaxz(i-1,j,k)+etaxz(i-1,j,k-1));

% Updating solution

vx(i,j,k)=vx(i,j,k)+resxcur/kfxcur*krelaxs;

end

if (i<ynum+1)

if (i==1 Il i==ynum Il j==1 Il j==xnum+1 Il k==1 II k==znum+1);

vy(i,j,k)=RY(i,j,k);

if (j==1)

vy(i,j,k)=2*RY(i,j,k) -vy(i,j + 1,k);

end

if(j==xnum+1)

vy(i,j,k)=2*RY(i,j,k)-vy(i,j-1,k);

end if(k==1)

vy(i,j,k)=2*RY(i,j,k)-vy(i,j,k+1);

end

if(k==znum+1)

vy(i,j,k)=2*RY(i,j,k)-vy(i,j,k-1);

end else

r_j=(j-1)*xstp+shift; r_j1=(j-2)*xstp+shift; r_j12=r_j-xstp*Q.5;

resycur=RY(i,j,k) + (1/r_j12)*(pr(i,j-1,k-1)-pr(i-1,j-1,k-1))/ystp; % dtau_r_phi/dr

dtau_r_phi=xkf1*(etaxy(i,j,k-1)*((1/r_j)*ykf1*(vx(i+1,j,k)-vx(i,j,k))+xkf1*(vy(i,j + 1,k) -vy(i,j,k))-(1/r_j)*(vy(i,j + 1,k)+vy(i,j,k))*Q.5));

dtau_r_phi=dtau_r_phi-xkf1*(etaxy(i,j-1,k-1)*((1/r_j1)* ykf1*(vx(i + 1,j-1,k) -vx(i,j-1,k))+xkf1*(vy(i,j,k) -vy(i,j-1,k))-(1/r_j1 ) *(vy(i,j,k)+vy(i,j-1,k))*Q.5));

resycur=resycur-dtau_r_phi;

% tau_r_phi

tau_r_phi=Q.5*(etaxy(i,j,k-1)+etaxy(i,j-1,k-1))* ((1/r_j12)*ykf1*(Q.5*(vx(i + 1,j,k)+vx(i + 1,j-1,k)) -Q.5*(vx(i,j,k)+vx(i,j-1,k))) +

xkf1*Q.5*(vy(i,j+1,k) - vy(i,j-1,k))-(1/r_j12)*vy(i,j,k)); resycur=resycur-2*(1/r_j12)*tau_r_phi; % dtau_phi_phi/dphi

dtau_phi_phi=2*ykf1*(etan(i,j-1,k-1)*((1/(r_j12))* (ykf1*(vy(i+1,j,k)-vy(i,j,k))) +

(1/(r_j12))*(Q.5*(vx(i + 1,j,k)+vx(i + 1,j-1,k)))));

dtau_phi_phi=dtau_phi_phi-2*ykf1*(etan(i-1,j-1,k-1)* ((1/(r_j12))*(ykf1*(vy(i,j,k)-vy(i-1,j,k))) + (1/(r_j12))*(Q.5*(vx(i,j,k)+vx(i,j-1,k)))));

resycur=resycur-(1/r_j12)*dtau_phi_phi; % dtau_z_phi/dz

dtau_z_phi=zkf1*(etayz(i,j-1,k)*(((zkf1*(vy(i,j,k+1)-vy(i,j,k))) + (1/(r_j12))*(ykf1*(vz(i + 1,j,k)-vz(i,j,k))))));

dtau_z_phi=dtau_z_phi-zkf1*(etayz(i,j-1,k-1)*(((zkf1*(vy(i,j,k)-vy(i,j,k-1))) + (1/(r_j12))*(ykf1*(vz(i + 1,j,k-1)-vz(i,j,k-1))))));

resycur=resycur-dtau_z_phi;

kphicur=-xkf*(etaxy(i,j,k-1)+etaxy(i,j-1,k-1))-(1/(r_j12A2))*ykf2*(etan(i,j-1,k-1)+etan(i-1,j-1,k-1)) -zkf*(etayz(i,j-1,k)+etayz(i,j-1,k-1));

kphicur=kphicur-xkf1*0.5*((1/r_j)* etaxy(i,j,k-1) -

(1/r_j1)*etaxy(i,j-1,k-1))-(1/(r_j12A2))*

(etaxy(i,j,k-1)+etaxy(i,j-1,k-1));

vy(i,j,k)=vy(i,j,k)+resycur/kphicur*krelaxs;

end

end

if (k<znum+1)

if (i==1 Il i==ynum+1 Il j==1 II j==xnum+1 jj k==1 jj k==znum); vz(i,j,k)=RZ( i,j,k);

if(i==1)

vz(i,j,k)=2*RZ(i,j,k) -vz(i + 1,j,k);

end

if(i==ynum+1)

vz(i,j,k)=2*RZ(i,j,k) -vz(i-1,j,k);

end if (j==1)

vz(i,j,k)=2*RZ(i,j,k) -vz(i,j + 1,k);

end

if(j==xnum+1)

vz(i,j,k)=2*RZ(i,j,k) -vz(i,j-1,k);

end else

r_j=(j-1)*xstp+shift; r_j12=r_j-xstp*0.5;

reszcur=RZ(i,j,k)+(pr(i-1,j-1,k)-pr(i-1,j-1,k-1))/zstp; % dtau_r_z/dr

dtau_r_z=(etaxz(i-1,j,k)*(xkf*(vz(i,j + 1,k)-vz(i,j,k))+xzkf*(vx(i,j,k+1) -vx(i,j,k)))-etaxz(i-1,j-1,k)*(xkf*(vz(i,j,k) -vz(i,j-1,k)) + xzkf*(vx(i,j-1,k+1) -vx(i,j-1,k))));

reszcur=reszcur-dtau_r_z; % tau_r_z

tau_r_z=Q.5*(etaxz(i-1,j,k)+etaxz(i-1,j-1,k))* (Q.5*xkf1*(vz(i,j + 1,k) - vz(i,j-1,k))+zkf1*(Q.5*(vx(i,j,k+1) + vx(i,j-1,k+1))-Q.5*(vx(i,j,k)+vx(i,j-1,k))));

reszcur=reszcur-(1/r_j12)*tau_r_z;

% dtau_phi_z/dphi

dtau_phi_z=ykf1*(etayz(i,j-1,k)*((zkf1*(vy(i,j,k+1)-vy(i,j,k))) + (1/(r_j12))*(ykf1*(vz(i + 1,j,k)-vz(i,j,k)))));

dtau_phi_z=dtau_phi_z-ykf1*(etayz(i-1,j-1,k)*((zkf1*(vy(i-1,j,k+1)-vy(i-1,j,k))) + (1/(r_j12))*(ykf1*(vz(i,j,k)-vz(i-1,j,k)))));

reszcur=reszcur-(1/r_j12)*dtau_phi_z;

% dtau_z_z/dz

reszcur=reszcur-(zkf2*(etan(i-1,j-1,k)*(vz(i,j,k+1)-vz(i,j,k))-etan(i-1,j-1,k-1)*(vz(i,j,k)-vz(i,j,k-1))));

kfzcur=-xkf*(etaxz(i-1,j,k)+etaxz(i-1,j-1,k))-

(1/(r_j12A2))*ykf*(etayz(i,j-1,k)+etayz(i-1,j-1,k))-

zkf2*(etan(i-1,j-1,k)+etan(i-1,j-1,k-1));

vz(i,j,k)=vz(i,j,k)+reszcur/kfzcur*krelaxs;

end

end

if (i<ynum && j<xnum && k<znum);

r_j=(j-1)*xstp+shift; r_j2=r_j+xstp*Q.5;

resc(i,j,k)=RC(i,j,k)- (1/r_j2)*(vx(i+1,j+1,k+1)+vx(i+1,j,k+1))*Q.5-(vx(i+1,j+1,k+1)-vx(i+1,j,k+1))/xstp-(1/r_j2)*(vy(i+1,j+1,k+1)-vy(i,j+1,k+1))/ystp-(vz(i+1,j+1,k+1)-vz(i+1,j+1,k))/zstp;

pr(i,j,k)=pr(i,j,k)+resc(i,j,k)*etan(i,j,k)*krelaxc;

end

end

end

end

end

end

ПРИЛОЖЕНИЕ 2. Реализация операции smoothing в программе Stokes3D-spher.

% return new approximation for velocity and pressure (vx,vy,vz,pr) % and distribution of residuals (resx,resy,resz,resc)

% (RX,RY,RZ,RC) is a distribution of right parts for equations % (etaxy,etaxz,etayz,etan) is a viscosity distribution

% notation (x,y,z) corresponds to (r,teta,phi) in spherical % coordinates

function[vx,resx,vy,resy,vz,resz,pr,resc]=

Stokes_ spherical_smoother(prnorm,etaxy,etaxz,etayz,etan,iternum, krelaxs,krelaxc,xnum,ynum,znum,xstp,ystp,zstp,RX,vx,RY,vy,RZ,vz,RC,p

r)

xkf1=1/xstp; ykf1=1/ystp; zkf1=1/zstp;

%offset from zero for r component

shift_r=1;

shift_teta=0.5;

for niter=1:1:iternum; for i=1:1:ynum+1;

for j=1:1:xnum+1;

for k=1:1:znum+1;

if (j<xnum+1)

if (i==1 || i==ynum+1 || j==1 || j==xnum || k==1 || k==znum+1 );

vx(i,j,k)=RX(i,j,k);

if(i==1)

vx(i,j,k)=2*RX(i,j,k) -vx(i + 1,j,k);

end

if(i==ynum+1)

vx(i,j,k)=2*RX(i,j,k) -vx(i-1,j,k);

end if(k==1)

vx(i,j,k)=2*RX(i,j,k) -vx(i,j,k+1);

end

if(k==znum+1)

vx(i,j,k)=2*RX(i,j,k) -vx(i,j,k-1);

else

r_j=(j-1)*xstp+shift_r;

teta_i=(i-1)*ystp+shift_teta;

teta_i1=teta_i-Q.5*ystp;

resxcur=RX(i,j,k)+(pr(i-1,j,k-1)-pr(i-1,j-1,k-1))/xstp; %tau_r_r

tau_r_r = Q.5*(etan(i-1,j,k-1)+

etan(i-1,j-1,k-1))*xkf1*(vx(i,j + 1,k) -vx(i,j-1,k));

resxcur=resxcur-2*(1/r_j)*tau_r_r;

%dtau_r_r/dr

resxcur=resxcur-(2*xkf1*(etan(i-1,j,k-1)*xkf1*(vx(i,j+1,k)-vx(i,j,k))-etan(i-1,j-1,k-1)*xkf1*(vx(i,j,k)-vx(i,j-1,k))));

%tau_r_teta

tau_r_teta = Q.5*(etan(i-1,j,k-1)+etan(i-1,j-1,k-1))* (xkf1*(Q.5*(vy(i,j + 1,k)+vy(i-1,j + 1,k)) -Q.5*(vy(i,j,k)+vy(i-1,j,k))));

tau_r_teta = tau_r_teta + Q.5*(etan(i-1,j,k-1) + etan(i-1,j-1,k-1))* (-(1/r_j)*(1/4)*(vy(i,j + 1,k)+vy(i-1,j + 1,k)+vy(i,j,k) + vy(i-1,j,k)) + (1/r_j)*Q.5*ykf1*(vx(i + 1,j,k)-vx(i-1,j,k)));

resxcur=resxcur-(1/r_j)*cot(teta_i1)*tau_r_teta;

%dtau_r_teta/dteta

dtau_r_teta = ykf1*(etaxy(i,j,k-1)*(xkf1*(vy(i,j+1,k)-vy(i,j,k))-(1/r_j)*Q.5*(vy(i,j+1,k)+vy(i,j,k))+ ykf1*(1/r_j)* (vx(i + 1,j,k) -vx(i,j,k))));

dtau_r_teta =dtau_r_teta - ykf1*(etaxy(i-1,j,k-1)*(xkf1* (vy(i-1,j + 1,k)-vy(i-1,j,k)) - Q.5*(1/r_j)*(vy (i-1,j + 1,k) + vy(i-1,j,k))+ ykf1*(1/r_j)*(vx(i,j,k)-vx(i-1,j,k))));

resxcur=resxcur-(1/r_j)*dtau_r_teta;

%dtau_r_phi/dphi

dtau_r_phi=zkf1*(etaxz(i-1,j,k)*((1/(r_j*sin(teta_i1)))* zkf1*(vx(i,j,k+1) -vx(i,j,k) )+xkf1*(vz(i,j + 1,k) -vz(i,j,k))-(1/r_j)*Q.5*(vz(i,j + 1,k)+vz(i,j,k)))); dtau_r_phi=dtau_r_phi-zkf1*(etaxz(i-1,j,k-1)* ((1/(r_j*sin(teta_i1)))*zkf1*(vx(i,j,k)-vx(i,j,k-1))+ xkf1*(vz(i,j + 1,k-1)-vz(i,j,k-1))-(1/r_j)* Q.5*(vz(i,j + 1,k-1)+vz(i,j,k-1))));

resxcur=resxcur-(1/(r_j*sin(teta_i1)))*dtau_r_phi; %tau_teta_teta

tau_teta_teta = (etan(i-1,j,k-1) + etan(i-1,j-1,k-1))* ((1/r_j)*ykf1*(Q.5*(vy(i,j + 1,k)+vy(i,j,k))-Q.5*(vy(i-1,j + 1,k)+vy(i-1,j,k))) + (1/r_j)*vx(i,j,k));

resxcur=resxcur+(1/r_j)*tau_teta_teta; %tau_phi_phi

tau_phi_phi=(etan(i-1,j,k-1)+etan(i-1,j-1,k-1))* ((1/(r_j*sin(teta_i1)))*zkf1*(0.5*(vz(i,j+1,k)+vz(i,j,k))-0.5*(vz(i,j + 1,k-1)+vz(i,j,k-1))) + (1/r_j)*vx(i,j,k) + (1/r_j)*(1/4)*(vy(i,j + 1,k) + vy(i-1,j+1,k)+vy(i,j,k)+vy(i-1,j,k))*cot(teta_i1));

resxcur=resxcur+(1/r_j)*tau_phi_phi;

kfxcur=-2*(xkf1A2)*(etan(i-1,j,k-1)+etan(i-1,j-1,k-1))-(1/(r_j)A2)*(ykf1A2)*(etaxy(i,j,k-1)+etaxy(i-1,j,k-1));

kfxcur=kfxcur-(1/(r_j)A2)*(zkf1A2)*(1/(sin(teta_i1))A2)*

(etaxz(i-1,j,k)+etaxz(i-1,j,k-1))-2*(1/(r_j)A2)*

(etan(i-1,j,k-1)+etan(i-1,j-1,k-1));

vx(i,j,k)=vx(i,j,k)+resxcur/kfxcur*krelaxs;

end

end

if (i<ynum+1)

if (i==1 Il i==ynum Il j==1 II j==xnum+1 jj k==1 jj k==znum+1); vy(i,j,k)=RY(i,j,k);

if (j==1)

vy(i,j,k)=2*RY(i,j,k) -vy(i,j + 1,k);

end

if(j==xnum+1)

vy(i,j,k)=2*RY(i,j,k)-vy(i,j-1,k);

end if(k==1)

vy(i,j,k)=2*RY(i,j,k)-vy(i,j,k+1);

end

if(k==znum+1)

vy(i,j,k)=2*RY(i,j,k)-vy(i,j,k-1);

end else

r_j=(j-1)*xstp+shift_r; r_j1=(j-2)*xstp+shift_r; r_j12=r_j-xstp*0.5; teta_i=(i-1)*ystp+shift_teta;

resycur=RY(i,j,k) + (1/r_j12)*(pr(i,j-1,k-1)-pr(i-1,j-1,k-1))/ystp; %tau_r_teta

tau_r_teta = 0.5*(etan(i,j-1,k-1)+

etan(i-1,j-1,k-1))*(0.5*xkf1*(vy(i,j + 1,k)-vy(i,j-1,k))-0.5*(1/r_j12)*(vy(i,j + 1,k)+vy(i,j-1,k)) + (1/r_j12)*ykf1*(0.5*(vx(i + 1,j,k)+vx(i + 1,j-1,k)) -0.5*(vx(i,j,k)+vx(i,j-1,k))));

resycur = resycur - 2*(1/r_j12)*tau_r_teta

%dtau_r_teta/dr

dtau_r_teta = xkf1*(etaxy(i,j,k-1)*(xkf1*(vy(i,j + 1,k)-vy(i,j,k))-(1/r_j)*0.5*(vy(i,j + 1,k)+vy(i,j,k)) + (1/r_j)*ykf1*(vx(i + 1,j,k)-vx(i,j,k))));

dtau_r_teta = dtau_r_teta - xkf1*(etaxy(i,j-1,k-1)*(xkf1*(vy(i,j,k)-vy(i,j-1,k))-(1/r_j1)*0.5*(vy(i,j,k)+vy(i,j-1,k)) + (1/r_j1)*ykf1*(vx(i + 1,j-1,k)-vx(i,j-1,k))));

resycur=resycur-dtau_r_teta;

%tau_teta_teta

tau_teta_teta = (etan(i,j-1,k-1)+etan(i-1,j-1,k-1))* ((1/r_j12)*ykf1*0.5*(vy(i+1,j,k)-vy(i-1,j,k))+

(1/r_j12)*(1/4)*(vx(i,j,k)+vx(i,j-1,k)+vx(i + 1,j,k)+vx(i + 1,j-1,k))); resycur=resycur-(1/r_j12)*cot(teta_i)*tau_teta_teta; %dtau_teta_teta/dteta

dtau_teta_teta = 2*ykf1*(etan(i,j-1,k-1)* ((1/r_j12)*ykf1*(vy(i+1,j,k) -vy(i,j,k)) + (1/r_j12)* 0.5*(vx(i + 1,j,k)+vx(i + 1,j-1,k))));

dtau_teta_teta = dtau_teta_teta-2*ykf1*(etan(i-1,j-1,k-1)* ((1/r_j12)*ykf1*(vy(i,j,k) -vy(i-1,j,k)) + (1/r_j12)*0.5*(vx(i,j,k)+vx(i,j-1,k))));

resycur=resycur-(1/r_j12)*dtau_teta_teta;

%dtau_teta_phi/dphi

dtau_teta_phi = zkf1*etayz(i,j-1,k)*((1/r_j12)*ykf1*(vz(i+1,j,k)-vz(i,j,k))-(1/r_j12)*cot(teta_i)*0.5*(vz(i + 1,j,k) +

vz(i,j,k))+(1/r_j12)*(1/sin(teta_i))*zkf1*(vy(i,j,k+1)-vy(i,j,k)));

dtau_teta_phi = dtau_teta_phi-zkf1*etayz(i,j-1,k-1)* ((1/r_j12)*ykf1*(vz(T+1,j7k-1)-vz(i,j,k-1))-(1/r_j12)*cot(teta_i)*0.5*(vz(i+1,j,k-1)+ vz(i,j,k-1))+(1/r_j12)*(1/sin(teta_i))* zkf1*(vy(i,j,k) -vy(i,j,k-1)));

resycur=resycur-(1/r_j12)*(1/sin(teta_i))*dtau_teta_phi;

resycur = resycur-(1/r_j12)*tau_r_teta;

%tau_phi_phi

tau_phi_phi = (etaxy(i,j,k-1)+etaxy(i,j-1,k-1))*

((1/r_j12)*(1/sin(teta_i))*zkf1*(Q.5*(vz(i+1,j,k)+vz(i,j,k))-Q.5*(vz(i + 1,j,k-1)+vz(i,j,k-1))) + (1/r_j12)*(1/4)*(vx(i,j,k) + vx(i,j-1,k)+vx(i + 1,j,k)+vx(i + 1,j-1,k) ) + (1/r_j12)* cot(teta_i)*vy(i,j,k));

resycur = resycur + (1/r_j12)*cot(teta_i)*tau_phi_phi;

% Koefficient at v_teta(i,j,k)

kphicur=-(xkf1A2)*(etaxy(i,j,k-1)+etaxy(i,j-1,k-1)) -xkf1*Q.5*((1/r_j)* etaxy(i,j,k-1) - (1/r_j1)* etaxy(i,j-1,k-1));

kphicur=kphicur - 2*(ykf1A2)*(1/(r_j12)A2)*(etan(i,j-1,k-1)+ etan(i-1,j-1,k-1))-(zkf1A2)*(1/(r_j12)A2)*(1/(sin(teta_i))A2)* (etayz(i,j-1,k)+etayz(i,j-1,k-1))-(1/(r_j12)A2)*((cot(teta_i))A2)* (etaxy(i,j,k-1)+etaxy(i,j-1,k-1));

vy(i,j,k)=vy(i,j,k)+resycur/kphicur*krelaxs;

end

end

if (k<znum+1)

if (i==1 Il i==ynum+1 Il j==1 Il j==xnum+1 Il k==1 II k==znum); vz(i,j,k)=RZ(i,j,k);

if(i==1)

vz(i,j,k)=2*RZ(i,j,k) -vz(i + 1,j,k);

end

if(i==ynum+1)

vz(i,j,k)=2*RZ(i,j,k) -vz(i-1,j,k);

end if (j==1)

vz(i,j,k)=2*RZ(i,j,k) -vz(i,j + 1,k);

end

if(j==xnum+1)

vz(i,j,k)=2*RZ(i,j,k) -vz(i,j-1,k);

end else

r_j=(j-1)*xstp+shift_r; r_j1=(j-2)*xstp+shift_r; r_j12=r_j-xstp*Q.5; teta_i=(i-1)*ystp+shift_teta; teta_i1=teta_i-Q.5*ystp; teta_i2=(i-2)*ystp+shift_teta;

reszcur=RZ(i,j,k) + (1/r_j12)*(1/sin(teta_i1))*(pr(i-1,j-1,k)-pr(i-1,j-1,k-1))/zstp;

%tau_r_phi

tau_r_phi = Q.5*(etaxz(i-1,j,k)+etaxz(i-1,j-1,k))* ((1/r_j12)*(1/sin(teta_i1))*zkf1*(Q.5*(vx(i,j,k+1)+

vx(i,j-1,k+1))-Q.5*(vx(i,j,k)+vx(i,j-1,k)))+Q.5*xkf1*(vz(i,j + 1,k)-vz(i,j-1,k))-(1/r_j12)*vz(i,j,k));

reszcur=reszcur-2*(1/r_j12)*tau_r_phi;

%dtau_r_phi/dr

dtau_r_phi = xkf1*etaxz(i-1,j,k)*((1/r_j)*

(1/sin(teta_i1))*zkf1*(vx(i,j,k+1) -vx(i,j,k))+xkf1*(vz(i,j + 1,k)-vz(i,j,k))-(1/r_j)*Q.5*(vz(i,j + 1,k)+vz(i,j,k)));

dtau_r_phi = dtau_r_phi-xkf1*etaxz(i-1,j-1,k)* ((1/r_j1)*(1/sin(teta_i1))*zkf1*(vx(i,j-1,k+1)-vx(i,j-1,k))+ xkf1*(vz(i,j,k)-vz(i,j-1,k))-(1/r_j1)*Q.5*(vz(i,j,k)+vz(i,j-1,k)));

reszcur=reszcur-dtau_r_phi;

%tau_teta_phi

tau_teta_phi = Q.5*(etaxz(i-1,j,k)+etaxz(i-1,j-1,k))* ((1/r_j12)*Q.5*ykf1*(vz(i+1,j,k)-vz(i-1,j,k))-

(1/r_j12)*vz(i,j,k)*cot(teta_i1)+(1/r_j12)*(1/sin(teta_i1))*zkf1* (Q.5*(vy(i,j,k+1)+vy(i-1,j,k+1))-Q.5*(vy(i,j,k)+vy(i-1,j,k))));

reszcur=reszcur-(1/r_j12)*cot(teta_i1)*tau_teta_phi;

%dtau_teta_phi/dteta

dtau_teta_phi = ykf1*etayz(i,j-1,k)*((1/r_j12)*ykf1*(vz(i+1,j,k)-vz(i,j,k))-( 1/r_j12)*cot(teta_i)*Q.5*(vz(i+1,j,k)+vz(i,j,k))+ ((1/r_j12)*(1/sin(teta_i)))*zkf1*(vy(i,j,k+1)-vy(i,j,k)));

dtau_teta_phi = dtau_teta_phi - ykf1*etayz(i-1,j-1,k)*

((1/r_j12)*ykf1*(vz(i,j,k)-vz(i-1,j,k))-

(1/r_j12)*cot(teta_i2)*Q.5*(vz(i,j,k)+

vz(i-1,j,k))+((1/r_j12)*(1/sin(teta_i2)))*zkf1*(vy(i-1,j,k+1)-vy(i-1,j,k)));

reszcur=reszcur-(1/r_j12)*dtau_teta_phi; %dtau_phi_phi/dphi

dtau_phi_phi = 2*zkf1*etan(i-1,j-1,k)*

(((1/r_j12)*(1/sin(teta_i1)))*zkf1*(vz(i,j,k+1)-vz(i,j,k)) + (1/r_j12)*Q.5*(vx(i,j,k+1)+vx(i,j-1,k+1)) + (1/r_j12)*cot(teta_i1)*Q.5*(vy(i,j,k+1)+vy(i-1,j,k+1)));