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

УДК 532.5+519.6

А.Ф. Воеводин, О.Н. Гончарова, О.А. Кондратенко

Численное моделирование переноса тепла в свободном слое жидкости при наличии термокапиллярных сил и дополнительных касательных напряжений

A.F. Voevodin, O.N. Goncharova, О.A. Kondratenko

Numerical Modeling of Heat Transfer in a Free Liquid Layer Under Action of Thermocapillary Forces and Additional Tangential Stresses

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

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

БО! 10.14258^3811(2013)1.2-02

The transfer processes in a plane layer of viscous incompressible fluid with free boundaries are studied under conditions of action of the additional tangential stresses on these boundaries. Modeling of the fluid flows is carried out on the basis of the exact solutions of the Navier- Stokes equations with use of the numerical algorithms to find the velocity field and position of the free boundaries. The problem of finding of temperature distribution of fluid in the infinite layer leads to numerical solution of the heat transfer equation in a rectilinear domain with movable boundaries. The different types of the «smooth» boundary conditions on the «lateral walls» are investigated.

Key words: deformation of plane layer,

thermocapillary forces, free boundaries, tangential stresses, heat transfer, numerical algorithms.

Введение. Моделирование нестационарных течений жидкостей, заполняющих области со свободными поверхностями, предполагает учет дополнительных касательных напряжений на свободной границе и тщательные исследования взаимодействия различных механизмов, вызывающих движение жидкости [1—3]. Из работ по изучению неустановившихся течений плоских слоев со свободными границами следует отметить [4, 5], посвященные математическому моделированию деформации вязкого слоя жидкости термокапиллярными силами (см. также [4-9]). Точные решения, построенные в [4,5,9], представля-

Работа выполнена в рамках проекта №7.3975.2011 Алтайского государственного университета (поддержан Министерством образования и науки РФ) и программы стратегического развития ФГБОУ ВПО «Алтайский государственный университет» на 2012—2016 годы «Развитие Алтайского государственного университета в целях модернизации экономики и социальной сферы Алтайского края и регионов Сибири», мероприятие «Конкурс грантов» (№2013.312.1.66).

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

В данной работе представлена общая схема численного решения полной задачи и численный алгоритм для исследования процесса переноса тепла в прямоугольной области с движущими-

ся «горизонтальными» границами. Данный алгоритм построен на основе продольно-поперечной конечно-разностной схемы, известной как метод переменных направлений. Представлены различные типы «мягких» граничных условий на «вертикальных» границах расчетной области, являющиеся следствием уравнения переноса тепла и условий для температуры на бесконечности. Приведены результаты численного исследования.

1. Постановка задачи о динамике бесконечного слоя вязкой, несжимаемой, теплопроводной жидкости. Пусть жидкость заполняет бесконечный плоскопараллельный слой И = {(х, г): —оо < х < +оо, —Z(t) < г < Z(t)}. Система координат выбирается таким образом, что ось Ох направлена вдоль свободных поверхностей iZ(t), остающихся недеформируемыми и параллельными во все моменты времени, а ось Ог — перпендикулярно к ним. При этом п = (0, ±1) — вектор нормали, в = (1,0) — касательный вектор к свободной поверхности.

Искомые функции (у = (и, го) — вектор скорости, р — давление, Т — температура жидкости) удовлетворяют системе уравнений Навье-Стокса и переноса тепла

Ut

1 ,

Рх Н- t-ч {^ХХ Н- У*Z

tie

wt + uwx + wwz = -pz + —(wxx + wzz),

Re

^X 'Wz 0,

Tt + uTx + wTz = (Txx + Tzz) (1)

KePr

и следующим условиям на свободных границах:

v • п = w\

dZ

dt ’

(2)

-p+ D(v)n\z=±Z(t) = -Pg, (3)

2s • D(v)n|z=±z(t) = t(x, t) - (4)

На свободных границах температура T и касательные напряжения r{x,t) — заданные функции времени и продольных координат:

T(x,±Z(t),t) = \аЦ)х2 + 0(t),

