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

МИКРОКОНВЕКЦИЯ В ОБЛАСТИ СО СВОБОДНОЙ ГРАНИЦЕЙ *

О. Н. Гончарова Новосибирский государственный университет, Россия Институт гидродинамики им. М. А. Лаврентьева СО РАН Новосибирск, Россия e-mail: [email protected]

The steady gravitational-thermocapillary convection in semicircular regions with free surfaces is considered. Cases when both mechanisms of the convection (gravity and thermocapillary) play an essential role are analyzed. Two models of the convection are used: the first one is the classical Oberbeck — Boussinesq model and the second model is the new model for description of the convection at low gravity forces, in small domains and fast changing temperature fields. The velocity field in the new model is not a solenoidal vector.

This paper presents the results of calculations of convective flows in the domains with free boundaries at the microgravity forces typical of a space orbital laboratory. Numerical solutions are obtained for a range of Prandtl number and at different Marangoni and Grashof numbers. Velocity and temperature fields are presented. Analysis of conditions of applicability of both models is given.

Термин “микроконвекция” появился в статье [1] для характеристики течений с малыми значениями параметра n = gl3/vx, где l — диаметр области, занятой жидкостью, v, х — коэффициенты кинематической вязкости и температуропроводности, g — ускорение силы тяжести. Еще один специфический термин — “изотермически несжимаемая жидкость” — впервые введен, видимо, в работе [2] для определения жидкости, плотность которой р зависит лишь от температуры Т. (Предположение о независимости плотности от давления P традиционно для теории конвекции; см., например, [3].)

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

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

1. Жидкость изотермически несжимаема и р = ро(1-вТ) (в — коэффициент объемного расширения).

2. Движение подобно движению несжимаемой жидкости, поэтому поле скоростей счи-

* Работа выполнена при поддержке Фонда Гумбольдта (Германия) и Российского фонда фундаментальных исследований, грант №99-01-00529.

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

14

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

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

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

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

Вместе с тем уравнение энергии упрощается с учетом гипотезы (3), и для простоты считается выполненной гипотеза (4). Гипотеза же (1) принимается в виде р = р(Т). Новая модель характеризуется несоленоидальностью поля скоростей. При использовании зависимости плотности от температуры вида р = р0/(1 + вТ) полученная система уравнений преобразуется к виду, в котором модифицированный вектор скорости Ш = V — вх^Т уже становится соленоидальным. Соленоидальность модифицированной скорости позволяет ввести аналог функции тока для плоских и осесимметричных задач и выполнять расчеты конвективных течений в переменных “функция тока — вихрь”.

Отметим, что с физической точки зрения обе зависимости плотности от температуры р = ро/(1 + вТ) и р = р0(1 — вТ) практически равносильны (в реальных конвективных течениях максимальные значения величины в|Т| не превышают 10-2).

Параметр п является дополнительным по отношению к числам Грасгофа Ог и Пран-дтля Рг критерием подобия, роль которого существенна для моделирования конвекции в слабых силовых полях, в микромасштабах, а также в жидкостях с большим произведением коэффициентов вязкости и температуропроводности. Как отмечено в [4], приближение Обербека — Буссинеска непригодно для описания конвекции, если п < 1. С физической точки зрения этот параметр равен отношению порядков скоростей, порожденных объемным расширением жидкости (вТ*х/1) и фактором плавучести (вТ*\2д/у).

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

В работе [7] проведено исследование конвективных процессов в областях со свободными границами, когда уже оба механизма конвекции играют существенную роль. Течение предполагается стационарным. Рассчитывается гравитационно-капиллярная конвекция в кольцевых областях. Одна из границ по-прежнему считается твердой, а другая является свободной, подверженной действию термокапиллярного эффекта. Расчеты проведены для различных значений чисел Прандтля, Марангони и Грасгофа. Исследована структура конвективных течений; предложены те ситуации, когда топология течения может существенно различаться. Количественные различия в величинах скоростей здесь также значительно меньше, чем выявленные в работе [6]. Результаты расчетов согласуются с аналитическим исследованием стационарной задачи микроконвекции [8], показавшим пригодность и классической модели для описания таких течений.

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

В работе приняты следующие дополнительные обозначения: г — радиальная координата; р — угловая координата; ш — вихрь скорости;

ф — функция тока, модифицированная функция тока;

и = —------тангенциальная компонента скорости;

дг

