Численное моделирование естественной конвекции в сферическом слое Текст научной статьи по специальности «Физика»

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2008 Математика и механика № 3(4)

УДК 536.25(32)

А.В. Гореликов, А.В. Ряховский ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЕСТЕСТВЕННОЙ КОНВЕКЦИИ В СФЕРИЧЕСКОМ СЛОЕ

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

Ключевые слова: численное моделирование, естественная конвекция, вращение, параллельная реализация.

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

Постановка задачи и математическая модель

Полагалось, что вязкая несжимаемая жидкость заполняет сферический слой толщиной H (рис. 1). Границы сферического слоя твердые и вращаются с постоянной угловой скоростью ю вокруг оси z. Температура нижней сферы Th, больше температуры верхней - Tc. Поле тяжести обладает центральной симметрией. Центр поля тяжести находится в центре сферического слоя, заполненного жидкостью. Модуль ускорения свободного падения линейно зависит от расстояния до силового центра:

g(r) = (g(r0 + H ) - g(r0))(r - r0) / H + g(r0).

В работе приняты следующие обозначения: (r, 0, ф) - сферические координаты; t - время; (ur, щ, иф) - компоненты скорости в сферических координатах; p -отклонение от гидростатического давления; T - температура; а - коэффициент температуропроводности; ß - коэффициент объемного расширения; рс - плотность при температуре Tc; v - кинематическая вязкость.

Безразмерные переменные определялись по формулам

t = vt /H2; R = r /H; U = Hur /v ; V = Hue /v ; W = Ниф /v ;

Ф = (T - Tc)/(TÄ - Tc) ; P = H2 p /pcv2.

Безразмерные параметры задачи:

А = г0 /Н - аспектное отношение; Gr = £(г0)Р(Гй - ГС)Н /у2 - число Грасгофа; Ле = юН2/у - число Рейнольдса, Рг = у/а - число Прандтля; у = (£(г0 + Н ) -я(г0)) / я(г0) - безразмерный параметр, характеризующий степень неоднородности модуля поля тяжести.

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

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

= -УР + Уа + (у(Я - А) +1 )Ог Ф ег; (1)

Б т

У У = 0; (2)

Бф 1 АЛ,

----= Дф ; (3)

Б т Рг

А < Л < А+1, 0 < 0 < п.

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

Граничные условия:

Л = А: Ф = 1, и = V = 0, Ж = Ле А Бт0, (4)

Л = А + 1: Ф = 0, и = V = 0, Ж = Ле (А + 1) 8т0.

Начальные условия (т = 0):

и = V = 0, = Ле Л Бт0, (5)

а) Ф = 0, б) Ф = 1, в) Ф = А(А+1)/Л - А.

Число Нуссельта вычислялось по формуле

1 л ЭФ

N = Nm /(Nm|& = 0), Nu =-----f-----dS ,

Pr S d R

S - верхняя или нижняя граница сферического слоя.

Поскольку число Грасгофа было определено через g(r0), то для интерпретации результатов определялось среднее по области число Грасгофа

Or '=— | Gr (у (R - A) +1 )dV,

VD D

где D - сферический слой, VD - объем сферического слоя.

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

W = W - Re R sin0.

Численная реализация

На языке Фортран 90 создан программный комплекс, позволяющий численно решать задачи конвекции вязкой несжимаемой жидкости в сферическом слое, как в осесимметричной (2D), так и в полной постановке (3D). При численной реализации использовался метод контрольного объема в сферических координатах [4]. Дискретные аналоги уравнений (1), (3) получены с использованием полностью неявной схемы и схемы со степенным законом для аппроксимации конвективных и диффузионных потоков на гранях контрольных объемов. Расчёт поля течения производится по модифицированному полунеявному методу для связывающих давление уравнений SIMPLER [4]. Полученные в результате дискретизации СЛАУ решались методом переменных направлений. Программный код распараллелен с помощью директив стандарта OpenMP [5, 6] для использования на многопроцессорных вычислительных системах с общей памятью.

Для использования директив распараллеливания OpenMP при решении алгебраических уравнений последовательный вариант метода переменных направлений потребовал модификации, поскольку выполнение первоначального кода параллельными потоками приводило к неадекватной работе программы и ошибкам. Ошибки возникали из-за состояния «гонки» (race condition) между параллельными потоками, решающими системы уравнений вдоль соседних сеточных линий. Состояние гонки возникает из-за того, что при решении СЛАУ вдоль одной сеточной линии для расчёта свободных членов уравнений используются значения искомой величины с предыдущей итерации в соседних сеточных линиях. Поэтому, если два потока обрабатывают две соседние линии узловых точек, то это может привести к тому, что они одновременно будут запрашивать на чтение и запись одни и те же переменные, что в итоге приводит к неопределенному результату.

