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

А.П. Васильев

ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ СИСТЕМЫ УРАВНЕНИЙ ДИНАМИКИ И ТЕПЛООБМЕНА ПРИ ПУЛЬСАЦИЯХ ГАЗОВОГО ПУЗЫРЬКА

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

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

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

_Э Эг 1 Э

(р>1 )=о,

1Г+7 Э

—1+—1 Э 1 Эг

эр2

Э1

1 Эр, 2 а-^Б2

Р0 Эг 3 р0

+у

)=о,

Э2-^ 2 Э—

+ 2 и— 2—1

р1 С]

ЭТ. ЭТ. —1+—1

Э 1 Эг

1Э

г2 Эг

Эг

Эг2 г Эг

— 2 то

+12-1 +—а-^Б2

г2 3 1

р2с2

ЭТ2 Э1

+—

Р2 () = *2р0Т2 ЭТ2 4 1 э/ Эг

1 _Э_

•2 Эг

2 ЭТ2 12г^

Л

жидкой фазы, В - индукция магнитного поля, Я2 -удельная газовая постоянная.

Для численного исследования выписанной системы уравнений перейдем к новым безразмерным переменным г * = гД0 - время, - = г/а (г) -радиус, где а(() - текущий радиус пузырька, а г0 - характерное время задачи, причем [2]:

а0 Э 1 Э Э/ч Э/ч а Э/ч ^/р0 Эг а Э- Э1 Э1 11 а Э-

и преобразуем ее к безразмерному виду.

Уравнение динамики пузырька

ё2а * 3 2 2

ёа * Л.

2

4 2 На2 ёа*

+--а *—* +

3 Яе1 Л,

4 1 ёа*

2 1

Яе а * & * We1 а *

= Р *-(к +1),

где к - параметр возмущения.

Уравнение энергии жидкой фазы

ЭЭ1 + а* Э * а *

—-

Э91 = 1 1 Э- = Ре1 а 2

+12

Ес1 Яе1

У

/. \2 а *

а*

'Э^ 2 Э91 4 -1 +---1

Э"2 - Э-

2 На2Ес1 а2 -6 3 Яе1 -4'

Уравнение для давления газа в пузырьке

Ф*

= 3

1 Эе,

- 3ур*

-=1

&* Ре2Ес2 а; Э-

Уравнение для скорости в газовой фазе

р* у—1 1 Эе 2 1 ф

—

,(-д *)=

Ре2Ес2 у р*а* Э- 3у &*

-'

Уравнение энергии для газовой фазы

Здесь индекс «1» относится к жидкой фазе, а «2»- к газовой; р0 - истинная плотность, с1 -удельная массовая теплоемкость, 11 - теплопроводность, —1 - скорость, р1 - давление, —1 - скорость, Т1 - температура 1-ой фазы, V, т - кинематический и динамический коэффициенты вязкости

Эе Эг

2 + -2*— -а* Эе2

Э-

у-1 р*

Э2е

у Ре2Ес2

Э-

2+ 2 Эе^ 2 - Э-

у—1 е2 ф* у р* &*

Система уравнений (1)-(5) является замкнутой. Неизвестными в ней являются: а * (г *) - при-

а

1

+

2

-

г

*

а

веденный радиус пузырька, р. (1 *) - приведенное давление, -№2*(т|,1 *) - приведенная скорость газа, 01 (т,1 *), 02 (т,1 *)- приведенные температуры в фазах, причем все неизвестные функции определены так:

а * =

а

р* =

^ * = -

w

)

Ро

0 = Т1

1 То , 2 То ,

где и0 = а0/10 - характерная скорость, а все другие величины с индексом «0» относятся к начальным параметрам, - - показатель адиабаты.

Краевые условия для системы уравнений (1)-(5) задаются в виде

1,= 0: а* = 1, а * = w, = 0, 81

Т = 0: ^ = 0, 82 <<

Э8 2

Эт

Э81

Т = 1: 12 ^ = 1, 82

т = ¥ : 81

Критерии и числа задачи определены так:

На = а0В _ - число Гартмана,

Vт

Яе1 =

а0 и0

- число Рейнольдса,

р0 и2

^е1 = ^ / - число Вебера,

Ч а0

где Ч - коэффициент поверхностного натяжения,

Ре1 =

а0^ х, Ре2 =:

а0 и0

- числа Пекле,

Д2а,

Д1

Др. Д1,

