CC BY
45
УДК 532.529
Е.В. Резанова
Численное исследование динамики сферической газосодержащей оболочки
Е. V. Rezanova
Numerical Investigation of Dynamics of a Spherical Gas-containing Layer
Численно исследуется задача о динамике жидкой оболочки со свободными поверхностями и процесс диффузии газа в ней. Представлен алгоритм расчета задачи в диффузионном приближении и результаты численных экспериментов, иллюстрирующие основные зависимости решения от параметров модели.
Ключевые слова: сферический слой, диффузионное приближение, численное исследование.
БОЇ 10.14258/І2Уави(2013) 1.2-07
The problem of dynamics of a spherical layer with free boundaries and a process of gas diffusion in the layer is investigated numerically. In the paper an algorithm of computing of the problem in a diffusive approach and the results of numerical experiments are presented. These results illustrate the main effects of the solution on the parameters of the problem.
Key words: spherical layer, diffusive approach, numerical investigation.
Введение. Исследование течений жидкости в областях со свободными границами, а также процессов тепло- и массопереноса в них продолжает оставаться сложной и актуальной задачей [1]. Интерес к подобным проблемам особенно возрастает, если требуется изучить поведение жидкостей в условиях тепло- и массопереноса на границах раздела [2]. Построению математических моделей динамики вязкой несжимаемой жидкости с учетом диффузионного потока газа, растворенного в жидкости и считающегося пассивной добавкой, посвящены работы [2-5]. В данных работах предложена модель формирования сферических микробаллонов в условиях кратковременной невесомости. Разрешимость задачи в полной постановке в малом по времени доказана в [3]. В [4] построен численный алгоритм решения задачи в диффузионном приближении, а в работе [5] доказаны теоремы существования и единственности гладкого решения для тепловой задачи.
Следуя [3], будет изучаться задача о динамике сферической оболочки со свободными поверхностями, содержащей внутри газовый пузырь. Динамика сферического слоя определяется инерционными, тепловыми, диффузионными факторами,
Работа выполнена в рамках проекта № 7.3975.2011 Алтайского государственного университета (поддержан Министерством образования и науки РФ) и программы стратегического развития ФГБОУ ВПО «Алтайский государственный университет» на 2012—2016 годы «Развитие Алтайского государственного университета в целях модернизации экономики и социальной сферы Алтайского края и регионов Сибири», мероприятие «Конкурс грантов» (№2013.312.1.66).
а газ, растворенный в жидкости, является пассивной примесью. Условия на свободных поверхностях будут представлять собой кинематические и динамические условия, условия, определяющие баланс энергии на внутренней границе и режим теплообмена с внешней средой на внешней границе. Учитывается диффузионный поток массы через внутреннюю поверхность и закон Генри [3, 6] как соотношение, связывающее концентрацию газа на границе раздела с давлением газа вне области.
В данной работе рассматривается сферически симметричный процесс. Тогда только радиальная составляющая скорости отлична от нуля, и все физические величины зависят от расстояния от начала координат и изменяются со временем. Зависящие только от времени давление, плотность и абсолютная температура внутри газового пузырька связаны уравнением Менделеева-Клапейрона. Основное внимание уделяется построению численного алгоритма решения задачи в условиях, когда именно диффузионные процессы считаются преобладающими и определяющими динамику сферической оболочки.
1. Постановка задачи. Пусть И\ и Й2 -внутренний и внешний радиусы сферического слоя Й1^) < г < Й2(^) (рис. 1); у(Ь,г) — радиальная составляющая скорости жидкости; Т — температура; С — концентрация газа в жидкости.
В безразмерном виде систему уравнений Навье-Стокса, переноса тепла и диффузии газа в жидкости можно записать следующим образом [3]:
Краевые условия на свободных границах также можно записать в безразмерном виде. При г = Р\Ц) имеем:
= Біг
<Ш-і
сЫ
(5)
Р.Р,.шт+г^, ,6)
т = т,
Я’
(7)
Рис. 1. Геометрия области течения (сферическая оболочка). Здесь Ді, Дз — внутренний и внешний радиусы оболочки; Рг,п, Рд — давление внешнее и внутри пузырька; Тг,„,Т, Тд — температура внешней атмосферы, температура жидкости и температура газа в пузырьке соответственно; С — концентрация газа в жидком слое
1 /7/?^ И
-БЪа2Рд + БЪа3 - (■гпТ) = (8)
_ т?°-дТ , с; гйгк.2^2 с1 ,алл , сШ1л _ Ді_ + 5/м4{_[ДіТ -(-)] + СТ—},
С = НА(Тд)Рд . При г = Р-2 {і) имеем:
(1Р2
V = 57г-
сЫ
(9)
(10)
Р = Рт-23і^- + ^р^, (И)
Ко Ь ог
ді> до ^ 1 дР 5/г^— +г’тт- = -Ей-——Ь ОІ ог р ог
Ке от ог
^(г21>) = 0,
(1)
(2)
спдТ дТ 113 , (т,дТ
61г-т+1’д7 = р-е^гх{т)д7)+ (3) +2±аМТ)[ф2 + 2(1)\
1 дг 1 дг Ре а г2 <9г(' ( ^ дг ^ (^
Здесь р — плотность жидкости; Р — давление; г/, П,а — коэффициенты кинематической вязкости, диффузии и поверхностного натяжения соответственно; х - коэффициент температуропроводности. Отметим, что предполагается зависимость всех коэффициентов переноса от температуры. В результате приведения системы уравнений к безразмерному виду возникли следующие безразмерные комплексы: число Струхала Б1г =
>’* „^ Р*
число Эйлера Ей =
число Рейнольдса
Ре = -------, число Пекле Ре = ---------, число Пекле
V*
диффузионное Ре а = ——, а также параметры
Р*
сТ,р,
характерные значения физических величин, а с теплоемкость жидкости.
^ + Ми(Т - Тип) = 0,
С = НА(Тип)Р:п.
(12)
(13)
Условие, определяющее изменение массы газа в пузырьке, имеет вид:
с, Ард с; ^ сІРі
5',л- = -5',яГ',»~+ (14)
1 3
с)С
Рес1 /?! ОГ
Здесь Тд(1) — температура газа в пузыре; Рд(^) и Pvn.it) — давление газа в пузырьке и внешнее 4 о
соответственно; т = —пР^рд — масса газа в пузырьке; А — коэффициент в законе Генри; п —
показатель в законе Генри; Ни = ^ — число
к
Нуссельта; Н
А*Р?
р*
БІ
СТ 4= Р4= 7 4= ’^4
«2
кТ1
-. Звездочкой обозначены
Су ТП^'О % ф
аз = ——^—, «4 = ——, где су — теплоемкость Агкг%к к±*
газа (при постоянном объеме); к — коэффициент теплопроводности; /3 — коэффициент теплоотдачи.
Формулы (5) и (10) представляют собой кинематические условия на свободных границах, (6) и (11) — динамические. Закон Генри представлен выражениями (9) и (13), соотношение (7) есть условие непрерывности температуры при переходе через внутреннюю границу оболочки г = Р-1^). Теплообмен с внешней средой определен общим условием третьего рода (12) на границе г = Ро^).
Равенство (8) выражает баланс энергии на внутренней границе [3].
При t = 0 задаются начальные значения всех искомых функций:
ДЮ = ^l(O) < Г < Р2 (0) = i?20, РдО = Рд(0),
Тдо = Тд{0), v(0,r) = v0(г), Т(0,г) = Т0(г), (15) С(0,г) = С0(г).
Внутри пузырька газа 0 < г < i?i(t) задано уравнение состояния:
Рд = R РдТд, (16)
где R! = Др*Т*/Р*, R — универсальная газовая постоянная.
2. Постановка задач для нахождения концентрации газа, температуры, скорости изменения объема оболочки, положения свободных границ, плотности газа в пузырьке. С помощью новой безразмерной пространственной переменной х = (г3 — i?i(t)) • (Д20 — Д?о) 1 может быть осуществлен переход в фиксированную область [0,1] на каждом шаге по времени. В новых координатах уравнение диффузии (4) имеет вид [3]:
ЭС д . ,дС\
Ж = 15т’’]) а?1’ (17)
D(t, х) = 9Ped 1(i?f0 — Рло) 2' ■[Pio-Rlox + RKt^DiTit)).
Поскольку имеют место следующие начальное и граничные условия для концентрации (см. (9), (13), (15)), для функции C(t,x) получим:
С(0, х) = С'о(х), C\x=1 = Cvn, С\х=0 = Си
где Cvn = A(T\r=R2)P^n.
Уравнение теплопереноса (3) без учета диссипативной функции в новых пространственных переменных примет вид [7,8]:
дТ д , дТл X(t,x) = 9Pe~1(R% - Д?оГ2-
(18)
■1р10-р10х + рШ/3хШ).
Начальное и граничные условия (15), (7), (12) в данном случае примут вид:
Т(0,ж) = Т0(х), Т\х=0=Тд,
\ , дТ (гр\ гр ч ^20 ^10
dx I-=i+Wm(TI-=i Т™> 3(дзо _ Д30 + д3)2/3
Тд — неизвестная функция, зависящая от времени, которая также будет найдена в ходе решения
задачи. Помимо этого, при х = 0 имеет место соотношение (8) в виде
1
dRl
Sha2Pg^ + Sha3^(mT) = (19)
ЗД1 дТ\ , ои Г d Го2^2 d /^м , (Ш1л
= й|0-й?0&и-”+“а4{й|Й1Т df(T)1 “л“}'
Заметим, что v = r~2V(t) (см. (2)). Тогда рассмотрим следующую задачу для системы обыкновенных дифференциальных уравнений для определения функций V(t), Ri(t) и pg(t):
f = ^(й + я;,№ + я1)^+ (20)
_ ±V(T)V
R2 — ^1 R&
•(i?2 + R2R1 + ^1) d2 d2 j ^ У(0) = Vo;
^2^1
T r 1 (IR2 _ 1
Ц' ~dT= Щ
i?i(0)=i?io, P2(0) = R2q,
dt V RV dt V RV t>0; (21)
dt ~ Ri dt Pg + Ped R20 - Д3П ^(T)' (22)
L10
0.
9(7
"9ж"^Ж = 0’ ^ > 0, Рд{ 0) = /ЭдО-
Здесь введены следующие обозначения: А(Т) =
НА{Т), Р'д = Рд Б, Ж = Si • 5.
В данной постановке осуществлен переход к новому безразмерному времени делением на число Струхала.
Если выбор характерных значений физических величин осуществлен согласно [4] для системы «жидкость — газ» типа «жидкое стекло — углекислый газ», то характерное время процесса выбирается равным 1.4 с, характерное значение радиуса жидкой оболочки г* = 0.05 см, характерная скорость г;* = 1 см/с, характерная температура разогрева внешней среды Т* = 1673 К, характерное внешнее давление Р* = 1013250дин/см2 =
1 атм, характерная плотность жидкого стекла в диапазоне температур Т = 1273К 4- 1673 К, р* = 2 г/см3, значения характерных коэффициентов таковы: V* = 36 см2/с, поверхностного натяжения <т* = а(Т*) = 280 дин/см. Показатель в законе Генри выбран равным п = 0.5. Безразмерные комплексы принимают следующие значения: Як = 0.036, Ре = 0.0014, Рел = 1219.5, Ре = 5, N4 = 0.681, вг = 0.0055, Р' = 6235, 5 = 703.65. В случае, если -у* = г*/£*, то вк = 1, Ре = 0.18, Рел = 44, 5 = 19545.7, Де = 5 • Ю^5.
3. Общая схема численного решения.
Общая схема решения задачи (17)—(22) состоит в осуществлении следующих этапов:
1. Переход на новый временной слой к + 1 начинается с расчета Rk+1 ,Vk+1, рк+1 методом Рунге-Кутта четвертого порядка точности для системы обыкновенных дифференциальных уравнений (20)—(22). Внешний радиус оболочки i?2+1 вычисляется исходя из закона сохранения объема оболочки R2(t) — R\(t) = Що ~ R-W ( см. кинематические условия на свободных границах (5), (10)).
2. С помощью неявной разностной схемы для уравнения (17) вычисляется концентрация газа на временном слое к + 1.
3. Температура жидкости определяется с помощью неявной разностной схемы второго порядка аппроксимации для уравнения (18).
Также был протестирован алгоритм, в котором метод Рунге-Кутта применялся для вычисления каждой из функций i?i(t), V(t), p(t) по отдельности. Численные эксперименты показали незначительную разницу в полученных значениях искомых величин.
Численная схема метода Рунге-Кутта расчета функций i?i(t), V(t), p{t) может быть представлена в виде:
Vk+1 = Vk + T-{kl + 2к2 + 2к3 + *4), ь Rk+1 = Rk + — (</i + 2q2 + 2q3 + (/4), pk+1 = рк + tt(«;i + 2 w2 + 2 w3 + w4),
D
где hi, qi, Wi (i = 1,..., 4) вычисляются следующим образом:
k1 = K(tk,Vk,Rk,pk),
т2 = + г/2, Ук + к\ — ,Р\ + (/1 —, рк +101—),
газ = \¥(^к + т/2,Ук + к2 — , + ц2~, рк + ги2—),
= \У{Ьк + г, Ук + к3т, Нк + д3т, рк + ии3т).
Для функций К,С2,\У имеют место зависимости вида:
К(^,Ук,К'1,рк) =
= +(«й2)^р^1+
1 1
( Т?к і Т?к^
.с7
+ R-ep*[P°-PL-2Si'aiT
RWi _ ±v{Tk)Vk (Rl)2 + R%Rhi + (Rl)2
Rk - Rk Re
(R%)2№)2
Q(tk,Vk,Rk1) = Vk
(ДЇГ
\¥{ік,Ук, Рк,рк) = -—ркук+
П Ґ^к г^к тук
—і——— В(Тк) _________-і-____^_______
Рег, К > ~к _ „к Е>3 оЗ
ґе<! х2 Х1 -^20 Л10
На каждом временном слое к вводится итерационный процесс нахождения концентрации газа С и температуры Т в оболочке. Для численного решения уравнения диффузии (17) строится неявная разностная схема второго порядка аппроксимации по пространственной переменной [4]:
S-1S+1 s-is 1 _ fS+l _____________ y^»S+l
1 г7=Г i+1 і
------------------- — ------------------------1.--------------
т~с hi hi^i
___ (7S+1 _ Cs+1
Di *'
(23)
І— 1
k2 = K(tk + + k\ — ,R,{ + qi —, pk +wi-),
&3 = K(tk + + k2 —, Rk + <72 2 > pk + w2 2 );
kA = К(tk + t, Vk + к3т, Rk + q3T, pk + w3t), qi = Q(tk,Vk,Rlpk),
q2 — Q(tk + —, Vk + ki~, Rk + <71-, pk + wi-), 43 = Q(tk + —, Vk + k2—, Rk + q2~, pk + w2^
q4 = Q(tk + r, Vk + к3т, Rk + q3T, pk + w3t), w1=W(tk,Vk,Rk,pk),
где С? — значения концентрации на в-той итерации, = хн - ж*_1, Ь = 0,5+ К+1), Di = 0, 5[Б(гк+1,хн^1) + Б(гк+1,хн)\, ts+1 =ts +тс.
Для реализации схемы (23) используется метод прогонки [1]. Применяется критерий сходимости итерационного процесса вида:
тах\С^+1 — С-\/тах\С-+1\ < е,
где е > 0 — малый параметр (см. [4]). Предполагается параболический профиль начального распределения концентрации с максимумом в центре слоя.
Для уравнения переноса тепла (18) строится аналогичная разностная схема
грТП-\-1 гргп 1 грТП-\-1 _ грТП-\-1
* ~Г< =и%+1~1±Чг^---------(24)
тт Ы т Ы+1
Здесь Хг = 0, 5[х(#й+1,ж*_1) + х(1к+1 > ж*)]; ^т+1 = 1т +тт- Для вычисления \ используется выражение (18).
Для реализации схемы (24) применяется метод прогонки с параметром, где в качестве параметра выбирается неизвестное значение температуры Т при х = О {Т\), совпадающее с Тд на каждом временном слое. Данной схеме (24) соответствует система линейных алгебраических уравнений
I ^7 ^ 7 /'} 1 7 . * 2, ...I.
Поиск Тг В виде Тг = 7-\-СцТг^\-\-рг ПОЗВОЛЯвТ выразить все Т.,, через Т:
Р = 7+ Дг, а уравнение (19) позволяет найти Т.
4. Результаты численных исследований задачи в диффузионном приближении.
Пусть температура жидкости и температура газа в оболочке равна температуре газа вне оболочки Туп{1), которая зависит только от времени.
Численные эксперименты проводились таким образом, чтобы выявить влияние внешних и внутренних факторов на динамику оболочки и диффузионные процессы в ней. Пусть Дю = 0.02 см, Д-2о = 0.05 см, и имеют место зависимости коэффициентов от температуры следующего вида [4]:
V = а,иехр(Ъи/Т), В = а,£>Т + Ьр,
а = ааТ + Ьо-, А = А/Т.
Для жидкого стекла имеют место следующие значения безразмерных коэффициентов: а„ = 0.18 • 10~5, Ь„ = 13.23, а0 = 1.195, Ъв = -0.195, аа = 0.299, Ъа = 1.299.
Сферическая оболочка, содержащая газ, помещается в разреженную атмосферу того же газа, поддерживаемую при определенной температуре.
Рис. 2. Графики зависимости внутреннего радиуса оболочки от времени при различном количестве газа внутри пузырька. Внешнее давление Р=0.03 атм.
Пусть внешняя атмосфера нагрета до определенной температуры = 1171 К, которая во время протекания процесса не меняется. Внешнее давление будет также постоянным, равным Р = 0.03 атм (рис. 2) или Р = 0.01 атм (рис. 3).
При больших значениях начальной плотности газа рдо = 0.92 • 10~3г/см3 (см. рис. 2 и 3, сплошные линии) динамика сферической оболочки и процессы в ней протекают более интенсивно. Так, например, при Р = 0.03 атм максимальное значение внутреннего радиуса оболочки Д^,гаж = 0.045 см. В случае, если начальная плотность газа в пузырьке имеет значение рдо = 0.46 • 10~3 г/см3, то наблюдается менее интенсивные процессы; /■'„ = 0.036 см.
Рис. 3. Графики зависимости внутреннего радиуса оболочки от времени при различном давлении вне пузырька. Внешнее давление Р=0.01 атм.
При одних и тех же значениях температуры внешней атмосферы (Т = 1171 К) и количестве газа в пузырьке (начальное значение плотности газа в пузырьке равно рд0 = 0.92 • 10~3г/см3 (рис.
2 и 3, сплошная линия) или рд0 = 0.46 • 10~3г/см3 (рис. 2, рис. 3, штриховая линия) и различных значениях внешнего давления (Р = 0.03 атм, Р = 0.01 атм) имеем обратный эффект: с увеличением давления интенсивность процесса снижается, конечные и максимальные значения внутреннего радиуса оболочки уменьшаются. Так, при рд0 = 0.46 • 10~3г/см3, Р = 0.03 атм максимальное значение внутреннего радиуса оболочки равно 0.036 см, значение на момент времени I = 0.5 с -
0.032 см. В случае, если давление имеет значение Р = 0.01 атм, то это максимальное значение больше 0.038 см. Значение внутреннего радиуса оболочки на момент времени I = 0.5 с равно 0.035 см.
Заключение. Описан алгоритм численного исследования полной задачи о динамике жидкой сферической оболочки со свободными поверхностями и процессов переноса тепла и газа в ней в условиях невесомости.
Представлены результаты численных исследований задачи в диффузионном приближении.
Проведены численные эксперименты для исследования зависимости динамики оболочки и процесса диффузии в ней от внешнего давления и количества газа в пузырьке.
Автор выражает искреннюю благодарность научному руководителю О.Н. Гончаровой за постановку задачи и обсуждение результатов.
Библиографический список
1. Андреев В.К., Гапоненко ЮА., Гончарова О.Н., Пухначёв В.В. Современные математические модели конвекции. — М., 2008.
2. Гончарова О.Н. Моделирование течений в условиях тепло- и массопереноса на границе // Известия АлтГУ. - 2012. - №1/2(73).
3. Гончарова О.Н. Математическая модель формирования сферических оболочек в условиях кратковременной невесомости //Г идродинамика быстропротекающих процессов. - 1987. - Вып. 82.
4. Гончарова О.Н., Пухначёв В.В. Диффузионное приближение в задаче формирования сферических микробаллонов в условиях кратковременной невесомости // Моделирование в механике. - 1990. - Т. 4 (21), №5.
5. Гончарова О.Н. Глобальная разрешимость задачи о формировании сферических микробаллонов / / Вычислительные методы прикладной гидродинамики. — 1993. — Вып. 106.
6. Пригожин И. Введение в термодинамику необратимых процессов. — Ижевск, 2001.
7. Гончарова О.Н., Закурдаева А.В., Резанова Е.В. Моделирование динамики и процессов тепло- и массопереноса в сферическом слое жидкости со свободными границами. — Пермь, 2013.
8. Закурдаева А.В., Резанова Е.В. Математическое моделирование процессов диффузии и теп-лопереноса в газосодержащей оболочке // МАК-2013. — Барнаул, 2013.