Модификация метода переменных направлений заключалась в том, что расчёт свободных членов в СЛАУ производится для всех контрольных объёмов перед выполнением прогонки, при этом используются значения с предыдущей итерации. При таком подходе во время работы параллельных потоков переменные, содержащие значения свободных членов, только считываются. Программный код с описанной модификацией без труда поддаётся распараллеливанию директивами стандарта OpenMP и исключает возможность состояния «гонки».

Программа тестировалась на задачах, имеющих аналитическое решение. Далее приведены результаты одного из тестовых расчетов. В области (0 < R < 2, 0 < 0 < п ) рассматривалась система (1) - (3) при Pr = 1, Gr = 0. Вращение отсутствовало (W = 0).

Граничные условия:

Ф|д = i = 0,5( e2 - exp(cos0) ) / sh2, Ф|д = 2 = 0,5( e2 - exp(2 cos0) ) / sh2,

U\r = 1 = U\r = 2 = cos0, V\r = 1 = V\r = 2 = -sin0.

Аналитическое решение:

Ф0 = 0,5( e2 - exp(R cos0) ) / sh2, U0 = cos0, V0 = -sin0, P0 = const.

Точность численного решения оценивалась по максимальным абсолютным и относительным ошибкам величин Ф, U, V:

Лф = max |Ф-Ф0|, 8Ф = -100%,

— 1...

где Ф =— J|Ф0|dV - среднее по расчетной области D абсолютное значение

VD D

аналитического решения Ф0. Аналогично для U и V.

Результаты тестов представлены в табл. 1 (n и m - количество расчетных точек по координатам R и 0 соответственно).

Таблица 1

n х m Ли §№ % Лр бр, % Лф 8ф, %

12 х 32 4,3949-10 4 8,7897-10 2 2,0212-10 3 2,5735-10 1 1,1921-10 3 1,4714-10 1

22 х 62 1,6454-10 4 3,2908-10 2 5,1711-10 4 6,5840-10 2 3,0853-10 4 3,8084-10 2

42 х 122 7,0334-10 5 1,4067-10 2 1,4147-10 4 1,8012-10 2 7,8400-10 5 9,6776-10 3

Расчеты проводились на сетке с 82 расчетными точками по радиусу и 242 по углу. В некоторых случаях использовались более подробные сетки со сгущениями у границ слоя. Скорость созданного многопоточного кода на компьютере с четы-рехядерным процессором Intel Core 2 Quad Q6600 с тактовой частотой 2,4 ГГц, при расчете на сетке 82*242 увеличивается в 2,2 раза по сравнению со скоростью выполнения однопоточного кода.

Результаты расчетов

Все расчеты проводились при фиксированных значениях числа Прандтля Pr, аспектного отношения A и степени неоднородности поля тяжести у: Pr = 1, A = 0,54, у = 16,7. При таких значениях A и у среднее по области число Грасгофа определяется соотношением Gr' ~ 11,8 Gr.

В первой серии численных экспериментов исследовалась зависимость стационарных решений от числа Грасгофа при Re = 10. Для Gr = 103, 104, 105 при расчетах использовались все варианты начальных условий (5). В результате, для Gr = 103 вне зависимости от типа начальных условий получалось одно и то же стационарное решение, представленное на рис. 2. В данном случае, течение в меридиональной плоскости имеет четырехвихревую структуру с восходящими потоками у оси вращения и в плоскости экватора (0 = п /2). При Gr = 104 для всех вариантов начальных условий также был получен только один тип решения (рис. 3), но с двухвихревой меридиональной циркуляцией с восходящим потоком в плос-

кости экватора. Качественно типы решений могут быть классифицированы по структуре температурного поля в меридиональной плоскости. В частности, полученные решения можно разделить на два типа: первый тип - с одним термиком в экваториальной области (рис. 3), второй - с тремя термиками, один в области экватора и два вдоль оси вращения (рис. 2). Для Gr = 105 при начальном условии б) в (5) получался первый тип решения, а при начальных условиях а), в) - второй тип решения (рис. 4), но при этом размер вихрей у оси вращения существенно уменьшается по сравнению с решением при Gr = 103.

а б в

Рис. 2. Поле температуры (а), структура течения (б) и осевая составляющая скорости Ж' (в) в сечении по ф (От = 103, Яе = 10)

а б в

Рис. 3. Поле температуры (а), структура течения (б) и осевая составляющая скорости Ж' (в) в сечении по ф (От = 104, Яе = 10)