= Ба (а,,а,

= (а ,,а , ,р , ,qД

где через а* обозначена производная q * =

Заменим производные по временной и пространственной переменным конечными разностями вида:

#

Э1

51

Эf

Эт

5т

Р20 , „0

р* =—— - приведенная плотность фаз, р20 - истин-

Р0

ная плотность газовой фазы в начальный момент времени.

Обыкновенные дифференциальные уравнения системы (1) и (3) в дальнейшем удобнее представить в виде

Э ^ (т,1) = f (т + 5т,1)- 2f (т,1)+f (т-5т,1) (13) Эт2 (5т)2

и введем сетку: т = т5т, 1 = п51, т = 0,1,...,М; п = 0,1,...,К, 5т = Ь/М, 51 = т/К - шаги по пространственной и временной переменным.

Уравнения (10) и (11) могут быть проинтегрированы по схеме Рунге-Кутта, например, четвертого порядка точности [3]. В этом случае для них можно написать следующие сеточные уравнения

(w * = а *):

ап+1 = ап + Wn51 +1 (( + к2 + кз )51, (14) 6

Wn+l = Wn +1

6

Рп+1 = Рп +1 (т1 + 2т2 + 2тз + т4), (16) 6

где коэффициенты Ц и ш1 (1=1,2,...,4) находятся через функции Ра (а* ,а* ,р*) и Рр (а* ,а* ,р* ,q*) согласно выбранной схеме [3], например,

к1 = Ра ОЧп^^Рп )51, т1 = Рр.^^Рп^п Ж" . Уравнение (4) на сетке ш5т, п51; записывается так

w & —- ^ Бр

Ре2Бе2 - Рпап 3-

а уравнения энергии для газовой и жидкой фаз в виде

8п+1,т = 8п,т + 51

wп2)- m5тwn 8п:т+1 -8 ап 5т

3

п,т

+--рр (an,wn,Рn,qn)+

- рп р п п п п

--1 р. 8п

(2) т

- Ре2БС2 Рпап

Х8(2) - 98(2) +8(2) 4

п,т+1 п,т п,т-1

(5т)2

2 8(2) - 8(2) 2 °п,т+1 °п,т

т5т 5т

т=1

V

+

+

(1) _ 0(1) п,ш

9п+1,ш _ 9

+ 51

(ш5л)2

- ш5л

9п,т+1 - 9

(1) п,ш

5л

?е1аП

9П1,)т+1 - 29

(1) п,ш

+ е:

(5л)2

еП1,)ш+1- е

(1) п,ш

+12

ш5л

Ес1 Яе1

5л

\2

(ш5л)6

2 Иа2Ес,

+1

3 Яе1 (ш5л)4

Начальные условия 1=0, п=0 для пузырька записываются так:

а0 _ 1, о _ о, р0 _1 q0 _ 0.

Начальное распределение температуры в газовой фазе

е02ш_ 1, ш _ 0,1,2,...,Мв, Мв _ 1/5л,.

Начальное распределение температуры в жидкой фазе

3(1) _ ]

0,ш

,М.

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

а) в центре пузырька

ш _ 0: е<2)_е<2);

б) на межфазной поверхности

е

qn _

е (2> -е1

п,МО п,МО-1

е(2) _е(1) е(1) _е(1) +^а 5л •

п,МО ип,Мв> п,МО+1 ип,МО т л Чпш1 >

Л1

5л

_е

в) на бесконечности е

_ 1.

На рисунке приведен текст ядра программы, написанной в среде программирования Ма1Са^ реализующей эту разностную схему.

В первой строке программы объявляется массив динамических характеристик изучаемого процесса Тм 5: вектор Тм 0 хранит текущее приведенное время, Тм 1 - текущий приведенный радиус пузырька, ТМ2 - текущую приведенную скорость межфазной границы, Тм 3 - текущее давление в газе, Тм 4 - приведенную плотность теплового потока на межфазной границе, Тм 5 - производную давления по времени.

Во второй строке объявляется двумерный массив О м М для хранения приведенных температур фаз. В вектор Q передается информация о начальном состоянии системы, которая хранится в

предварительно заданном векторе 8. Далее программа выполняет два вложенных цикла: внутренний - по пространственной переменной ]=1,2,.. ,,М, и внешний - по временной переменной ¡=0,1,.. На каждом проходе временного цикла программа обращается к подпрограмме с передачей в

нее динамических характеристик из вектора Q. Подпрограмма - это численная реализация

одношагового метода Рунге-Кутта (в данном случае четвертого порядка точности); она рассчитывает динамические характеристики текущего состояния системы по ее предшествующему состоянию. Связь между динамической частью системы уравнений и теплофизической осуществляется через ячейку Тм 4. При каждом обращении к подпрограмме происходит заполнение одной строки матрицы Тм 5: эта информация используется при вычислении температур в узлах сетки. Данные о начальном распределении температуры хранятся в массиве 0, откуда она считывается в локальный массив е.