д — модифицированное давление;

Т* — относительная температура;

Т0 — характеристическая температура;

Я — радиус полукруга;

р0 — характеристическая плотность;

о = о0 — 7(Т — Т*) — коэффициент поверхностного натяжения;

7 — температурный коэффициент поверхностного натяжения.

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

Стационарная гравитационно-капиллярная конвекция исследуется в полукруге

Свободной границей является диаметр полукруга (р = п, р = 2п, 0 < г < Я), а полуокружность (г = Я, п < р < 2п) представляет собой твердую границу с потоком тепла через нее.

Уравнения конвекции рассматриваются в переменных “функция тока — вихрь” в безразмерной форме. При этом характерный размер, скорость, температура и давление выбираются, согласно [9], равными соответственно Я, 7Т0/^, Т0, и 7Т0/Я.

2. Классическая модель Обербека — Буссинеска

Уравнения рассматриваются в полярной системе координат:

1 дф

радиальная компонента скорости;

0 < г < Я < +то, п < р < 2п.

(1)

Дф + ш = 0,

(2)

і дТ пдТ\ п АТ — Ма ( V——|--------— | — 0,

дг г др /

(3)

где Ив — —, Ма — ИвРг, Рг

V

, Ог

9вТ0—3

числа Рейнольдса, Марангони,

V” X Vх

Прандтля и Грасгофа соответственно; произведение вТ0 можно представить следующим

Ог дК3

образом: вТ0 = —, где ц = ------ .

П "X

Для того, чтобы при рассмотрении стационарной конвекции реализовать быстропеременные температурные поля, создаются локальные особенности теплового потока как через свободную, так и через твердую границу. В связи с этим рассматриваются два вида граничных условий для обеих моделей. Обозначим их условно: I — “всплеск” на свободной границе, II — “всплеск” на твердой границе. (Для сравнения будет рассмотрен и базисный вариант, когда “всплески” отсутствуют.)

Вариант I. Граничные условия могут быть представлены в виде 0 < г < К, р = п, р = 2п :

ф — 0, и

ТГ, р — 2п,

дг ’

дТ

~дr,

р — п,

дТ

др

0, р — п, р — 2п (г — — *), —Тв, р — 2п (г — — *).

На твердой границе г — — :

дф дТ

Ф = 0, — = 0, — = Тс сое р.

дг дг

Вариант II. Граничные условия могут быть представлены в виде 0 < г < К, р = п, р = 2п :

ф — 0, и

дТ

~дr,

дТ

~дт

р — 2п, р — п,

дТ

др

0.

На твердой границе г — — :

ф — 0,

дф

дг

дТ

дг

Тс сое р, р — р*,

Тс сое р, р — р„.

3. Модель микроконвекции

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

|1 + (вТ„)Т! Дш - Ие ^) +

0

+(в-о)< ~Тт — +

г др дг г дг др

