CC BY
17
А.П. Васильев
ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ СИСТЕМЫ УРАВНЕНИЙ ДИНАМИКИ И ТЕПЛООБМЕНА ПРИ ПУЛЬСАЦИЯХ ГАЗОВОГО ПУЗЫРЬКА
Рассматривается система дифференциальных уравнений, описывающая теплофизические и динамические процессы при осцилляциях газового пузырька. Предложен алгоритм численного решения и приведены результаты решения.
В работе [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)
где к - параметр возмущения.
Уравнение энергии жидкой фазы
ЭЭ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)
1 Эе,
- 3ур*
-=1
&* Ре2Ес2 а; Э-
Уравнение для скорости в газовой фазе
р* у—1 1 Эе 2 1 ф
(2)
(3)
—
,(-д *)=
Ре2Ес2 у р*а* Э- 3у &*
-' (4)
Уравнение энергии для газовой фазы
Здесь индекс «1» относится к жидкой фазе, а «2»- к газовой; р0 - истинная плотность, с1 -удельная массовая теплоемкость, 11 - теплопроводность, —1 - скорость, р1 - давление, —1 - скорость, Т1 - температура 1-ой фазы, V, т - кинематический и динамический коэффициенты вязкости
Эе Эг
2 + -2*— -а* Эе2
Э-
у-1 р*
Э2е
у Ре2Ес2
Э-
2+ 2 Эе^ 2 - Э-
у—1 е2 ф* у р* &*
(5)
Система уравнений (1)-(5) является замкнутой. Неизвестными в ней являются: а * (г *) - при-
а
1
+
2
-
г
*
а
веденный радиус пузырька, р. (1 *) - приведенное давление, -№2*(т|,1 *) - приведенная скорость газа, 01 (т,1 *), 02 (т,1 *)- приведенные температуры в фазах, причем все неизвестные функции определены так:
а * =
а (() = р(()
р* = ■
^ * = -
w
)
Ро
0 = Т1 (г,1) 0 = т2(г-')
1 То , 2 То ,
где и0 = а0/10 - характерная скорость, а все другие величины с индексом «0» относятся к начальным параметрам, - - показатель адиабаты.
Краевые условия для системы уравнений (1)-(5) задаются в виде
1,= 0: а* = 1, а * = w, = 0, 81 (т|,0)=1, 82 (|,0) = 1, (6)
Т = 0: ^ = 0, 82 <<
Э8 2
Эт
Э81
(7)
Т = 1: 12 ^ = 1, 82(1,1 .) = 8,(1,1 ,), (8)
т = ¥ : 81 (~,1 *)= 1. (9)
Критерии и числа задачи определены так:
На = а0В _ - число Гартмана,
Vт
Яе1 =
а0 и0
- число Рейнольдса,
р0 и2
^е1 = ^ / - число Вебера,
Ч а0
где Ч - коэффициент поверхностного натяжения,
Ре1 =
а0^ х, Ре2 =:
а0 и0
- числа Пекле,
Д2а,
Д1
Др. Д1,
= Ба (а,,а,
= (а ,,а , ,р , ,qД
(10) (11)
где через а* обозначена производная q * =
Э|
Заменим производные по временной и пространственной переменным конечными разностями вида:
# (т,1) = f (т,1 + 51)- f (т,1)
Э1
51
Эf (т,1) = f (т + 5т,1)-f (т,1)
(12)
Эт
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 ( + 2к2 + 2кз + к4), (15)
6
Рп+1 = Рп +1 (т1 + 2т2 + 2тз + т4), (16) 6
где коэффициенты Ц и ш1 (1=1,2,...,4) находятся через функции Ра (а* ,а* ,р*) и Рр (а* ,а* ,р* ,q*) согласно выбранной схеме [3], например,
к1 = Ра ОЧп^^Рп )51, т1 = Рр.^^Рп^п Ж" . Уравнение (4) на сетке ш5т, п51; записывается так
w & —- ^ Бр (an,Wn,Рn,qn )5т, (17)
Ре2Бе2 - Рпап 3-
а уравнения энергии для газовой и жидкой фаз в виде
8п+1,т = 8п,т + 51
wп2)- m5тwn 8п:т+1 -8 ап 5т
3(2)
п,т
+--рр (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т
(18)
т=1
V
+
+
(1) _ 0(1) п,ш
9п+1,ш _ 9
+ 51
(ш5л)2
- ш5л
(1)
9п,т+1 - 9
(1) п,ш
5л
?е1аП
9П1,)т+1 - 29
(1) п,ш
+ е:
(1)
(5л)2
еП1,)ш+1- е
(1) п,ш
+12
ш5л
Ес1 Яе1
5л
\2
(ш5л)6
2 Иа2Ес,
+1
(19)
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);
б) на межфазной поверхности
е(2) -е(2)
qn _■
е (2> -е1
п,МО п,МО-1
е(2) _е(1) е(1) _е(1) +^а 5л •
п,МО ип,Мв> п,МО+1 ип,МО т л Чпш1 >
Л1
5л
_е()
в) на бесконечности е
(1)
_ 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 (t , e)
Рисунок: текст ядра программы, написанной в среде программирования 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>
(()_.[ _ 8 р°в2та0 и2(0^ (()
где
^(()_ |а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(t*)
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 с.