Во внутреннем цикле рассчитывается текущее значение пространственной переменной л, значе-

эе2

ние производной а _ —2, которая участвует в качестве аргумента функции FWG(an,pn,FP ,а , л), вычисляющей локальную скорость газа в пузырьке. Эта функция предварительно должна быть описана (правая часть уравнения (17)). Положение узла сетки по пространственной переменной отслеживается двумя операторами условного сравнения если л<1 (ш<МО), то расчет температуры ведется для газовой фазы по уравнению (19); при ш=МО (межфазная граница) программа выполняет условия на границе и подсчитывает плотность теплового потока, посылая результат для хранения в ячейку Тм 4. При ш>МО (жидкая фаза) температура в узлах рассчитывается по формуле (20). После перебора пространственных узлов сетки

эе„

программа удовлетворяет требованию

Эл

_ 0,

при

Л=0 в центре пузырька и начинает новый цикл по временной переменной.

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

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

1

а

п

1

х

+

х

+

2

1

а

п

S 1ST := TN, 5 — 0 eN, M — 0

е — 0

Q — S

for i е 0 .. N - 1

Ti 0 — DS (Q 0

Ti 1 — DS (Q 1

Ti 2 — DS (Q 2

Ti 3 — DS (Q 3

Ti 4 — DS (Q 4

Ti 5 — DS (Q 5

a — Ti, 1

w — Ti , 2

p — Ti , 3

q — Ti, 4

5p — Ti , 5

for j е 1 .. M - 1

h — j -5h

е

q —

if j < MG

1

i, j + ^"i, j

5h

ei+1, j — ei, j ...

+ 5t

FWG (a , w , p , FP , q ,h)-h'w ( ) ei, j+1 - ei, j

5h

P

g- 1___

g Pe 2'Ec 2

+ I-i '^j '5p

Ji, j

p' a

ei,j+1 - 2'ei,j + ei,j-1 2 -+ —

5h

i, j+, j 5h

if j = MG - 1

1

e

i , MG - ei , MG-1

5h

ei+1, MG+1 — ei+1 , MG + Ti , 4 '5h ' "

if j > MG

1

ei+1, j — ei, j + 5t

(-1) 'a

--h

2 1

3i , j + ^"i , j

5h

ei , j +1 - 2'ei, j + ei , j-1 2 - + —

+ 12

Pe 1' a

Ec 1 / w Re 1 la

5h

i, j+1 _ i, j 5h

2

2

1 2 Ha 'Ec 1 w

6 + 3 Re 1 4

h Re 1 h

( t ))i> Q — It j

augment

Рисунок: текст ядра программы, написанной в среде программирования MatCad

+

h

l

l

+

h

ei+1, о — ei+1 , 1

В качестве несущей фазы задавался жидкий галлий с теплофизическими свойствами при температуре 1000С, а дисперсной - воздух. Начальный радиус пузырька назначался равным 3 мм, начальное давление в газе p(0)=100 кПа. Исследовалась реакция системы на мгновенный сброс давления, поэтому параметр возмущения принимался равным k=-0,5, что соответствовало конечному давлению в системе 50 кПа. Показатель адиабаты воздуха задавался равным у=1,4, индукция магнитного поля составляла B=1 Т.

Критерии задачи при этих данных были: Re1=39000, Pe1=11360, Pe2=393, Ш=135, We1=616, We2=0,096, Ec1=1,2x10-4, Ec2=1,54x10-4. Параметры сетки задавались равными: M=20, N=1000, ^2, 1=50 и подбирались экспериментально из условия устойчивости конечно-разностной схемы. Контроль точности осуществлялся удвоением числа шагов по пространственной переменной.

На рисунке 1 показана зависимость давления р* в пузырьке от времени г *. Анализ кривой показывает, что время релаксации пузырька (время перехода системы к равновесному состоянию) можно разбить на две фазы: на первой фазе процесс характеризуется затухающими колебаниями всех динамических и тепловых параметров. На второй фазе пульсации динамических параметров исчезают, но тепловое равновесие еще не достигнуто, пузырек продолжает переход в равновесное состояние регулярным образом, так что в конце этой фазы между жидкостью и газом в пузырьке устанавливается термическое равновесие. Эта особенность привнесена магнитным полем, при отсутствии которого динамические и тепловые параметры системы переходят к новому состоянию равновесия синхронно. Данную особенность можно объяснить высокой интенсивностью гашения колебаний электромагнитными силами через механизм джоулевой диссипации. В результате чего резко сокращается осциллирующая фаза процесса, но температурные поля еще не успевают за это время прийти в состояние равновесия. С увеличением индукции магнитного поля длительность первой фазы, как показывают расчеты, сокращается. С быстрым исчезновением пульсаций пузырька термическое равновесие в системе не успевает возникнуть, и начинается вторая фаза -выравнивание температур между газом и жидкостью при уже достигнутом механическом равновесии пузырька. Необходимо отметить, что на этой фазе механизмы вязкой и джоулевой диссипации не работают по причине отсутствия пульсаций скорости.

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

На рисунке 2 показаны графки зависимостей температуры во времени в центре пузырька и на его поверхности.

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

1 0.8 0. 6 0.4 0.2

Р*

1/ V/ 1

0

200 400 600 800 1000

Рисунок 1. Зависимость приведенного давления р, от временной переменной 1. Переход к реальному времени осуществляется умножением 1 на г0 /20.

е.

0.9

0.85

0.8

0.75

0.7

2

1

г*

10

20

30

40

50

Рисунок 2. Зависимость приведенной температуры е2 от приведенного времени г, при расширении газового пузырька: кривая 1 - центр пузырька, кривая 2 -поверхность пузырька. Переход к реальному времени осуществляется умножением г* на го.

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

Графики, изображеные на рисунке 4, иллюстрируют пространственное распределение температуры в газе и жидкости в различные фазы первого периода колебаний пузырька.

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

Межфазный теплообмен полностью характеризуется числом Нуссельта

№ _-

2а () ЭТ2

Та-(Т2) Эг

где Т—температура на межфазной границе, а (Т2) -среднемассовая температура газа в пузырьке.

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

Ми _-

2

эе,

еЛ2_)1 - р*а3 Эл

На рисунке 5 показан график функции №=№(1;*): на первой фазе процесса число Нус-сельта испытывает пульсации, в целом оставаясь величиной положительной, а на второй - принимает почти постоянное значение, равное 12. Знакопо-ложительность числа Нуссельта говорит о направлении теплового потока из жидкой в газовую фазу: т. е. при расширения пузырька он непрерывно подогревается теплотой окружающей жидкости.

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

Под воздействием вязких сил к моменту времени 1 будет превращена в теплоту кинетическая энергия, равная

1 ¥ 2 4

дД()_ | _ 16^010и01 т((), V _ — ра3(1),

где

1Ц()_ |а* *.

Кинетическая энергия, превращенная в теплоту за счет джоулевой диссипации, будет равна

^2>

где

^(()_ |а3w2&*.

0

Теплота через межфазную границу определится следующим выражением

(»

0х()_ |^(ах_-4ра0Т012(0! ),

где

11(( На *Э02

л *

Л_1

Л_1

0 ЭЛ

Здесь через 1т ,1^1^ обозначены диссипатив-ные интегралы, на рисунке 6 показаны их графики.

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

0

0 Е

.995

0.99

.985

01 2

1

t*

1.1

10 20 30 40 50 Рисунок 3. Зависимость температуры в жидкости 91 от приведенного времени 1*: кривая 1 - на расстоянии одного шага (ЛЬ), кривая 2 - двух шагов от поверхности пузырька.

0.9

0.3

0.7

е

2 /

i

1 h

1

0.5

1.5

Рисунок 4. Распределение приведенной температуры по пространственной переменной Ь в системе «пузырек -жидкость» в различные фазы первого периода колебаний: кривая 1 - -^=0, 3 - —=Т, 2 - ^=Т/2. Кривые 1 и 3 относятся к двум последовательным максимальным расширениям пузырька.

30

Nu

20

10

t*

10 20 30 40 50 Рисунок 5. График функции Nu

30 28 26 24 22 20 18 16 14 12 10 8 6 4 2 0

0

I

2

t*

1

200

400

600

800

1000

Рисунок 6. Зависимость диссипативных интегралов ¡¡(Ч) от времени процесса: 1т (1) - кривая 1, I j (1) - кривая 2, Г (1) - кривая 3. Значения интегралов 1т (1) и I ^ (1) увеличены в 50 раз.

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

Проведенное исследование позволяет сделать вывод о том, что наличие магнитного поля

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

Список использованной литературы:

1. Нигматулин Р.И. Динамика многофазных сред. М.: Наука, 1987. - Т. 1, Т. 2. - 464 и 360 с.

2. Корн Г., Корн Т. Справочник по математике (для научных работников и инженеров). М.: Наука, 1974. - 831 с.