т(х, t) = XT {t).

Здесь A(t), 0(t), r(t) — произвольные функции, зависящие от времени; D(v) - тензор скоростей ______________________ 2

деформации; Рд= Рд - ' D(vE)n\x=±Z(ty,

Рд - внешнее давление; т(х, t) — дополнительные касательные напряжения, вызываемые внешней

средой (т(х,г) = 2Кер1УЪ • £>(у3)п|2=±г(4)). Предполагается, что действиями нормальных напряжений со стороны внешней среды возможно пренебречь. Здесь возникли следующие безразмерные комплексы: Не — число Рейнольдса (Де =

——); Рг — число Прандтля (Рг = —); Ма — число V X

(7 Т I

Марангони (Ма =----------); р,Т/ — отношения плот-

ах

ностей и коэффициентов кинематической вязкости газа и жидкости, где V — коэффициент кинематической вязкости; х — коэффициент температуропроводности; I — характерный размер обла-

I

сти течения, V* — характерная скорость; I* =-----

характерное время; Т* - характерная температура; р* = ру2 — характерное давление; р — плотность жидкости; гид — скорость газа (внешней среды). В случае линейной зависимости коэффициента поверхностного натяжения от температуры имеем: а = ао — <тт(Т — То) (<тт > 0 в случае нормального термокапиллярного эффекта).

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

Решение, определяющее динамику слоя П, имеет вид [2,5]:

и(х, z, t) = fit, z)x, 'jix,z,t) = — / f(t,a)da.

Jo

(6)

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

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

Функции fit,z) и Zit), которые определяют поле скоростей и положение свободной границы, удовлетворяют интегро-дифференциальным уравнениям:

ft + f2 ~ fz ! f(t, a)da - = 0, (7)

<иг гг[г)

dt

f{t,z)dz,

(8)

а также граничным и начальным условиям [9,11]:

fziZ(t),t) =r(t) - -^^-A(t), /z(0,t) = 0, (9)

f(z, 0) = fo{z) = 0 (0 < z < Z0); Z(0) = Zo iZo > 0).

Численные алгоритмы для вычисления /(#, г) и Z(t) являются алгоритмами типа «предиктор-корректор» [9, 12, 13]. Если функция fk(z) = /(ік,г) задана, то положение свободной границы гк+1 можно найти согласно алгоритму «предиктор-корректор» следующего вида:

гк+1 = г*-1 - 2Аі

Jo

предиктор и

Zk+1 =Zk_At

,.zk+1 ,-zk

/ fk+1{z)dz+ / fk(z)dz I о Jo

— корректор [9, 11]. На каждом временном слое +к+1

определяется шаг сетки /?г = Z/M (Z = гк+1). при этом гт = (т — 1) /? -, (т = 1 ,...,М;М = М + 1), = 0, £-ду = Z. Для ап-

проксимации производных функции / используются конечно-разностные аналоги второго порядка. Алгоритм «предиктор-корректор» второго порядка точности используется для решения уравнения (7) с реализацией следующих трех этапов:

/к+1/4 = + о.5Д#[-Л1/й+1/4 + Ф*],

/к+1/2 = /к+1/4 + 0.5Аг[-А2/к+1/2],

/к+1 = /к + АЦ-А/к+1/2 + Фй].

Здесь Л = Лх + Л2, а Л* — конечно-разностные операторы, соответствующие дифференциальным операторам ( —1 /Не)д2/дг2 и — Рд/дг (см. (8) для определения Ф). Интегралы, входящие в (7), (8), вычисляются с помощью квадратурных формул. Тестирование построенного в данной работе численного алгоритма выполняется на основе результатов исследования, представленных в [3,9]. Поскольку на каждом временном слое расчеты осуществляются на новой пространственной сетке (гт), то применяются интерполяционные формулы (см.: [2,9]).

