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

Вычислительные технологии

Том 5, № 5, 2000

РАСЧЕТ МИКРОКОНВЕКЦИИ В ДЛИННОМ ПРЯМОУГОЛЬНИКЕ *

О. Н. Гончарова Новосибирский государственный университет, Россия Институт гидродинамики им. М. А. Лаврентьева СО РАН

Новосибирск, Россия e-mail: olga@person4.math.dcn-asu.ru

Non-stationary convection is studied numerically. The method of alternating directions is used for the solution of the equations written for the variables “stream function” and “vortex”. The calculations are carried out for the real Prandtl, Rayleigh and Boussinesque numbers.

Введение

Система уравнений Обербека — Буссинеска является традиционной математической моделью для описания тепловой гравитационной конвекции и представляет собой результат упрощения полных уравнений механики сплошных сред при условиях [1]:

1. Плотность жидкости р зависит лишь от ее температуры Т (жидкость изотермически несжимаема), и р = ро(1 — вТ) (в — коэффициент объемного расширения, р0 — некоторое характерное значение плотности).

2. Движение подобно движению несжимаемой жидкости, так что поле скоростей считается соленоидальным. Однако, в уравнении импульса приближенно учитывается изменение плотности.

3. Вклад диссипации и сил давления в описание притока тепла пренебрежимо мал.

4. Динамический коэффициент вязкости р, коэффициенты теплопроводности k и удельной теплоемкости cp предполагаются постоянными.

Для исследования конвекции в областях малой протяженности, в слабых гравитационных или быстропеременных температурных полях В. В. Пухначев [2, 3] предлагает использовать новую, уточненную, модель. Новая модель получила название модели микроконвекции для характеристики течений с малыми значениями параметра n = gl3/vx, где l — характерный размер области, занятой жидкостью, v, x — коэффициенты кинематической вязкости и температуропроводности, g = \g\, g — ускорение силы тяжести. Модель микроконвекции исходит из точных уравнений неразрывности и импульса. Вместе с тем, уравнение энергии упрощается в духе гипотезы (3), и для простоты считается выполненной гипотеза (4). Гипотеза же (1) принимается в виде р = р(Т). Новая модель характеризуется

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №99-01-00529, и Совета Программы поддержки ведущих научных школ, грант №00-15-96162.

© О. Н. Гончарова, 2000.

несоленоидальностью поля скоростей V. Однако, в случае зависимости плотности от температуры р = р0/(1 + вТ) полученная система уравнений преобразуется к виду, в котором модифицированный вектор скорости V = V — вхУТ уже становится соленоидальным. Соленоидальность модифицированной скорости IV позволяет ввести аналог функции тока для плоских и осесимметричных задач и выполнять расчеты конвективных течений в переменных “функция тока — вихрь”.

Зависимости плотности от температуры, принятые в модели Обербека — Буссинеска и в модели микроконвекции, практически равносильны, т. к. в реальных конвективных течениях максимальные значения величины в|Т| не превосходят 10-2. Параметр же п является дополнительным по отношению к числам Рэлея И,а и Прандтля Рг = и/х критерием подобия. При этом И,а = п ■ Ви, где Ви = вТ* — так называемое число Буссинеска, а Т* — характерный перепад температуры. С физический точки зрения параметр микроконвекции равен отношению порядков скоростей, порожденных объемным расширением жидкости (вТ*х/1) и фактором плавучести (вТ*12д/и).

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

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

Задача математического моделирования микроконвекции состоит в нахождении модифицированной скорости IV, модифицированного давления д, температуры Т, удовлетворяющих в области течения системе уравнений [2]:

#*+#-V#+вх(УТ -V#—V#-УТ )+в 2Х2(АТ-УТ—VIУТ |2) = (1+вТ )(—Уд+и Д#)— вдТ,

= 0,