дТ ( п — (Ап —

дг

1 дТ / V — — (Аv —

г др

|

Ог /дТ

1 дТ

Ма V д~ С08р — 'гір 81Пр

вТо ( АТ + дТди + 1 дТди +

( и АТ |—-—-—|—^ —| А(и,Т,і,п) I —

Рг

дг дг г2 др др

в 2То2 Ма Рг

1

дТ дАТ дТ дАТ

+

г дг др др дг

Аф + и — 0,

(4)

(5)

Здесь

[1 + (РТо)Т] АТ — Ма. ( V3- + — ) - (вТо) | V Т|2 — 0.

Группа слагаемых

Л 1 д д 1 д2 А — -—г— +

г дг дг г2 др2

V — | —, -1 —

дг г др

(6)

„ . „ д2Т Зі,? , д2Тдь2 д2Тдь,

А(и, і ,і,п) — 2 ^ ^ —----+

Зхду ду дх2 дх ду2 ду также записывается в полярных координатах, при этом і\, і2 — компоненты вектора скорости в декартовой системе координат Оху.

Вариант I. Граничные условия имеют вид 0 < г < —, р — п, р — 2п:

ф — 0, и

ЗЗг, р — 2п,

дг дТ

ддг, р — п,

дТ Г 0, р — п, р — 2п (г — — *),

3р 1 —*Тв, р — 2п (г — — *).

На твердой границе г — — :

, пв-от . дф 1 вТо зт зт

ф — —— мі-с 81,1 р аг — — мізр- зг — То С08 р

Вариант II. Граничные условия имеют вид 0 < г < —, р — п, р — 2п:

дТ

-з-, р — 2п, зт

ф — 0, и — < ^ — — 0.

^ р = п.

дг

+

0

На твердой границе г = Я :

ф

_Яш

-Я

РТ0 гр Ма Тс

дТ

дг

р, Р = р*, дф

дг

р, Р = р*,

Та сое р, р = р*

% сое р, р = р*

дФ = I вТ0 дТ

дг Я Ма др ’

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

Численное исследование поставленных задач для систем уравнений (1) - (3) и (4) - (6) проводится методом установления с использованием продольно-поперечной конечно-разностной схемы, формально имеющей второй порядок аппроксимации [10]. В связи с тем что в уравнениях (1), (3), (4), (6) конвективные слагаемые берутся с нижнего “временного” слоя с использованием направленных разностей, порядок аппроксимации становится первым. Предлагаемая для исследования методика была апробирована на тестовых задачах и расчетах нестационарной свободной конвекции в кольцевой области [11].

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

цк+1/2 _ цк

Цк+1 _ цк+1/2

= Ли [Л^к + л2ик+1/2 + тк], = ~Ли [Л1Цк+1 + Л2Цк+1/2 + тк]

(7)

где и

0.5т

ик = и(¿к), Л1 и Л2 — разностные операторы для аппроксимации

соответствующих дифференциальных, Ли = Хц для модели Обербека — Буссинеска, Лц = Хц(1 + (вТо)Тк) для модели микроконвекции, Хц — итерационный параметр, а Тк включает в себя все, начиная со второго, слагаемые в левых частях уравнений (1), (3), (4), (6), насчитываемые на предыдущем слое.

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

ф

«+1/2 _

фЯ

0.5т

ф8+1 — ф^+1/2

0.5т

= Х^ (Л1ф*+1/2 + Л2ф* + ик+1)

Хф (Л1ф5+1/2 + Л2ф3+1 + ик+^,

с параметром Х^.

Для реализации вышеизложенной схемы расчета вводится разностная сетка гп = (п _ 1)Н, (п = 1,...,М + 1), Н = Я/М,

Рт = (т _ 1)а, (т = т,...,М + 1), а = 2п/М, (та = п), /(гп,рт) = /п,т.

При этом имеет место следующее представление для Лl¡ и Л2/:

Лі/ = -

тп

тп+1/2

/n+1,m /

h2

n,m /n,m /n-1,m

— Tn—1/2 ' _

h2

/

n,m+1

-2/n ,m + ,fn ,m-1

т2 a2

n

Для аппроксимации конвективных слагаемых используется идея аппроксимации против потока:

v/ + дт т др

/n+1 ,m /n-1,m + un,m

/

n,m+1

/n,m-1

2h ' Tn 2a

/n+1 ,m - 2/n ,m + /n- 1,m ^n^1 /n,m+1 - 2/n ,m + in ,m-1

+

Смешанные производные типа

2h

д 2/

дтдр

+

2a

входящие в А(ш,Т, v,u), и первые производные вну-

три расчетной области традиционно аппроксимируются симметричными конечно-разностными аналогами со вторым порядком точности [12]. Первые производные на границе расчетной области аппроксимируются односторонними разностями.

Для задания граничного условия для вихря на твердой границе выводятся условия типа Тома [13]:

+1,т = —То Фм,т,

П,2

2 . вТо дТ

ШN + 1,m — — Т2 WN,m —

h2

Ma др

1 + 2 \ вТо Т . Л+2^

r2 + ш)- маТоsinрпr + к2)'

для классической модели и модели микроконвекции соответственно. При этом в случае II вместо Tg следует писать Tg для р = р*.

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

1. Внешний итерационный процесс состоит в последовательном расчете функций Тk+1, uk+l из систем уравнений (7). При этом на промежуточном (k+1/2)-м слое осуществляется прогонка в направлении р, а на основном (k + 1)-м слое — прогонка в направлении r (стартовые значения определяются состоянием покоя: Т := Т0 = const, ш := 0, ф := 0).

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

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

[13]:

max - fn 1 < Ef max Ijn+mil,

n,m n,m

где i — номер итерации, Ef — заданная точность расчета fг+1.

Помимо того что осуществляется дополнительная проверка точности выполнимости граничных условий по величине [13]

є — max |wN+1i,m(^Ntm) - ^N+W^N.m^

v

n,m

т

n

стационарное течение считается достигнутым, если выполняется не менее K внешних итераций и max k — fnm1 < £f •

n,m

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

max ^П^1 [13]).

n,m

Вопросы, связанные с поправкой свободных границ и устойчивостью течений, можно в случае необходимости решать согласно [9, 14]. Если при этом H(x) — отклонение свободной границы от положения у = 0, — R < x < R, то уравнение для этой поправки может быть записано в следующем виде:

о 1 л R

5P — 2 = — H", H (±R) = 0, Hdx = 0.

ау Са J-R

Здесь штрих означает производную по x, Са = ^Т0/о0 , v2 выражается через радиальную и тангенциальную компоненты скорости v, и, а 8P представляет собой отклонение давления от равновесного уровня.

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

Для жидкостей типа глицерина (Glyc.1, Glyc.2, Glyc.3), расплавов кремния (Silic.1, Silic.2) и стекла (Glass 1, Glass 2) при действии микроускорений, достижимых на орбитальной станции, проводятся расчеты поставленных задач с использованием как традиционной модели, так и модели микроконвекции.

Расчеты выполняются на сетках 41 х 41 и 81 х 81. Радиус полукруга R = 1 см. Значения параметров в граничных условиях для варианта I: Тв = 70, 700, Tq = 35; для варианта II: Tg = 35, Tq = A ■ Tq, A = 10.

В таблице приведены безразмерные параметры.

Параметры задачи

Вещество Pr Ma Re Gr п вТо Gr Ma вТо Pr в2Т2 Ma Pr

Glyc.1 104 3 ■ 102 2 - 0 1 3 .1 5 1 О 1 со 10-1 1.5 ■ 10-2 СО - 0 1 5 1.5 ■ 10-6 10-10

Glyc.2 104 1 10-4 .0 5 1 О 1 Сл 10-1 .0 5 1 О 1 СО - 0 1 5 .0 5 1 О 1 00 10-12

Glyc.3 104 1 10-4 .1 5 1 О 1 со 10-1 .1 5 1 О 1 ю .1 5 1 О 1 со 6 - 0 1 .5 1 10-8

Glass 1 104 10 10-3 со - 0 1 .5 .0 10-2 .4 5 1 О 1 Сл 5 1 О 1 00 4.5 ■ 10-9 10-14

Glass 2 104 1 10-4 .0 5 1 О 1 10-2 СО - 0 1 .5 .4 5 1 О 1 00 4.5 ■ 10-10 10-15

Silic.1 4 1 о 1 со 1 2 0 1 Ю .0 3 1 О 1 1 .0 3 1 О 1 .0 3 1 О 1 .0 8 1 О 1 Сл 10-13

Silic.2 4 1 О 1 со 1 2.5 ■ 102 2 1 о 1 1 2 1 О 1 2 1 О 1 5 1 О 1 ю 10-5

Расчеты, проведенные по двум альтернативным моделям в случае отсутствия особенности теплового потока на границах, показывают лишь некоторые количественные различия и практически одинаковую качественную картину. Про этом стационарное решение имеет одновихревую структуру (см., например, рис. 1, а), а для жидкости 01уе.1 имеет место двухвихревое течение типа “два маленьких вихря в одном”.

При расчетах с особенностью теплового потока на твердой границе качественные различия в топологии течения уже наблюдаются для жидкостей 01уе.1, 01уе.3. При этом для первой жидкости проявляется ярче двухвихревая (типа “два в одном”) структура течения, когда расчеты ведутся по модели микроконвекции (рис. 2). Для второй жидкости одновихревая структура при расчетах по модели микроконвекции сохраняется, но центр вихря смещен вправо вниз к точке всплеска (ср. рис. 1, а, б). Для остальных жидкостей получаются, по-прежнему, одновихревые структуры. Следует заметить лишь сдвиг центра влево для кремния при использовании модели микроконвекции (рис. 3).

Рис. 1.

Рис. 2.

1.00 - 0.80 -0.60 -0.40 - 0.20 0.00 0.20 0.40 0.60 0.80 1.00

Рис. 3.

■loo Н--------1-------1-------1-----1 ~ - -— Г-1-1-1--

-1.00 - 0.80 -0.60 - 0.40 - 0.20 0.00 0.20 0.40 0.60 0.80 1.00

Рис. 4.

При расчетах с особенностью теплового потока на свободной границе можно говорить либо об одновихревом (как на рис. 1, а), либо о двухвихревом (рис. 4) течении для жидкостей типа расплава стекла (Glass 1, Glass 2), глицерина (Glyc.2, Glyc.1), кремния (Silic.1). При этом распавшаяся двухвихревая структура наблюдается при более интенсивной особенности (а для Glyc.1 при более слабой особенности, как и в предыдущих случаях, наблюдаются два малых вихря в одном). Отличия же при расчетах по разным моделям проявляются при исследовании течения жидкостей Glyc.3 и Silic.2. Модель микроконвекции позволяет говорить о двухвихревой распавшейся структуре со смещением второго вихря к свободной границе, к зоне всплеска (см. рис. 4), в то время как течение, рассчитанное по классической модели, имеет одновихревую структуру (как на рис. 1, а).

Порядки безразмерных скоростей для кремния, глицерина и стекла равны соответственно на свободной границе: ~ 100 — 101; ~ 10-1 — 101; ~ 100 — 101 ; внутри области: 100; ~ 10-1 — 101; ~ 10°. Размерные скорости (см/с) — порядка ~ 10-2 — 10° (кремний Silic.1, Silic.2); ~ 10-3 — 10-1 и ~ 10-4 — 10-2 (стекло Glass 1, Glass 2 соответственно); ~ 10-2 — 10° и ~ 10-4 — 10-2 (глицерин Glyc.1 и Glyc.2, Glyc.3 соответственно).

На рис. 5-7 приводятся типичные для рассмотренных задач семейства изотерм. Изотермы имеют характерный для конвекции вид. При исследовании по двум альтернативным моделям качественных либо количественных различий в температурных полях не наблюдается. Вместе с тем очевидным образом проявляются отличия при рассмотрении вариантов I, II (см. рис. 6, 7) и базового (в случае отсутствия на границах особенностей теплового потока, см. рис. 5).

Следует заметить, что, как и в работе [7], качественные отличия между расчетами, про-

Рис. 5.

веденными по двум моделям, возникают при таком изменении характеристик свободной поверхности, когда, например, 7 ~ 10-4 г/(е2-К), а Са ~ 10-5.

Рис. 6.

Рис. 7.

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

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

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

[2] НАДОЛИН К. А. Конвекция в горизонтальном слое жидкости при инверсии удельного объема. Там же, №1, 1989, 43-49.

[3] ДЖОЗЕФ Д. Устойчивость движений жидкости. Мир, М., 1981.

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

[5] ГОНЧАРОВА O. Н. Численное исследование уравнений конвекции изотермически несжимаемой жидкости. В “Вычисл. технологии”, 2, №5, 1993, 95-101.

[6] ГОНЧАРОВА O. Н. Микроконвекция в слабых силовых полях. Сравнение двух моделей при численном исследовании. ПМТФ, 38, №2, 1997, 58-63.

[7] ГОНЧАРОВА O. Н. Численное исследование микроконвекции в областях со свободными границами. Там же, 38, №3, 1997, 64-68.

[8] ПУХНАЧЕВ В. В. Стационарная задача микроконвекции. Динамика сплошной среды, 111, 1996, 109-116.

[9] ZEBIB A., HOMSKY G. M., MEIBURG E. High Marangoni number convection in a square cavity. Phys. Fluids, 28, No. 12, 1985, 3467-3476.

[10] САМАРСКИЙ А. А., ВАБИЩЕВИЧ П.Н., Матус П.П. Разностные схемы с операторными множителями. Ин-т матем. моделир. РАН, Ин-т матем. НАНБ, Минск, 1998.

[11] ВОЕВОДИН А. Ф., ГОНЧАРОВА О. Н. Расчет свободной конвекции при изменяющемся поле сил тяжести. Динамика сплошной среды, 67, 1984, 21-28.

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

[13] ТАРУНИН Е. Л. Вычислительный эксперимент в задачах свободной конвекции. Изд-во Иркутского гос. ун-та, 1990.

[14] БАБСКИЙ В. Г., КОПАЧЕВСКИЙ Н. Д., МыШКИС А. Д., СЛОБОЖАНИН Л. А., Тю-ПЦОВ А. Д. Гидромеханика невесомости. Наука, М., 1976.

Поступила в редакцию 22 января 1999 г., в переработанном виде 30 апреля 1999 г.