2. Численный алгоритм исследования процесса переноса тепла в прямоугольной области с движущимися границами. Компоненты скорости (см. точные решения вида (6) уравнений Навье-Стокса) находятся численно [10, 11] и зависят только от А(1), 0(#), т(#). Задача сводится к численному нахождению распределения температуры в слое. Для этого введем в рассмотрение прямоугольную область

£1ь = {—1^ < х < Ь, —Z < £ < Z}, (11)

где Z = Z(t). На «длинных» сторонах £ = ±Z прямоугольника О, ь должны быть заданы граничные условия (2)—(4). На искусственно введенных «вертикальных» торцах х = ±Ь данной области должны быть поставлены «мягкие» условия

для температуры, являющиеся следствием условий на бесконечности и уравнения переноса тепла [2,14,15]. При этом ограничимся случаем растекания слоя жидкости (см. рис. 1). Если пред-

Рис. 1. Пример растекания слоя

положить, что на бесконечности (при (х —> ±оо)) поставлено условие Т —> Т,х, где = const, то в предположении Тх —)■ 0 имеем, что на вертикальных торцах может быть поставлено условие вида [2]:

". (12)

Если на бесконечности (при (х —> ±оо)) выполняются все предыдущие предположения, приводящие к условию (12), то на «вертикальных» торцах предполагаем выполнение следующего условия:

Tt + иТх — 0.

(13)

Это условие представляется нам наиболее естественным с физической точки зрения.

Заметим, что на «вертикальных» торцах может быть поставлено также условие вида:

