CC BY
33
Гималтдинов И.К.
Стерлитамакск ая государственная педагогическая академия
О ЗАТУХАНИИ ИМПУЛЬСНЫХ ВОЗМУЩЕНИЙ В ТРУБЕ С ПУЗЫРЬКОВОЙ ЖИДКОСТЬЮ СО СТУПЕНЧАТЫМ РАСПРЕДЕЛЕНИЕМ
ПУЗЫРЬКОВ ПО СЕЧЕНИЮ
Изучена эволюция волн давления в трубе, заполненной газожидкостной средой, при ступенчатом распределении пузырьков по сечению трубы. Результаты расчетов сопоставлены с экспериментальными данными. Показано, что из-за появления поперечных течений при неравномерном распределении пузырьков происходит более интенсивное затухание импульсного давления по сравнению со случаем гомогенного распределения.
Постановка задачи.
Рассмотрим двумерные осесимметричные волновые возмущения в трубе с газожидкостной смесью при неравномерном распределении пузырьков по сечению трубы. Такое волновое движение может реализоваться, например, при торцевом воздействии поршнем давления или скорости на систему с осесимметричным распределением пузырьков в канале. В частности, такая ситуация имеет место, когда пузырьки находятся вблизи стенки трубы в кольцевом объеме или расположены в приосевой зоне в виде газожидкостного ядра, или расположены в виде сферического кластера. На рис. 1 представлено схематичное изображение для первой ситуации, когда осевая зона с радиусом Rp заполнена «чистой» жидкостью, а периферийная часть (цилиндрический слой шириной ДО = Rc - Rp) - смесью с однородной объемной концентрацией газа а^ (аgo << 1) и одинакового радиуса a0 пузырьками.
1. Основные уравнения
Полагая общепринятые допущения для пузырьковых жидкостей, запишем систему макроскопических уравнений масс, числа пузырьков, импульсов и давления в пузырьках в приближении цилиндрической симметрии [1]
d рі иг
-р- + Рі — + Рі dt г
Эиг + Эи z Эг Эz
Эи г Эи 7 —- + —7 Эг Эz
= 0(1 = и), = 0,
dn и г
— + п—1- + п dt г
V /
d иг Эpl А d и7 Эр1 А
Р—L + ^ = 0, р—- + = 0,
dt Эг dt Э7
=
ар1
dt
da
Зур
g 3( у-1)
w - ——- q,
'а э э э
— = — + иг — + и7 — dt Эt Эг Э7
(1.1)
dt ’
«і +аg =1 аg = 43Рпа3, р1 =р0а1, P = Pg +Pl,
где а - радиус пузырьков, у - показатель адиабаты для газа, р{ - давления фаз, р0 - истинные плотности фаз, а{ - объемные содержания фаз, q - интенсивность теплообмена, п - число пузырьков в единице объема, ’ - радиальная скорость пузырьков, и г и - радиальная и осевая составляющая скорости. Нижними индексами {= 1,р отмечены параметры жидкой и газовой фаз. При описании радиального движения в соответствии с уточнением, предложенным в [2], будем полагать, что w = + ’А, где wR опре-
деляется из уравнения Релея - Ламба, а ’А определяется из решения задачи о сферической разгрузке на сфере радиуса а в несущей жидкости в акустическом приближении
3 2 „ ’о Рр - Р1 ^- р1
а 0 + 4у, —0 = - ё 1 ё
2
dt
^ =-
Рисунок 1. Схема задачи 1 - поршень, 2 - рабочая область
а
а
где V1 - вязкость жидкости, С1 - скорость звука в «чистой» жидкости.
Будем полагать, что жидкость является линейно сжимаемой, а газ калорически совершенным
р1 = Ро + С1 (р1 - рю), РБ = Р80ТБ , где Я - газовая постоянная. Здесь и в дальнейшем индексами 0 внизу снабжены параметры, относящиеся к начальному невозмущенному состоянию.
Тепловой поток д задается приближенным конечным соотношением
q = Nul
Tg - T0 2a
T
g
Pg
Po
Nu= VPe, Pe> 100, Nu = 10, Pe< 100,
Pe = 12(g-1)
aw
kg =■
, _ T0| “g ' CgPg ’
где T0 = const - температура жидкости, cg и 1 g - теплоемкость и теплопроводность газа, Nu и Pe- числа Нуссельта и Пекле.
Из этой математической модели в частном случае при a go = 0 следует волновое уравнение для линейно сжимаемой жидкости. При исследовании взаимодействия волн в «чистой» жидкости с пузырьковой средой это обстоятельство в свою очередь позволяет использовать сквозные методы расчета.
2. Метод численного расчета. Для численного анализа задачи об эволюции волн давления в трубе, заполненной газожидкостной средой, при неравномерном (ступенчатом) распределении пузырьков по сечению трубы и по взаимодействию волны давления со сферическим пузырьковым кластером удобнее пользоваться системой уравнений, приведенной в разд. 1, записанной в лагранжевых координатах. Это, в частности, связано с тем, что в лагранжевых координатах первоначальные границы неоднородностей остаются неподвижными. Из уравнений, приведенных в разд. 1, после некоторых преобразований можно получить следующую систему в лагранжевых переменных:
dt
диг
“дТ
_1_L Г0 Jp
\
- _1_L Г0 Jp
dpj дг dpj дг
дг0 dz0 dz0 дг0
dpj дг dpj дг
дг0 dz0
dz0 дг0
dz
-----= uz,
dt z
дг
dt = у г;
др1_
dt
C2P0 CJ pj
1 -a
g0
3a gw
dJ
dt
pi0 , ag
^2 J
p0J
w
//
dt
3gp
g 3(g-1)
-w - ——- q,
(2.1)
да
at"
5wr
dt
wA =
pj
Pj
pj 3 2 w
----2wR -4v
a
P0CJ agg3
dag. = 3a!w ag dJ
dt
g
J dt
J=
dz дг dz дг
dz0 дг0 дг0 dz0
dJ
dt
duz дг + dz диг duz дг
dz0 дг0
dz0 дг0
дг0 dz0
dz диг дг0 dz0
у J
где 70 и г0 - лагранжевы переменные, в качестве которых берутся начальные эйлеровы координаты, I - якобиан перехода от лагранжевых к эйлеровым переменным [2]. Система (2.1) решалась численно по явной схеме. Приведенные уравнения, из-за учета меж-фазного теплообмена и акустической разгрузки пузырьков, являются системой с достаточно сильной естественной диссипацией, поэтому не требуется вводить искусственную вязкость.
Приведем принцип построения разностной схемы, принятой для решения данной задачи. Для аппроксимации дифференциальных уравнений используем равномерную шахматную сетку
г0Г Tk )
50i+1 = z0i + hz0,
j = г0і + h^
z0i+X VУг
znn = 0
г00 = 0,
= z0i + 0,5hz0,
= г0і + 0,5hгo,
i = 0,1, j = 0,1,.
,Mj -1
,Мг -1
‘0M.
= Rc
tk = kx, к = 0,1,2,...
где И70, Иг0 и х - соответственно шаги по координатам 7 , г и времени. К узлам сетки (7о1,Г)Ик ) будем относить сеточные функ-
ции скоростей
ij
и эйлеровых пере-
к к
менных zij и г- , к «полуцелым» точкам
сеточные функции всех
а
а
а
= w = wr + wA
P
g
R
P
g
а
а
0
г
г
0
г
г
0
z
г
t
z
г
k
Oi+
остальных параметров. Такая аппроксимация обеспечивает устойчивость решения волновых задач в однофазных системах (жидкостях и газах) конечно-разностным методом [2].
Для проведения численных экспериментов примем следующие начальные и граничные условия. Условия при (t=0), z > 0, соответствующие исходному состоянию покоя неоднородной пузырьковой смеси в трубе, запишутся в виде
Р = Po , Ur = Uz = 0 , P = Po .
Кроме того, для газожидкостного кольца необходимо добавить:
К < ro < Rc:ag =ag0> P = Pl0 (-ag0 )>Pg = Po,a = ao,W = °
[o < Г0 < Rp: ag0 = 0, p=p°,
а для газожидкостного ядра:
Rp < r0 < Rc : ag0 = 0 P = Pi0,
0 < 10 < Rp :ag = ag0,P = P°)(-ag0),Pg = P0,a = a0,W = 0,
В случае однородной газожидкостной смеси в трубе радиуса Rc имеем:
{о<r<Rc: ag =ago, P = P0(-ag0),Pg = Po, a = ao, w=0, На оси симметрии ( r0 = 0 ) и стенке трубы ( r0 = Rc ) принимаются условия непротекания жидкости ur = 0 ,
На торцевой границе приняты два вида граничных условий, а именно условие на жестком поршне
z0 = 0: 0 < r0 < RC Uz (r0,t )=u0z (t) (2.2)
или условие на поршне давления
z0 = 0: 0 < r0 < RC Pl (r0,t )= Pl0 (t) (2.3)
Отметим, что задание скорости поршня, в частности, соответствует воздействию на среду метанием жесткого ударника, а задание давления на границе z0 = 0 - разрыву мембраны между камерой высокого давления, заполненной газом, и рабочей камерой, заполненной исследуемой системой.
На другой торцевой границе (z0 = Lz) условие ставится в зависимости от конкретной задачи. Если на z0 = Lz среда граничит с газом, то для волн, распространяющихся по пузырьковой смеси, эта граница эквивалентна свободной поверхности и здесь в качестве граничного условия задается условие постоянства давления P(Lz,r )= p0 . Если же на границе z0 = Lz движение жидкости ограничено жесткой стенкой, то для скорости ставится условие uz = 0 . Возможен третий тип условия. Для того, чтобы воз-
мущения «уходили» из области расчетов, не отражаясь от границы 70 = , на этой границе
необходимо поставить неотражающие условия, в частности этого можно добиться, используя импедансное соотношение, связывающее амплитуды давления и скорости [3] Др1 = р0С1,
где Др1 и Д^ - текущие значения возмущений давления и осевой скорости для лагранжевой системы для границ расчетной области.
3. Результаты расчетов.
Проведены расчеты применительно к экспериментальным данным [4], полученным в ударной трубе длиной 1,5 м с внутренним диаметром 0,053 м. Внутри рабочего участка трубы располагалась тонкостенная лавсановая трубка диаметром 0,0375 м. Диаметр лавсановой трубки был подобран таким образом, чтобы площадь поперечного сечения внутри нее была равна площади кольца между лавсановой трубкой и стеной рабочего участка. Рабочий участок заполнялся жидкостью и насыщался пузырьками газа через генератор, расположенный в нижней части трубы. Опыты проводились для трех различных структур пузырьковой среды. Пузырьки подавались равномерно либо по сечению всего рабочего участка (гомогенная среда), либо в кольцевой зазор между лавсановой трубкой и твердой стеной рабочего участка (пузырьковое кольцо), либо внутрь лавсановой трубки (пузырьковое ядро). После создания в трубе пузырьковой смеси с заданной кон-рисурацией, по системе производился торцевой удар с помощью жесткой пластины, площадь которой равна сечению трубы. Метание пластины организовывалось с помощью устройства, состоящего из электромагнитной катушки, которая отталкивает пластинку при прохождении в ней импульса электрического тока. Среднее по сечению трубы объемное содержание газа во всех трех случаях составляет аБ0 = 0,005, при этом истинные объемные содержания пузырьков и их радиусы для обоих случаев расслоенной структуры были равны а б0 = 0,01, а0 = 5,3 10-3 м . В качестве рабочей жидкости использовался 50 %-й (по массе) раствор глицерина в воде, в качестве газовой фазы - фреон-12 или азот. Волны давления регистрировались пьезоэлектрическими датчиками D1, D2, D3 и D4, расположенными вдоль рабочего участка на расстояниях z0 = 0 , 0,25, 0,76 и 1,25 м от места входа волн давления соответственно. Численные расчеты применительно к экспери-
ментам проводились при следующих теплофизических параметрах:
^ = 595 Дж/кг • К , рg0 = 5,06 кг/м3 ,
Xg = 0,0097 м • кг/к • р3 , г = 1,14 (^ = 3,2 10-6 м2 /с) (для фреона-12);
^ = 1041 Дж/кг • К , р ^ = 1,15 кг/м3
1 g = 0,0256 м • кг/к • р3 , у = 1,4 (^ = 2,2 10-5 м2/с) (для азота); р 1 = 1115 кг/м3 , V1 = 6 10-6 м2/р, то = 293К , р0 = 0,1 МПа . При расчетах для И70, Иг0 и т применялись сле-
дующие значения: И70 = Иг0 = 10 3м, т = 10 7 р
3.1 Гомогенное распределение пузырьков. На рис. 2, б) приведено сравнение расчет-
0.49т
0.51т
0.25т
1.5 т
ных осциллограмм (правые картины) с экспериментальными из [4] (левые картины) для гомогенной смеси (схема рабочего участка представлена на рис. 2, а)) с пузырьками из фреона-12. На рис. 2, в) и г) представлены экспериментальные и расчетные данные о затухании импульса давления (амплитуды первого всплеска) в пузырьковой жидкости с фреоном (в) и азотом (г) при двух значениях амплитуды начального сигнала. Импульс, полученный в экспериментах метанием жесткого ударника, в расчетах моделировался заданием скорости поршня на границе 70 = 0 в виде «треугольника»
б)
&Р=Р,-Ро 0,272-
МПа
ОН
т
1000 2000 0 1000
0,117
4000 5000 4000 5000
Ар
ЛРо
1
0,8
0,6
0,4
0,2
0
в)
О 1 ♦ 2 А 3 А 4
0,2
0,4
0,6
0,8
1,2
1,4
г. м
Ар
Л/?о
1
0,8
0,6
0,4
0,2
г)
О 1 ♦ 2
Д 3 □ 4
0,2
0,4
0,6
0,8
1,2
1,4
г. м
Рисунок 2. Затухание волн давления по длине рабочего участка для гомогенной структуры среды. а) схема задачи; б) осциллограммы для датчиков В1-Э4; в) фреоновые пузырьки 1 - эксперимент, 2 - расчет для Др0 /р0 = 2,5 3 - эксперимент, 4 - расчет для Др0 /р0 = 9,5 г) азотные пузырьки 1 - эксперимент, 2 - расчет для Др0 /р0 = 1,9 3 - эксперимент, 4 - расчет для Др0 /р = 14 Расчеты проводились при р0 = 0,1 МПа, Т0 = 300 К .
г, м , 0.02
0.01
-0.01
-0.02
0.02
0.01
t = 0,4 мс .жттп,,, дМ1 ,ппи, „п, л
* |||| II >11 Щ| ||/ II I I I
«V 'и }■т1ипУ V/
Ш \ А Д^Г Д1\ Л ? Т\ Л
щщ * ч Ш Ш и щ ^ ^
_^#Шгтт .^пл м/.. л /л г /\
4П
;Г***шииш
^^ИНИШ'' ‘ "Ш'7 V ‘¿Ч||^ VII' ■ ‘IV V V чу ^^ииншит иЛ \п О Ш |п ЛЧ
^21’Г, 7 ”,ТИТ7,71ТГ,7ТГГ; ~ /ПТ, ~ПТ7, ТТ~ТГ,Т1Г,
К
Ц1ИИШ1
^Ц-Цццц
/III П Ш1 III и и и к
м
0.25
0.5
-0.01
-0.02
^ / = 1.2 мс
- - -
\тп>’Ч -? ' ^й^^л.А-^,тппм^м«гмт-'VI
^п^^*^-^,от',т,м'йгичгитиу •?■
41 у V1*
^иииип^ и А/!
Рисунок 3. а) показания датчиков Б1-Б4; б) затухание волн давления по длине рабочего участка (фреоновые пузырьки):
1 - эксперимент, 2 - расчет - гомогенная структура среды (Др0 /р0 = 2,5)
3 - эксперимент, 4 - расчет - газожидкостное кольцо вблизи стенки (Др0 /р0 = 2,1);
в) поле скоростей в различные моменты времени;
г) эпюры давления в различные моменты времени.
0.25
0.5
20 = 0:
(1 )=
Ди —, 0 < г < г.
-Ди •-
11 12
г1 < г < г2 + г1, (3.1)
0, г > г2 + г1,
где Ди2 - амплитуда возмущения скорости, ^ и г3 - периоды подъема и спада скорости ударника соответственно. При расчетах для них приняты следующие значения Ди2 = 0,7м/р, г! = 1,8 -10-4 с , г2 = 2,2 • 10-4 с . Как следует из рис. 2, (б), (в) и (г) наблюдается неплохое согласование расчетных данных с экспериментальными. В случае пузырьков азота, у которых коэффициент теплопроводности больше, чем у фреона, затухание происходит более интенсивно. Видно также, что с ростом амплитуды первоначального сигнала интенсивность затухания усиливается (точки 3 и 4 на рис. 2, в) и г)).
3.2. Кольцевая структура На рис. 3 представлены экспериментальные и расчетные данные по динамике волн давления в случае кольцевого расположения фреоновых пузырьков в трубе (схематическое изображение приведено на рис. 1, а)). Рисунки (а) соответствуют осциллограммам для датчиков D1-D4, левые картины иллюстрируют экспериментальные данные, правые - расчетные осциллограммы. В расчетах параметры, определяющие граничный импульс, а также все прочие параметры смеси задавались так же как для рис. 2, кроме объемного содержания газа в области пузырькового кольца, которое составляло а ё0 = 0,01. На рис. (б) приведены расчетные и экспериментальные данные о затухании амплитуды первого всплеска давления в волновом пакете в системе с различным распределением пузырьков фреона. Из сравнения данных, приведенных на рис. 3, (а) и (б) следует, что, когда пузырьки расположены в кольцевом слое вблизи стенки (точки 3 и 4 на рис. 3, б) ), имеет место неплохое соответствие между экспериментальными и расчетными данными. При неравномерном по сечению распределении пузырьков (б) затухание происходит более интенсивно, чем в случае равномерного (или гомогенного) распределения. Но вместе с тем следует отметить, что для кольцевого расположения пузырьков, на начальных расстояниях (z » 0,25 м) амплитуда сформировавшегося всплеска давления значи-
тельно превышает соответствующие пиковые значения давления, возникающие при гомогенном распределении. Это связано с тем, что на начальном от границы участке формирования волны, происходит её поперечное преломление в зону пузырькового слоя, расположенную вблизи стенки жидкости и являющуюся более сжимаемой средой. На рис. 3, в) приведены поля скоростей (цифры у рисур соответствуют моментам времени). Поля скоростей построены на основе значений ^ и иг (для более наглядного представления поля скоростей представлены не для всей оси z). Видно, что при эволюции волны образуются поперечные течения. Более детальная картина динамики волны в трубе, с кольцевым расположением пузырьков, приведено на рис. 3, г) в виде эпюров давления от координат г и z. Здесь и в последующих рисунках распределения давления представлены для всего осевого сечения трубы (-Яс < г < Яс). Числа, около эпюр давления, соответствуют моментам времени. Из рис. 3, г) видно, что возмущения, распространяющиеся в виде предвестника вдоль осевой зоны по чистой жидкости, значительно опережают основную волну. В частности, к моменту времени г = 0,4 мр возмущения по «чистой» жидкости еле заметной амплитудой дошли до сечения с координатой z0 = 0,6 м, в то время как основная волна дошла лишь приблизительно к отметке с z0 = 0,04 м. Вследствие того, что акустический импеданс р0с для «чистой» жидкости, значительно выше, чем аналогичный параметр для пузырьковой смеси, возмущения, распространяющиеся вдоль центральной зоны «чистой» жидкости быстро съедаются из-за разгрузки на боковой границе с пузырьковой жидкостью. Эта граница (г0 = Яр) для возмущений, распространяющихся в центральной зоне, фактически играет роль свободной поверхности. А для основных возмущений, распространяющихся по газожидкостному кольцу, эта граница (г0 = Кр) действует аналогично твердой стенке. Кроме того, из-за сильного различия акустических импедансов «чистой» и пузырьковой жидкостей, профили давления по сечению трубы вблизи поршня в период его активного воздействия ( 0 < г < г1 + г2 0) сильно неоднородны. Из расчетных данных, приведенных на рис. 3, г), следует также, что в процессе эволюции определяющую роль играет газожидкостное кольцо. На фоне распространяющихся в пузырьковом кольце со скоростью ~ 160 - 180 м/р
волн давления, быстрые возмущения, распространяющиеся по зоне «чистой» жидкости, как уже отмечено, практически незаметны.
3.3. Пузырьковое ядро
На рис. 4 представлены результаты численных расчетов по динамике волн давления в системе, когда в приосевой зоне находится пузырьковая жидкость из азотных пузырьков, в пристенной - «чистая» жидкость (рис. 4, а ) ). На рисурах, отмеченных буквой б) приведены расчетные данные, иллюстрирующие затухание амплитуды первого всплеска давления в жидкости с различным распределением пузырьков по сечению. На рис. 4, г) приведены эпюры давления, соответствующие моментам времени 0,4 мс,
1,2 мс, 2 мс. Граничный импульс задавался как для предыдущего случая, представленного на рис. 3. Для величин параметров, определяющих воздействие поршнем, приняты следующие значения Д^ = 0,38 м/р, г1 = 0,25мр, г2 = 0,25 мр. При одинаковом воздействии жестким поршнем на гомогенную и расслоенную среду, датчик D1 с координатой z = 0,001м, расположенный на стенке, регистрирует импульсы в случае пузырькового ядра с амплитудой Др = 2,7 • 105 Па , в случае пузырькового кольца и гомогенного распределения с амплитудами - Др = 3,2 105 Па и Др = 2,1 • 105 Па соответственно. Причем (см. рис.
4, б) ) на начальном участке (Zо < 0,6 м) интенсивность затухания для неоднородной смеси сильнее, чем для гомогенной. Но, как следует из показаний более дальнего датчика (z0 > 0,8м), наоборот, затухание становится более интенсивным в случае гомогенного распределения. Характер затухания волны в этой ситуации такой же, как и в случае пузырькового кольца. Аналогичным образом происходит преломление волн в зону сильно сжимаемой пузырьковой среды.
На рис. 4, в) приведены поля скоростей (цифры у рисур соответствуют моментам времени). На рис. 4, г) представлены эпюры давления в различные моменты времени. Здесь также видно, что к моменту времени г = 0,4 мр слабые возмущения, распространяющиеся со скоростью звука в «чистой» жидкости, дошли до
отметки z0 = 0,6 м, в то время как основные волны давления в зоне пузырькового ядра только начинают формироваться. По эпюрам давления, соответствующим временам 1,2 мс и 2 мс (рис. 4, г) ) видно, что волна, образованная ударом жесткого поршня по всему сечению трубы, вдоль приосевой пузырьковой зоны движется со скоростью ~ 100 м/с, приблизительно равной скорости волны в соответствующей гомогенной пузырьковой смеси.
Как видно из рис. 3, г) и 4, г), соответствующих моменту времени 0,4мр , амплитуда давления в областях «чистой» и пузырьковой жидкостей в зоне контакта с твердым поршнем одинакова. Отметим, что для этого момента времени воздействие скоростью жесткого ударника продолжается. Этот факт является парадоксальным, потому что, как было отмечено выше, акустические импедансы «чистой» и пузырьковой жидкости отличаются более, чем в 15 раз, и при одинаковом воздействии скоростью ударника на границу z0 = 0 при продолжении воздействия амплитуда давления в областях «чистой» и пузырьковой жидкости должна быть значительно неоднородной. Для более детального анализа этой ситуации на рис. 5, в) приведены профили давления для сечения z = 0,001м при кольцевом распределении пузырьков фреона по сечению трубы с объёмным содержанием аЁ0 = 10-2 (штриховая и штрихпунктирная линии соответствуют случаям «чистой» и гомогенной пузырьковой жидкости (а Ё0 = 0,005 ) ). Параметры системы и первоначального импульса для рис. 5 такие же, как для рис. 3. Видно (сплошная линия), что на начальном этапе (г < 0,001 мр) воздействия жестким поршнем может наблюдаться сильно неоднородное «П»-образное распределение давления. Но постепенно, для больших времен, значительно превышающих времена прохождения волн давления расстояний порядка линейных масштабов поперечной неоднородности пузырьковой системы (в данном случае они определяются величинами Яс = 27 • 10-3 м и Яр = 19 • 10-3 м) происходит на поршне установление более однородного профиля давления. На рис. 5, б) представлены рас-
Список использованной литературы:
1. Нигматулин Р.И. Динамика многофазных сред. М.: Наука, 1987. Ч. 1. 464 с.; Ч 2. 360 с.
2. Самарский А.А, Попов Ю.П. Разностные методы решения задач газовой динамики. М.: Наука, 1980. 352 с.
3. Ильгамов М.А., Гильманов А.Н. Неотражающие условия на границах расчетной области. - М.: ФИЗМАТЛИТ, 2003. -240 с.
4. Донцов В.Е., Накоряков В.Е. Волны давления в газожидкостной среде с расслоенной структурой жидкость - пузырьковая
смесь // ПМТФ. 2003. Т.33. С. 102-107.