Здесь Е — граница области течения, д = р — в(и + и' — х)хДТ, р = р/р0 — д ■ X, р — давление, и' — так называемая вторая вязкость.

Начально-краевая задача для системы Обербека — Буссинеска состоит в определении физической скорости V, модифицированного давления рр, температуры Т, удовлетворяющих системе уравнений [1]:

Т + IV ■УТ + вх|УТ|2 = (1 + вТ)хДТ,

(1)

начальным условиям (при £ = 0)

V = #о(ж), Т = То (х)

У + V - УV = —Ур + иДУ — вдТ,

&уУ = 0,

Т + у ■ УТ = хЛТ, (2)

начальным

V = 0, Т = То(х)

и граничным условиям

дТ

у = 0 дП =1 (х^

Нестационарная конвекция исследуется в области 0 < х < х0, 0 < у < у0 с твердыми непроницаемыми границами, две из которых (х = 0, х = х0) теплоизолированы, а на оставшихся задан периодический тепловой поток. Сила тяжести направлена по оси Ох. Уравнения конвекции и микроконвекции (1), (2) рассматриваются в переменных ш — вихрь скорости, ф — функция тока или модифицированная функция тока (при этом VI, г>2 — компоненты физической или модифицированной скорости):

VI = фу, ^2 = — фх,

Ш + VI + ^2^у = />Дш + вдТу + , (3)

Дф = —ш, (4)

Ті + ^1Тх + ^2Ту = ХДТ + ^Т • (5)

При этом для модели микроконвекции имеют место следующие соотношения:

/>=(1+ вт )^, х = (1 + вт )х,

= в ( — Тх?у + Ту?х) + Vе (Д^2Тх — Д^1Ту) +

+ ( —вХ) (ш ДТ + УТ ■ Уш + 2^2уТху + ^1уТхх — ^2хТуу) + ( —^2Х2) (ДТхТу — ДТуТх) ,

Рт = -вх|УТ |2.

Для модели Обербека — Буссинеска имеем, соответственно

> = V X = X, = 0 рт = 0-Граничные условия для модели микроконвекции:

х = 0 : ф = 0, фх = вхТу, Тх = 0, (е;

х = хо : ф = вххоА від 7і, фх = вхТу, Тх = 0, (7;

У = 0: ф = вххА від 7І, фу = —вхТх, Ту = А від 7^, (8

У = 1 : ф = вххА від 7І, фу = —вхТх, Ту = А від т*; (9

для модели Обербека — Буссинеска:

0, Т=х 0, х ф 0, ф= 0 х= (іо;

х х О = 0, ф 0, Т 0, (іі

ід Ту 0, у ф 0, ф= 0 У= (12;

40 ід Ту ,0 у ф ,0 ф= 1 У= (із;

Начальные условия для обеих моделей:

Т 0, ф= 0, 3 0 (14;

2. Численное исследование

Численное исследование задач (3) - (9), (14) и (3) - (5), (10) - (14) проводится с использованием продольно-поперечной конечно разностной схемы (метод переменных направлений), формально имеющей второй порядок аппроксимации [5, 6].

Для уравнений (3), (5) схема расчета записывается в следующем общем виде:

фк+^2~ = [>Фх - «1Ф]Х + [>Ф„ - «2Ф]Г/2 + Рк,

Фк+1 -/ф*+1/2 = [>Фх - «1Ф]Х+1 + [>Фу - «зФ]^2 + ^• (15)

Для модели микроконвекции, если Ф = ш, то

Р = вдТу + - в^УТ ■ Уш,

если Ф = Т, то

рг = Рт - вх|УТ|2.

Для модели Обербека — Буссинеска, если Ф = ш, то

Р = вдТу + ,

если Ф = Т, то

Р = Рт •

Для решения уравнения (4) на каждом временном шаге ^ = кт, к = 1, 2, ... применяется итерационная схема:

ф*+1/2 - ф*

т/2

ф* + 1 — ф* + 1/2

Л [ФХ+1/2 + фуу + шк+^,

Л[ф*+1/2 + ф*+1 + шк+1], (16)

т/2

с итерационным параметром Л.

Для реализации вышеизложенной схемы расчета вводится разностная сетка:

х„ = (п - 1)Ьх, п = 1, 2, ••., N1, ^х = N, N1 = N +1,

Ут = (т - 1)йу, т = 1, 2, •••, М1, Л,у = , М1 = М + 1,

Уо

М;

\к

Ф(^к, хп, ут) — Ф га,т •

Системы линейных алгебраических уравнений (15), (16) представляются в виде — п Фк+1/2 + Ь Фк+1/2 _ с Фк+1/2 = А п = 2 N

'Л'га,т^П)т- 1 "Т" п,т °га,та^п,т+1 и'га,т; •••)

-пга,тФП+1,т + Ьга,тФП+ог - Сга,тФП+1,т = т = 2 •••, M,

и решаются методом прогонки [7]. При этом используется следующая разностная аппроксимация для входящих в (15) слагаемых

[>Фх - фуФ] х ~

^n+1

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

При этом для аппроксимации граничных условий для температуры, задающих тепловой поток (8), (9), используются дополнительно сами уравнения переноса тепла (5).

Граничные условия для вихря задаются соотношениями типа Тома или Вудса [8], которые в модели микроконвекции имеют вид при у = 0:

туры Tfc+1, а потом — вихря скорости wfc+1 по схеме (15). При этом на промежуточном (k + 1/2)-м слое осуществляется прогонка в направлении у, а на основном (k + 1)-м слое — прогонка в направлении х. (Стартовые значения определяются начальными данными

T := T0 = const, и := 0, ф := 0.)

2. На каждом k-м временном слое вводится внутренний итерационный процесс расчета ф5+1 из систем уравнений (16) с измененной очередностью прогонок. После окончания итераций при s = S считается, что с заданной точностью определены значения ф на (k + 1)-м временном слое: фк+1 = ф5.

Итерационный процесс считается сошедшимся, если выполнен критерий сходимости [8]:

где в — номер итерации, — заданная точность расчета ф5+1.

Точность выполнимости граничных условий для вихря определяется, например, для условий (17), (18), по величине [8]:

k+1/2 1 k+1/2 и ^n,l +2 Wn,2 :

а для модели Обербека — Буссинеска:

(1S)

соответственно.

Общая схема решения задачи состоит в осуществлении следующих этапов.

1. Алгоритм перехода на новый временной слой ік+1 начинается с расчета темпера-

max ^n+m - "0n,m1 <^max ад,

є = max І^дС^1) - ^п,і(ФП,2)і-

n

3. Результаты численного анализа

Расчеты выполняются на сетке 200 х 20 для области 0 < х < х0, 0 < у < у0 при х0 = 10 см, уо = 1 см.

Устойчивость же алгоритма и порядок сходимости проверяются путем вычислительного эксперимента на последовательностях сеток (г = 1 : 200 х 20, г = 2 : 400 х 40, г = 3 : 800 х 80), при этом ведутся наблюдения за величиной г^, характеризующей для каждой сетки с номером г интенсивность движения тах |фга,т| или являющейся интегральной нор-

мой

= . Экспериментальный порядок сходимости г и приближенно

определенная относительная погрешность е вычисляются согласно правила Рунге (см., например, [9]):

1п(|г2 - Г1 |/|Гз - Г2 |) 1 | Г1 - Г2 |

1п 2

1 -

Г 2

По этой формуле вычисляются г 1.8 и е 3 %, в том случае, если основной величиной

наблюдения является интенсивность движения. Если же таковой является интегральная Ь2-норма для функции тока, то г ~ 2, е ^ 3 %.

Основные значения параметров в граничных условиях для температуры: Т0 = 35, А = 70, 7 =2 или 7 = 0.5. В табл. 1 приведены основные параметры задачи. Размерные величины даны в системе СГС.

Таблица 1

Основные параметры

Рг Ие П є V X 9 в 7

N.1 0.75 0.014 1 0.01 0.15 0.2 0.03 0.0003 2(0.5)

N.2 0.01 2.1 0.4 0.02 0.015 1.5 0.009 0.0006 2(0.5)

N.3 0.1 0.21 0.4 0.02 0.15 1.5 0.09 0.0006 2(0.5)

N.4 (эИ) 0.0054 0.049 0.77 0.000026 0.00265 0.49 0.001 0.0000075 2(0.5)

Таким образом, расчеты проводятся для нескольких модельных жидкостей (N.1, N.2, N.3, N.4) при действии микроускорений, достигаемых на орбитальной станции, в условиях малости параметра п, рассматриваются при различных значениях д. Модели различаются по числу Рг и И,е = V*//V. При этом в качестве характерной скорости выбирается характерная скорость объемного расширения V* = вТ*х//, / = у0 = 1- В качестве характерного времени процесса естественным образом может быть выбрано значение £* = //V*. Следует отметить, что определение жидкостей, как модельных, довольно условно. Модель N.4 представляет собой расплав кремния, а модель N.1 — воздух (по значениям основных параметров), часто используемый для проведения реальных экспериментов по конвекции.

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

И = ^1’ И = ^2’ х(0) = х’ у(0) = у>

рассчитываются по улучшенному методу Эйлера второго порядка точности [9]. Здесь ^1, г>2 — компоненты физической скорости.

п, т

г

На рис. 1-5 представлены траектории жидкой частицы вплоть до момента времени і = 240 с, за исключением рис. 1, г, представляющего траектории на момент времени і = 60 с. При этом на осях отложены значения, умноженные на 104.

Рисунки 1, а, 2, 3, 4 представляют качественное различие в траекториях для разных моделей N.1, N.2, N.3, N.4 для частицы, находящейся в начальный момент времени в одной и той же точке.

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

Рис. 1. Траектории жидкой частицы. Модель N.3.

а — начальные координаты — (5,0.6), параметры расчета — Гг = 0.1, е = 0.02, п = 0.4, 7 = 2; б — начальные координаты — (5,0.6), параметры расчета — Рг = 0.1, е = 0.02, п = 0.4, 7 = 0.5; в — начальные координаты — (5,0.1), параметры расчета — Рг = 0.1, е = 0.02, п = 0.4, 7 = 2.

Рис. 1 (окончание). Траектории жидкой частицы. Модель N.3. г — начальные координаты — (5,0.1), параметры расчета — Рг = 0.1, е = 0.02, п = 0.4, 7 = 2.

На рис. 1, г представлены траектории частицы, находящейся в начальный момент в точке (5, 0.1), для сравнения влияния амплитудного коэффициента А в граничном тепловом режиме. Естественно, чем больше коэффициент А, тем более интенсивным является движение.

Траектории, рассчитанные с использованием модели Обербека — Буссинеска, заполняют отрезки прямых, параллельных Ох, и отмечены жирно. Траектории, рассчитанные по модели микроконвекции, очевидно, имеют более сложную структуру.

Итак, помимо отличия, диктуемого различными математическими моделями, и отличия в поведении траекторий в зависимости от угловой частоты (7) и амплитудного коэффициента А в граничном тепловом потоке наблюдаются различия в траекториях частиц, находящихся в начальный момент времени в разных точках (см. рис. 1, а, 1 в, 5, а, 5, б). При этом на рис. 5, а и 5, б (нижние диаграммы соответствуют модели Обербека — Буссинеска) для точки, расположенной вблизи теплоизолированного торца, наблюдаются качественно похожие траектории, но существенна разница в амплитудах, вызванная различием в порядках скоростей, в направлении движения (см. рис. 5, а) и в угле наклона (см. рис. 5, б).

Кроме качественного отличия в траекториях, получаемых с использованием разных

Рис. 2. Траектория жидкой частицы. Модель N.1. Начальные координаты — (5,0.6), параметры расчета — Pr = 0.75, е = 0.01,

П =1, Y = 2.

Рис. 3. Траектория жидкой частицы. Модель N.2. Начальные координаты — (5,0.6), параметры расчета — Рг = 0.01, е = 0.02,

П = 0.4, 7 = 2.

Рис. 4. Траектории жидкой частицы. Модель N.4 (зг/).

Начальные координаты — (5,0.6), параметры расчета — Рг = 0.0054, е = 0.000026, п = 0.77, 7 = 2.

математических моделей, наблюдаются и количественные различия для скоростей. Данные результаты представлены в табл. 2 наряду с указанием характерных скоростей и времени движения для различных моделей. Скорости, рассчитанные по модели микроконвекции (ММ), на 1-2 порядка превышают те, что предписываются моделью Обербека — Буссинеска (ОВМ). Значения характерных скоростей подтверждают расчетные величины.

Рис. 5. Траектории жидкой частицы. Модель N.3. а — начальные координаты — (9.75, 0.5), параметры расчета — Рг = 0.1, е = 0.02, п = 0.4, 7 = 2; б — начальные координаты — (0.25, 0.5), параметры расчета — Рг = 0.1, е = 0.02, п = 0.4, 7 = 2.

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

Поля температуры, рассчитанные по обеим моделям, практически одинаковы по своим количественным и качественным характеристикам (табл. 3). Изотермы имеют характерный для конвекции вид и практически параллельны длинным сторонам прямоугольника (рис. 6), что подтверждает предположение о доминирующем изменении температуры в направлении Оу и незначительном — в направлении Ох.

На рис. 6,7, а, б (две верхних диаграммы) представлены также линии тока для модели N. 3, рассчитанные по модели Обербека — Буссинеска и по модели микроконвекции соот-

1.0

0.0 5.0 10.0

Рис. 6. Поля температуры (изотермы). Модель N.3.

Параметры расчета — Рг = 0.1, е = 0.02, п = 0.4, 7 = 2.

0.0 5.0 ю.о

Рис. 7. Линии тока. Модель N.3.

Параметры расчета — Рг = 0.1, е = 0.02, п = 0.4, 7 = 2. Модель Обербека — Буссинеска (а),

модель микроконвекции (б).

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

Таблица 2

Абсолютные величины скоростей (максимальные порядки):

(ММ/ОВМ)

і* V* N ы

N.1 476 0.0021 1 0 1 1 О 1 10-3/10-5

N.2 31 0.0315 1 0 1 со 1 О 1 1 0 1 ю 1 О 1

N.3 31 0.0315 1 0 1 ьс 1 О 1 со 1 0 1 ьс 1 О 1

N.3 (7 = 0.5) 31 0.0315 1 0 1 со 1 О 1 1 0 1 ю 1 О 1

N.4 (ей) 7774 0.000129 1 0 1 Сл 1 О 1 1 0 1 1 о 1

Таблица 3

Изменение температуры

Т (ММ) Т (ОВМ)

N.1 11.02 - 59.83 10.65 - 59.30

N.2 11.27 - 60.69 10.24 - 59.76

N.3 11.27 - 60.70 10.24 - 59.76

N.3 (7 = 0.5) 15.75 - 54.52 7 .3 4. 5 - 3 .6 5. 1

N.4 (ей) 7.62 - 62.41 7.61 - 62.39

Заключение

Для изучения конвекции в малых масштабах и под действием микроускорений применяется, наряду с классической моделью Обербека — Буссинеска, альтернативная модель микроконвекции. Главной отличительной чертой этой модели является учет несоленои-дальности поля скоростей.

Микроконвекция изучалась в длинном прямоугольнике (с соотношением сторон 1 к 10), вытянутом по направлению силы тяжести, в случае теплоизоляции торцов прямоугольника и периодического потока тепла через длинные стороны.

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

Автор выражает искреннюю благодарность А. Ф. Воеводину за консультации при выполнении расчетов.

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

[1] Джозеф Д. Устойчивость движений жидкости. Мир, М., 1981.

[2] Пухначев В. В. Модель конвективного движения при пониженной гравитации. Моделирование в механике, 6(23), №4, 1992.

[3] Андреев В. К., Капцов О. В., Пухначев В. В., Родионов А. А. Применение теоретико-групповых методов в гидродинамике. Наука, Новосибирск, 1994.

[4] Пухначев В. В. Микроконвекция в вертикальном слое. Известия РАН, МЖГ, №5, 1994, 76-84.

[5] Марчук Г. И. Методы расщепления. Наука, М., 1988.

[6] Самарский А. А., Вабищевич П.Н., Матус П.П. Разностные схемы с операторными множителями. Ин-т матем. моделир. РАН, Ин-т матем. НАНБ, Минск, 1998.

[7] Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. Наука, М., 1978.

[8] Роуч П. Вычислительная гидродинамика. Мир, М., 1980.

[9] КАлиткин Н. Н. Численные методы. Наука, М., 1978.

Поступила в редакцию 14 декабря 1999 г.