Т( +иТх = \ТХХ. (14)

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

Для численного исследования процесса переноса тепла используется конечно-разностная схема второго порядка аппроксимации. После дискретизации уравнения переноса тепла (1) по времени численная схема записывается в виде:

rpk-\-l/2 _ rj-\k

(15)

0,5Д#

+ХтТ^1/2 - (К2Т)к+V2,

гпк-\-1 _ грк-\-1/2

---------------= ЛтТк+1хх - (К1Т)к+1 +

0,5Д* 1 ;

+ХтТ^1/2 - (К2Т)к+V2.

Здесь Ат = 1/(ДеР?’), Тк(х, г) = Т(1к,х, г).

Конвективный член КТ = иТх-\-гиТ~ представляем здесь в виде суммы двух слагаемых {К{Т = 0.5(иТх + (иТ)) и К2Т = 0.5(гиТг + (м/Г)г).

11\с 11, = {—Ь < х < Ь, —Z < z < Z}

выбрана в качестве расчетной области. Для реализации указанной выше схемы вводится сетка (хп, гт): хп = (п - 1) 1гх (п = 1,..., Й, N = N + 1), 1гх — шаг в направлении оси Ох, 1гх = Ь/М. Данная сетка является подвижной по вертикальной

координате: гт = (то — 1)Ьг (то = 1, ...,М,М = М + 1), где Ъ,х — шаг в направлении оси Ог, ]гг = Z/M, а Z = Zk+1 — новое положение свободной границы на временном слое (к + 1). При построении решения в половине расчетной области, симметричной относительно Оу, должно быть поставлено условие при х = 0 следующего вида: Тх= 0.

Пусть Ткт = Т(гк,хп,гт). Тогда уравнения

(15) представляют собой системы линейных алгебраических уравнений, которые можно записать в виде

— а грк+1/2 т грк+1/2 _ грк+1/2 _ ,

и'П,т-і-п..т.— 1 ' ип,тл- п,т ^п,т п,т+1 ип,ті

—а ггк+1

7 Грк+1 _ ГрК^І _ 1

ип,гп-1-п,т (-п,т-/-п+ 1?т

пк+1

Т,

п+ 1,т

- т

-I- п.

1ъ1

__ Тп тп Тп—і т ип-^іТлЦ-і уп ^п — іТп—І т

КІ

2 К

3. Результаты численного исследования. Тестирование численного алгоритма проводится с использованием (формального) точного решения (ТР) уравнения теплопроводности вида:

Т

-Хі

(16)

и решить методом прогонки [16]. Для аппроксимации вторых производных применяются разностные аналоги второго порядка и центральные разности для аппроксимации конвективных слагаемых. В качестве примера приведем разностную аппроксимацию слагаемого иТх + (иТ)х:

Здесь ип+1 = 0.5(м„+1 + ип).

Первые производные внутри расчетной области аппроксимируются традиционно симметричными разностными аналогами со вторым порядком. Первые производные на границе расчетной области будут аппроксимироваться несимметричными конечно-разностными аналогами также второго порядка. Так как пространственная сетка подвижна в направлении Ог, а переход на новый временной слой (/г + 1) начинается с расчета новой пространственной сетки, то следует найти значение температуры Тк, рассчитанной на предыдущем слое, в узлах новой сетки (гт). Мы используем интерполяционные процедуры. С помощью интерполяционных формул Ньютона [12]

гм = ТЫ + (т(г')~г(го)) (*-*„)+

+ (ты-2ты + ты)(г-гоКг-г1)+

Т(г3) - ЗТ(г2) + ЗТ(г1) - Т(г0)

6/г|

■(г - г0)(г - - г2)

вычисляем на новой пространственной сетке значения искомых функций, известных на предыдущем временном слое.

Здесь Л = 2/(ДеРг) — коэффициент температуропроводности. Пусть Де = 1, Рг = 17. Численное решение находится в случае, когда положение границ фиксировано и определяется значением Z =

0.5. Для численного исследования процесса теплопроводности используется конечно-разностная схема второго порядка аппроксимации вида (15). Пусть Ь = 10 (см. (11) для определения области О-ь, выбираемой в качестве расчетной).

При проведении тестовых расчетов в расчетной области О-ь с шагом Нх = = 0.002 по про-

странственным переменным различие между аналитическим (16) и численным решениями наблюдается только в седьмом знаке после запятой, что демонстрирует хорошее совпадение численного и точного решений (г = 0.0001).

При проведении тестовых расчетов с шагом = 0.02 по пространственным переменным различие между аналитическим (16) и численным решениями наблюдается в шестом знаке после запятой (г = 0.0001). Отметим, что при увеличении шага по пространству результаты становятся менее точными. Тестовые расчеты с помощью точного решения (16) проводились в случае, когда на границе области заданы значения самой функции. Вместе с тем проведены тестовые расчеты и в случае, когда на торцах заданы условия типа (12). На основе точного решения неоднородного уравнения переноса тепла с правой частью, продиктованной видом полного уравнения переноса тепла (см. (1)) и точного решения (16), проведены тестовые расчеты согласно полной схеме (15), которые также продемонстрировали хорошее совпадение результатов. Различия в численном и точном решениях в случае задания на всех границах области «температуры» согласно (16) наблюдаются в третьем знаке после запятой, и в случае задания на «вертикальных» торцах условия (12) также в третьем знаке после запятой. Расчеты выполнены при ]гх = ]гг = 0.002, г = 0.0001.

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

В данной работе проведено сравнение результатов о распределении температуры для различных типов мягких граничных условий на торцах (см. (12)—(14)). На рисунках 2—7 показано распределение температуры в слое жидкости в случае различных зависимостей от времени функций А{Ь) и г(£) (величины представлены в безразмер-

ном виде). Исследования проводились в случае растекания слоя (см. рис. 1). В начальный момент времени положение свободной границы жидкости определяется как Zq = 0.5. Для зависимости A{t) = {А0 = -0.1) и r{t) = Ату^-р

[Ат = 400, го = —0.1) приведено сравнение результатов численного построения решения для момента времени t = 0.05 (рис. 2—4). Для всех типов условий в результате наблюдается качественно одинаковая картина, но имеют место некоторые количественные различия. Эти отличные для разных типов граничных условий особенности проявляются вблизи «вертикальных» торцов. Это может быть связано с аппроксимацией первых и вторых производных для температуры. Наибольшее значение на торце температура имеет при использовании условия (14), наименьшее — при использовании (12).

А

Рис. 2. Случай A(t) = ------ (Ао = —0.1), т(£) =

1 +

Агт^т (то = -0.1, = 400), Тхх = 0,to = 0.05

1 + tA

An

Рис. 3. Случай A(t) = ------------ (Ао = —0.1), r(t) =

1 -\~ t

Ат—^— (то = -0.1, Аг = 400), Tt+uTx = 0, t0 = 0.05 1 + tA

Ап

Рис. 4. Случай А{1) = -------- (Ао = —0.1), т(£) =

1 ~\~ £

Ат Т° (т0 = -0.1, Ат = 400), Г4 + иТх = \ТХХ, to = 0.05

Аналогичные результаты получены и в случае экспоненциальных зависимостей А(1) =

А01ехр(1) (Лз = 0.1), (Лз = -0.1), т(1) = АтТ()1 ехр(1) (Ат = 400, то = 0.1) при 1 = 0.05 (рис. 5-7).

Рис. 5. Случай A(t) = Aotexp(t) (Ао = —0.1), т(£) = ATT0t exp{t) (то = 0.1, Ат = 400), Txx = 0, to = 0.05

Рис. 6. Случай A(t) = Aotexp(t) (Ao = —0.1), ?(f) = Atto texp(t) (to = 0.1, AT = 400), Tt + uTx = 0, to = 0.05

0 0.6 1 1.6 2 25 3 3.5 4 4.5

Рис. 7. Случай А(£) = Aotexp(t) (Ао = —0.1), т(£) = Aттotexp(t) (то = 0.1, Ат = 400), Г( + иТх = \ТХХ, to = 0.05

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

Библиографический список

1. Гончарова О.Н., Кабов О.А. Гравитационно-термокапиллярная конвекция жидкости в горизонтальном слое при спутном потоке газа // ДАН. - 2009. - Т. 426, №2.

2. Goncharova O.N., Kabov О.A. Deformation of a viscous heat conducting free liquid layer by the thermocapillary forces and tangential stresses: Analytical and numerical modeling // Microgravity sci. technol. — 2010. — Vol. 22, №3.

3. Goncharova O.N., Kabova Yu.O., Kabov O.A. Thermocapillary convection in a free liquid layer in the presence of an adjacent gas flow // J. Computational Thermal Sci. — 2011. — №3(5).

4. Pukhnachov V.V. On a problem of a viscous strip deformation with a free boundary C.R. // Acad. Scien. — Paris, 1999. — №328(1).

5. Pukhnachov V.V. Model of a viscous layer

deformation by the thermocapillary forces // Max-Planck-Institut fuer die Mathematik in den Naturwissenschaften: preprint. — Leipzig, 2000. —

№50.

6. Longuet-Higgins M.S. A class of exact, time depemdent free-surfece flows // J. Fluid Mech. — 1972. - №55 (3).

7. Андреев В.К. Устойчивость неустановив-шихся движений жидкости со свободной границей. — Новосибирск, 1992.

8. Galaktionov V.A., Vazquez J.L. Blow-up of a cass of solutions with free boundary for the

Navier-Stokes equations / / Advances in Differential equations. — 1999. — V. 1.

9. Пухначева Т.П. Численное решение задачи о деформировании вязкого слоя термокапиллярными силами // Симметрия и дифференциальные уравнения. — Красноярск, 2000.

10. Гончарова О.Н., Кондратенко О.А. Деформация вязкого теплопроводного слоя в условиях дополнительных касательных напряжений / / Известия АлтГУ. — Барнаул, 2011. — №1/2(69).

11. Гончарова О.Н., Кондратенко О.А. Численное моделирование переноса тепла в свободном слое жидкости на основе точных решений уравнений Навье-Стокса // Известия АлтГУ. — Барнаул, 2013. - №1/1(77).

12. Самарский А.А., Гулин А.В. Численные методы. — М., 1989.

13. Марчук Г.И. Методы вычислительной математики. — М., 1977.

14. Воеводин А.Ф., Овчарова А.С. Численное решение задачи о качестве воды в открытом русловом потоке // Водные ресурсы. — 1977. — №4.

15. Роуч П. Вычислительная гидродинамика.

- М., 1980.

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