CC BY
54
НАУЧНО-ТЕХНИЧЕСКИИ ВЕСТНИК ИНФОРМАЦИОННЫХ ТЕХНОЛОГИИ, МЕХАНИКИ И ОПТИКИ январь-февраль 2016 Том 16 № 1 ISSN 2226-1494 http://ntv.i1mo.ru/
SCIENTIFIC AND TECHNICAL JOURNAL OF INFORMATION TECHNOLOGIES, MECHANICS AND OPTICS January-February 2016 Vol. 16 No 1 ISSN 2226-1494 http://ntv.ifmo.ru/en
УДК 517.95
ЭТАЛОННЫЕ РЕШЕНИЯ УРАВНЕНИЙ СТОКСА С ПЕРЕМЕННОЙ ВЯЗКОСТЬЮ В ЦИЛИНДРИЧЕСКИХ И СФЕРИЧЕСКИХ КООРДИНАТАХ ДЛЯ ТЕСТИРОВАНИЯ АЛГОРИТМОВ И.В. Макеев", И.Ю. Попов", И.В. Блинова"
a Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация Адрес для переписки: ilya.makeev@gmail.com Информация о статье
Поступила в редакцию 13.10.15, принята к печати 26.10.15 doi:10.17586/2226-1494-2016-16-1-161-167 Язык статьи - русский
Ссылка для цитирования: Макеев И.В., Попов И.Ю., Блинова И.В. Эталонные решения уравнений Стокса с переменной вязкостью в цилиндрических и сферических координатах для тестирования алгоритмов // Научно-технический вестник информационных технологий, механики и оптики. 2016. Т. 16. № 1. С. 161-167.
Аннотация
Рассматриваются стоксовы течения жидкости в областях со сферической и цилиндрической формой границы. Подобный выбор параметров системы является естественным при моделировании течений в задачах геофизики. Выводятся аналитические решения уравнений Стокса и неразрывности для частных случаев зависимости вязкости и плотности от цилиндрических координат. Данные решения задают класс осесимметричных течений, для которых вязкость является функцией радиуса. Выводятся аналитические решения уравнений Стокса и неразрывности в сферической системе координат для частного случая сферически симметричной вязкости. В работе показано, как на основе данных эталонных решений могут быть построены тестовые задачи для оценки качества работы численных алгоритмов. Приведены примеры тестирования многосеточных методов численного решения уравнений Стокса с переменной вязкостью в цилиндрической и сферической системе координат. Преимуществом данного подхода к использованию эталонных решений является возможность тестировать численные алгоритмы при различных перепадах вязкости и плотности. В работе предложена численная схема для многосеточного метода решения уравнений Стокса с переменной вязкостью в сферической системе координат. При построении решения используется последовательность ортогональных смещенных сеток. Качество численной схемы было проверено путем сравнения численного решения с аналитическим решением тестовой задачи. Ключевые слова
уравнения Стокса, переменная вязкость, метод Гаусса-Зейделя, эталонные решения, многосеточные методы
BENCHMARK SOLUTIONS FOR STOKES EQUATIONS WITH VARIABLE VISCOSITY IN CYLINDRICAL AND SPHERICAL COORDINATES FOR ALGORITHM TESTING I.V. Makeev", I.Yu. Popov", I.V. Blinova"
a ITMO University, Saint Petersburg, 197101, Russian Federation Corresponding author: ilya.makeev@gmail.com Article info
Received 13.10.15, accepted 26.10.15 doi:10.17586/2226-1494-2016-16-1-161-167 Article in Russian
For cit"tion: Makeev I.V., Popov I.Yu., Blinova I.V. Benchmark solutions for Stokes equations with variable viscosity in cylindrical and spherical coordinates for algorithm testing. Scientific and Technical Journal of Information Technologies, Mechanics and Optics, 2016, vol. 16, no. 1, pp. 161-167.
Abstract
Stokes flows in cylindrical and spherical geometry are considered. Such flows are rather natural for geophysics. We derive some exact particular solutions of Stokes and continuity equations for particular dependence of viscosity and density on cylindrical coordinates. These solutions correspond to axisymmetric flows for the case when viscosity is a function of radius. We suggest exact particular solutions of Stokes and continuity equations with variable viscosity and density in spherical coordinates for the case of spherically symmetric viscosity and density distributions. We demonstrate how these solutions can be used for creation of test problems suitable for benchmarking numerical algorithms. Examples of such benchmarking are presented. The advantage of this benchmarking approach is the ability to test numerical algorithms for variable viscosity and density gradients. We suggest numerical scheme of multigrid algorithm for solving Stokes and continuity equations with variable viscosity in a spherical coordinate system. Calculations are performed on a sequence of orthogonal staggered grids. The quality of the numerical scheme was verified by comparing the numerical solution with the analytical solution of the test problem. Keywords
Stokes equations, variable viscosity, Gauss-Seidel method, benchmark solutions, multigrid methods
Введение
Решения уравнений математических моделей, применяемых в геофизике для описания течения земной мантии, как правило, могут быть получены только при помощи численных алгоритмов [1, 2]. Если произвести сравнение с прямыми экспериментальными данными невозможно, то появляется необходимость в использовании методов, позволяющих оценивать качество работы численных алгоритмов. Одним из подходов является сравнение нескольких алгоритмов между собой [3, 4]. Другим способом проверки является сравнение результата работы алгоритма с эталонными аналитическими решениями некоторой тестовой задачи. В ряде работ рассматриваются такие решения для декартовой системы координат [5-10]. Значительно меньше внимания уделялось эталонным решениям в криволинейных системах координат [11, 12], в то время как с вычислительной точки зрения такие случаи являются наиболее сложными.
Целью настоящей работы является вывод решений системы уравнений Стокса и неразрывности в цилиндрической и сферической системах координат для случая, когда вязкость и плотность являются функциями радиуса. Следует отметить, что мы находим класс частных решений именно системы уравнений, а не краевых задач. На основе этих решений ставятся тестовые краевые задачи для оценки качества работы численных алгоритмов.
Среди конечно-разностных алгоритмов решения уравнений Стокса одним из наиболее эффективных является многосеточный метод [13-16]. В данной работе мы также предлагаем реализацию и производим тестирование численной схемы для многосеточного метода решения уравнений Стокса с переменной вязкостью в сферической системе координат.
Формулировка задачи
Уравнения Стокса и неразрывности для случая переменной вязкости и плотности имеют вид
У-ст = -рО, (1)
У(ру) = 0, (2)
где ст - тензор напряжений; р - плотность жидкости; О - напряженность гравитационного поля, V -вектор скорости. Задача построения эталонных решений сводится к нахождению частных решений системы уравнений (1)-(2) в цилиндрической и сферической системах координат для таких параметров системы, при которых вязкость, плотность и напряженность гравитационного поля являются функциями координаты г : п = п(г), Р = р(г), О = О(г).
Решения уравнений в цилиндрической системе координат
В цилиндрических координатах (г, ф, г) тензор напряжений ст имеет следующее компоненты:
ду
°гг =- Р + 2П-дГ", °ФФ =- Р + 2П
Г 1 ^ г дф г
ду,
о,, =-Р + 2Л—^, дг
(
о гф = п
1 дг.+.дуф_ ^
г дф дг г
ду, дуг
а" = п ЦТ
оф, = п
Г 1
дг г дф
Л
где п - вязкость жидкости, Р - давление жидкости.
Будем искать частные решения системы уравнений (1)-(2) (т.е. компоненты скорости и давление) в следующем виде: уг = уг (г), уф = уф (г), V, = у, (г), Р = Р(г). В этом случае уравнения (1) примут
вид
2пг-1 у'г + 2п'у; + 2пу'- 2пг-2уг - Р' = -рОг,
п' уф - г-1п'уф + пуф + г-1пуф - пг-2уф = -РОф,
п г-1у, + п' у, + пу, = - РО,.
Уравнение неразрывности (2) может быть представлено в форме рг-Ч + р'уг + ру'г = 0 .
Проинтегрировав уравнения, находим компоненты скорости и давление:
у г = С (гр)-1,
уф = С2/ (г) + Сз г + С1 (г) / (г) + С2 (г )г, (3)
г Г г, А
у, =-| (пг2 ) 1 | г\рОг^г1 + С
ёг2 + с5,
Р(г) = | (рОг + 2цг-1у' + 2ц'у'г + 2п\у"г - 2цг-2уг)"г ,
(4)
где
( (
/ (г) = ехр
1
V V
1 1
- + —
г2 Чг2
2 ,
I ъ"Г1 1 ПТ1
л-1 А I
+ С
"г2
у
С
1(г ) = |
грОф
П( / - / Т)
-(Зг , С
(Г) Ч
/рО
ф "г.
П( / - 1'г)
Формулы (3)-(4) задают осесимметричные течения для произвольных видов функций вязкости и плотности.
В качестве примера рассмотрим течение с линейной зависимостью вязкости от радиуса п = а1 г и постоянной плотностью р = а2 в области, заданной неравенствами 1 < г < 2,0 < ф < 1,0 < г < 1. Выберем частное решение (3)-(4) с параметрами а1 = 1, а2 = 1:
-1 -2 -1 -1 уг = г , уф = г + г , у2 = г , Р(г) = 2г , (5)
На рис. 1 показано поле скоростей для решения (5) в декартовой системе координат. Заметим, что построенные решения имеют особенность на оси цилиндра, значит, для тестирования можно брать лишь
области, не содержащие ось.
1
0,8 0,6 2 0,4 0,2 0
1,5
Vх \ Л
1,
У
0,5
1
0 0 0,5 я
Рис. 1. Поле скоростей, декартовы координаты Решения уравнений в сферической системе координат
В сферических координатах (г, 9, ф) компоненты тензора напряжений а выглядят следующим об-
разом:
агг = -Р + 2п ^, а99 = -Р + 2п (- ^ + - I, Чф = -Р + 2п дг V г д9 г
(
1
дуф у V
+ +-9(9)
г Бт(9) дф г г
(
аг9 = П
а9ф = П
г |+1 дУг дг V г ) г д9
( ат(9) д ( Уф
д9
(
агф = п
1
дуг
Л, \\
- + г-
ат(9)
1
г ат(9) дф дг
ду9
г
V ))
г Бт(9) дф
Будем искать частные решения системы (1)-(2) в следующем виде: Уг = уг(г), у9 = у9(г,9) уф = уф (г, 9), Р = Р(г). В этом случае уравнения (1) значительно упрощаются:
1 д ( о ду
I г22п— 1 +
г 2 дг V
— 2п
1
д
дг ) г аш(9) д9
1 дуе_
г д9
- 12П
Пг дг (у I ®ш(9)
дР
Л
с1§(9) | —— = -рО,
дг
д_
г 2 дг
1
+-п
г
г пг
д ( у,
дг V г
1
д
г аш(9) д9
2п|1 ^ + ^ |а^) г д9 г
\
д ( у.
дг V г
\
--2п| + -£■ с1м(9) |с1м(9) = 0,
г V г г
2
У
1 5
r п
5r
\\
уу
1
5
r sin(0) 59
5 ( V
+П —I — 5r ( r
+ n
sin(9) 5
59
( ( п
V v
sin(9) 5
59
ЛЛ
sin(9)
sin(9)
/у
sin(9)
ctg(9) = 0.
Уравнение неразрывности (2) примет вид
1 д-(pv0 sin(9)) = 0.
1 5 . 2 .
(pr 2Vr ) + -
г- дг г 8т(6) д9
Для случая, когда вязкость является степенной функцией радиуса п = С га , находим следующие частные решения уравнений (1)-(2):
V = -Т i P^M^ V9 = V91 ctg(9X Vp = r | c1 f-A
pr J V W
A
• + c2
1 5 (
о -I r22n
r2 5rV '
5 ( v
sin(9),
л
n—4^—2 + 2^-— v91
5r ) 5r V r ) r r
(6)
(7)
где v91 =
P(r) = f| pG,
ArC + BrC2 , C12 = 2|-(a +1) +7(a+1)2 + 4a j, a <-3 - 2^2 или a >-3 + 2^2 , A, B, С -
произвольные константы.
Ввиду особенности решений (6)-(7) возникает необходимость исключать точку r = 0 из области, в которой рассматривается течение.
Многосеточный метод
Численная схема для решения уравнений Стокса с переменной вязкостью в декартовой системе координат описана в [1]. Мы предлагаем аналогичную численную схему для решения уравнений Стокса с переменной вязкостью в сферической и цилиндрической системах координат. Данная схема позволяет находить численное решение уравнений (1)-(2) в общем случае, а не только ограничиваясь радиальной зависимостью параметров.
Как правило, итерации цикла построения решения многосеточных методов включают три основные операции: smoothing (сглаживание), restriction (проекция), prolongation (интерполяция). Сферические координаты являются ортогональными, поэтому реализация операций проекции и интерполяции не отличаются от их реализации для декартовых координат. Рассмотрим более подробно операцию сглаживания.
Метод Гаусса-Зейделя для сферических координат
Операция сглаживания может быть реализована на основе метода Гаусса-Зейделя. Введем трехмерную сетку для давления, вязкости, компонент скорости и плотности жидкости. Для уравнения Стокса и неразрывности получим следующую итерационную процедуру вычисления давления и компонент скорости в узлах сетки:
pnew _ р P(i, j ,l ) = P(i, j ,l ) +
ьр continuity p. continuity n(i, j,l) j,l 0relaxation ,
щопйшг^ = Rcontitrrnty - J_ ( r 2vr )--L_ 5- (v9 sin 9)--^ 5P (vp )
X j ,l
r 2 5r
r sin 9 59
r sin 9 5ф
AR
vr(i, j,l)- vr(i, j,l) ■
ve(i, j,l)- "e(i, j,l) ■
r -Stokes j,l
CV, (i, j,l)
a d9-Stokes
ARi, j,l ----I
с
Stokes ^ relaxation >
0 Stokes relaxation >
vnew = V +
vф(i, j,l ) = vф(i, j,l ) +
v9 <Л Л1 )
д рф- Stokes ARi,],l
с
0 Stokes relaxation '
Vp(ii jJ)
где 9 CXx7tZ 1 9 S'aation - коэффиЦиенты Релаксации; Cvr (i, j,lCv9 (i, jlCvp (i, j ,l)
(8) (9) (10)
являются коэффици-
ентами, на которые умножаются значения компонент скорости уг^,j^), у9(г-,j^), уф(г-,^) в конечно-
разностной записи уравнений Стокса. Для случая с переменной вязкостью остатки в выражениях (8)-(10) примут вид
2 5r
r
* d r- Stokes ARi,j,i _ or- Stokes 2 --T - rr r ~ (Trr )-1ctg(9)Tre-or r
-1-5- (t V r se(T re) 1 5T rp . Tee + Tpp + 5P
r sin(e) Sp r 5r
* ne-Stokes ARi, i ,i _ D e- Stokes Ri, j,l 2 _ Tre r -4-(Tre)-1ctg(e)Tee - or r
1 5 ( ) r se(Tee) 1 5Tep T re-T pp ctg(e) 1 5P
r sin(e) Sp r r 0e '
ARtp- Stokes _ nffl- Stokes Rr, j,l 2 --T Vp r (Trp )-1ctg(e)Tep-or r
1 5 ( ) r se(Tep ) 1 5Tpp Trp + Tep ctg(e) , 1 5P
r sin(e) Sp r r sin(e) 5p
где т - девиаторныи тензор напряжении.
Аналогичным образом может быть построена операция сглаживания для многосеточного метода в цилиндрической системе координат.
Проверка сходимости численных алгоритмов
Для проверки численной схемы используем сравнение численного решения с аналитическим решением уравнений (1)-(2). Рассмотрим частное решение (6)-(7): уг = г"^^/л/У, у9 = г^(9),
уф = г~58ш(9), Р(г) = 1,176г^ + рОгг (р = 1, Ог = 10, 09 = 0, Оф = 0, п = г3) в области, граница которой определяется неравенствами 1 < г < 2, 0,5 < 9 < 1,5, 0 < ф < 1. На границе области зададим значения для скоростей и давления, равные значениям аналитического решения. Внутри области решим уравнения Стокса численно, используя предложенную схему для многосеточного метода в сферической системе координат. В силу единственности решения численный результат должен совпадать с аналитическим. Выполнив сравнение этих двух решений, оценим качество численной схемы. Зависимость относительной ошибки е, вычисленной с нормой от величины шага сетки d (в логарифмическом масштабе) показана на рис. 2.
Аналогичным образом произведем тестирование подобной численной схемы многосеточного метода для цилиндрических координат. Сравним численное решение с аналитическим решением (5) в области 1 < г < 2, 0 <ф< 1, 0 < г < 1. На рис. 3 показана зависимость относительной ошибки е, вычисленной с нормой ^2, от величины шага сетки d (в логарифмическом масштабе). Положительный наклон линий графиков на рис. 2, 3 свидетельствует о сходимости численного решения к аналитическому.
Разумеется, предложенный метод тестирования не является всеобъемлющим (как, впрочем, и все остальные способы тестирования), а дает только необходимое условие правильности работы алгоритма. Проверка производится лишь на радиальной зависимости параметров и в областях, не содержащих начало координат в сферической системе и ось цилиндра в цилиндрической.
-1
ш ад о
-2
-3
-1,15 -1,05 -0,95 -0,85 -0,75 -0,65 log d
Рис. 2. Зависимость величины относительной ошибки s от величины шага: давление (P) - 1; vr - 2;
ve - 3; vp - 4
1
4
-1
-2
ш ьв о
-3
-4
-5
-1,55 -1,45 -1,35 -1,25 -1,15 -1,05 log d
Рис. 3. Зависимость величины относительной ошибки s от величины шага: давление (P) - 1; vr - 2;
Vp - 3; vz - 4
Заключение
В настоящей работе были выведены точные аналитические решения для ряда частных случаев системы уравнений Стокса и неразрывности с переменными вязкостью и плотностью в цилиндрической и сферической системах координат. Данные решения могут использоваться в качестве эталонов при тестировании численных алгоритмов решения уравнений Стокса с переменной вязкостью. Предложена и протестирована численная схема для многосеточного метода решения уравнений Стокса с переменной вязкостью в сферической системе координат.
Литература
1.
Gerya T. Introduction to Numerical Geodynamic Modeling. Cambridge, Cambridge University Press, 2010. 358 p.
Ismail-Zadeh A., Tackley P. Computational Methods for Geodynamics. Cambridge, Cambridge University Press, 2010. 313 p.
Deubelbeiss Y., Kaus B.J. Comparison of Eulerian and Lagrangian numerical techniques for the Stokes equations in the presence of strongly varying viscosity // Physics of Earth Planetary Interiors. 2008. V. 171. N 1-4. P. 92-111. doi: 10.1016/j.pepi.2008.06.023
Duretz T., May D.A., Gerya T.V., Tackley P.J. Discretization errors and free surface stabilization in the finite-difference and marker-in-cell method for applied geodynamics: a numerical study // Geochemistry, Geophysics, Geosystems. 2008. V 12. N 7. Art. Q07004. doi: 10.1029/2011GC003567
Blankenbach B., Busse F. et al. A benchmark comparison for mantle convection code // Geophysical Journal International. 1989. V. 98. N 1. P. 23-38.
Busse F.H., Christensen U., Clever R., Cserepes L., Gable C., Giannandrea E., Duillou L., Houseman G. et al. 3D convection at infinite Prandtl number in Cartesian geometry - a benchmark comparison // Geophysical and Astrophysical Fluid Dynamics. 1993. V. 75. P. 39-59. doi: 10.1080/03091929408203646 Popov A.I., Lobanov I.S., Popov I.Yu., Gerya T.V. Benchmark solutions for nanoflows // Nanosystems: Physics, Chemistry, Mathematics. 2014. V. 5. P. 391-399.
Popov I.Yu., Lobanov I.S., Popov S.I., Popov A.I., Gerya T.V. Practical analytical solutions for benchmarking of 2-D and 3-D geodynamic Stokes problems with variable viscosity // Solid Earth. 2014. V. 5. N 1. P. 461476. doi: 10.5194/se-5-461-2014
Tosl N., Stein C., Noack L., Huttig C., Maierova P., Samuel H., Davies D.R., Wilson C.R., Kramer S.C., Thieulot C., Glerum A., Fraters M., Spakman W., Rozel A., Tackey P.J. A community benchmark for viscoplastic thermal convection in a 2-D square box // Geochemistry, Geophysics, Geosystems. 2015. V. 16. N 7. P. 2175-2196. doi: 10.1002/2015GC005807.
10.Popov A.I., Gerya T.V., Lobanov I.S., Popov I.Yu. On the Stokes flow computation algorithm based on Woodbury formula // Nanosystems: Physics, Chemistry, Mathematics. 2015. V. 6. P. 140-145. doi: 10.17586/2220-8054-2015-6-1-140-145
11.Tackley P. Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid // Physics of the Earth and Planetary Interiors. 2008. V. 171. N 1-4. P. 7-18. doi: 10.1016/j.pepi.2008.08.005
2.
3.
4.
5.
6.
7.
9.
1
12.van Keken P.E., Currie C., King S.D., Behn M., Cagnioncle A., He J., Katz R.F., Lin S.-C., Parmentier E.M., Spiegelman M., Wang K. A community benchmark for subduction zone modeling // Physics of the Earth and Planetary Interiors. 2008. V. 171. P. 187-197. doi: 10.1016/j.pepi.2008.04.015
13.Brenner S.C. A nonconforming multigrid method for the stationary stokes equations // Mathematics of Computation. 1990. V. 55. P. 411-437. doi: 10.1090/S0025-5718-1990-1035927-5
14.Kang K., Kwak D., Yon Y. The analysis of multigrid method for nonconforming method for the stationary stokes equations // Bulletin of the Korean Mathematical Society. 1996. V. 33. N 3. P. 343-357.
15.John V., Tobiska L. A coupled multigrid method for nonconforming finite element discretizations of the 2D-stokes equation // Computing. 2000. V. 64. N 4. P. 307-332.
16.Wang M., Chen L. Multigrid methods for the stokes equations using distributive Gauss-Seidel relaxations based on the least squares commutator // Journal of Scientific Computing. 2013. V. 56. N 2. P. 409-431. doi: 10.1007/s10915-013-9684-1
Макеев Илья Владимирович
Попов Игорь Юрьевич
Блинова Ирина Владимировна
аспирант, инженер-исследователь, Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация,
ilya.makeev@gmail.com
доктор физико-математических наук, профессор, заведующий кафедрой, Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация, popov1955@gmail.com кандидат технических наук, доцент, доцент, Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация, ит-a@yandex.ru
Ilya V. Makeev Igor Yu. Popov Irina V. Blinova
postgraduate, engineer, ITMO University, Saint Petersburg, 197101,
Russian Federation, ilya.makeev@gmail.com
D.Sc., Professor, Head of Chair, ITMO University, Saint Petersburg,
197101, Russian Federation, popov1955@gmail.com
PhD, Associate professor, Associate professor, ITMO University,
Saint Petersburg, 197101, Russian Federation, irin-a@yandex.ru