а.2 б.2 в.2

Рис. 4. Поле температуры (а), структура течения (б) и осевая составляющая скорости Ж' (в) в сечении по ф (От = 105, Ле = 10)

Для всех остальных значений числа Грасгофа в диапазоне от 5102 до 5-105 стационарные решения получались методом продолжения по параметру, т.е. в качестве начальных условий использовалось стационарное решение при близком значении Gr. Решения первого типа были получены для всех чисел Грасгофа в исследованном множестве значений. Решения второго типа не удалось получить только в промежутке 104 < Gr < 4-104. Таким образом, при медленном вращении (Ле = 10) возможно существование, по крайней мере, двух различных типов стационарных режимов конвекции в достаточно широком диапазоне значений числа Грасгофа.

Конвекция во всех рассмотренных случаях имеет сложную трехмерную структуру, на вращение жидкости вокруг оси г накладывается меридиональная цирку-

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

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

N = 0,15^' °’25, (6)

которому соответствует кривая 3 на рис. 5. При Gr > 5,5-105 не удалось получить стационарного режима конвекции.

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

N 7

6

5

4

3-К

Л

2

1-106 2 1 06 3 106 4 106 5 106 вт'

Рис. 5. Зависимость числа Нуссельта N от среднего по области числа Грасгофа От': 1 -решения первого типа, 2 - решения второго типа, 3 - соотношение (6)

Рис. 6 Зависимость числа Нуссельта N от Ле (От = 105).

1

В третьей серии расчетов исследовалась зависимость стационарных решений от числа Рейнольдса при фиксированном числе Грасгофа Gr = 105. Использовался метод продолжения по параметру, т.е. в качестве начальных условий использовалось решение при предыдущем значении Ле. В диапазоне 0 < Ле < 45 для каждого значения Ле получены оба типа стационарных решений, при 50 < Ле < 85 - только решения второго типа. Зависимость числа Нуссельта от числа Рейнольдса представлена на рис. 6. Тепловой поток уменьшается при увеличении числа Рейнольдса, поскольку при увеличении Ле усиливается центробежная сила и, как следствие,

уменьшается интенсивность конвекции. Результирующее силовое поле складывается из гравитации и центробежной силы, и при достижении некоторого значения Ле (в данном случае при Ле = 95) становится сильно неоднородным, и, как следствие этого, конвекция становится нестационарной. При больших числах Рейнольдса (в данном случае при Ле > 3000) конвекция прекращается, поскольку центробежная сила полностью компенсирует гравитацию. Аналогичное исследование было проведено для Gr = 103. В этом случае стационарное решение получается при любом значении Ле и уже при Ле = 100 конвекция прекращается.

Заключение

Разработан комплекс программ, позволяющий численно решать задачи конвекции в сферическом слое. Проведено распараллеливание программного кода с помощью директив стандарта ОрепМР для использования на многопроцессорных вычислительных системах с общей памятью. Ускорение выполнения типичного для данной работы расчета при использовании многопоточного кода на процессоре с четырьмя ядрами составило 2,2 раза. Исследована зависимость стационарных решений задачи о естественной конвекции в равномерно вращающемся сферическом слое от чисел Грасгофа и Рейнольдса. Установлено, что при медленном вращении (Ле = 10) возможно существование, по крайней мере, двух различных типов стационарных режимов конвекции, соответствующих одному значению числа Грасгофа. Установлено также, что при увеличении интенсивности вращения (Ле = 50) снимается вырождение стационарных решений, т.е. каждому числу Грасгофа соответствует только одно стационарное решение. Определены зависимости числа Нуссельта от среднего по области числа Грасгофа и числа Рейнольдса.

ЛИТЕРАТУРА

1. Яворская И.М., Беляев Ю.Н. Конвективные течения во вращающихся слоях // Механика жидкости и газа. Т. 17. М.: ВИНИТИ (Итоги науки и техники), 1982. С. 3 - 85.

2. Кирдяшкин А.Г., Добрецов Н.Л, Кирдяшкин АА. Турбулентная конвекция и магнитное поле внешнего ядра Земли // Геология и геофизика. 2000. Т. 41. № 5. С. 601 - 612.

3. Аплонов С.В. Геодинамика. СПб.: Изд-во С.-Петерб. ун-та, 2001. 360 с.

4. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984. 152 с.

5. Немнюгин МА., Стесик О.Л. Современный Фортран. Самоучитель. СПб.: БХВ-Петер-бург, 2004. - 496 с.

6. Горелик А.М. Программирование на современном Фортране. М.: Финансы и статистика, 2006. 352 с.

Статья принята в печать 27.10.2008 г.