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

УДК 532.5

раздел ФИЗИКА

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

© Е. В. Волкова1,2*, Э. Ш. Насибуллаева1’2

1 Центр «Микро- и наномасштабная динамика дисперсных систем»,

Башкирский государственными университет Россия, Республика Башкортостан, 450076 г. Уфа, ул. Заки Валиди, 32.

Тел./факс: +7 (347) 229 96 70.

E-mail: ekaterina. volkova@ bashedu. ru 2Институт механики им. Р. Р. Мавлютова Уфимского научного центра РАН Россия, Республика Башкортостан, 450054 г. Уфа, пр. Октября, 71.

Тел./факс: + 7 (347) 235 52 55.

E-mail: elvira. nasibullayeva @ bashedu. ru

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

Ключевые слова: акустическое поле, одиночныгй пузыгрек, направленная диффузия, численные методы .

Введение

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

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

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

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

Постановка задачи и основные уравнения

Рассматривается одиночный сферически-симметричный газовый пузырек, который совершает радиальные колебания под действием акустиче-

ского поля, изменяющегося по закону синуса: pa (t) = —APsin ®t (AP - амплитуда акустического

поля; ш - круговая частота).

Моделирование задачи основывается на следующих допущениях: теплообмен в жидкости отсутствует; изменение размеров пузырька за период колебания происходит только за счет диффузии; разница между скоростью жидкости на поверхности и скоростью изменения радиуса пузырька незначительно для уравнения Келлера-Миксиса; длина волны звукового поля много больше размера пузырька: шa << C (a - радиус пузырька; C - скорость звука в жидкости).

Скорость переноса массы газа mg через подвижную границу пузырька определяется через градиент массовой концентрации газа, растворенного в жидкости, с по следующей формуле:

2 т^дс

mg = 4na plD— ,

g l dr r==a

где p; - плотность жидкости; D - коэффициент диффузии; r - сферическая координата.

Полагая поведение пузырька политропным, получим, что динамика одиночного пузырька будет описываться следующим нелинейным дифференциальным уравнением Келлера-Миксиса:

1—— \aa +—| 1 —

=| 1+-

2 V ЗС P a dP

I------+-------------’

С) Pi PiC dt

в котором

P =

Po +:

2o

-*o У

\-3У

2o 4pa

- Ро - Ра (?), а а

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

m

a

m

v g o У

v ao У

* автор, ответственный за переписку

сти. Нижний индекс «0» обозначает значение параметра в начальный момент времени. Начальные условия для уравнения (1) имеют вид:

а(0) = а0, а(0) = 0.

Уравнение диффузии растворенного в жидкости газа в сферических координатах имеет следующий вид:

дс а2а дс = ^ 1 д ^ 2 дсЛ

дг г2 дг

дг

■ + -

г

г

дг

2 • / 2

где а а / Г - радиальное поле скоростей.

Граничное условие на поверхности пузырька вычисляется по закону Г енри:

(

С

= Са = НР;

* 0

т„

У У3г а

V т, 0у

V ао у

где ^ - давление газа; Н -константа Г енри.

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

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

вид: Л = С = с Из уравнения диффузии

I г=0 I г=~ ~'

(2) вычисляется поток массы из пузырька в жид-.дс

кость

І = ~Рі°Т

дг

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

д = г3 - а3. Обозначая ~ = с - См, ~а = Са - С, и

вводя новую переменную для обозначения времени

г

т(г) = 9В$а\г')йг\ т(0) = 0 уравнение диффу-

0

зии (2) вместе с начальными и граничными условиями примут вид:

д~ =_э_

дт дд

1 +

д

.4/3

а

э~

дд

Р = 0, Р 0 = Ра,

1д=^ 1д=0 а

~1 0 = 0,

1т=0

(3)

где а и Са рассматриваются как функции времени т.

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

задачи на две вводится следующее осреднение по периоду акустической волны Т = 1/ю:

1 т(Т)

1 {('г) т=Т(Т) 11 {('Т)ЛТ =

ТТ

= | /(д, г)а\г)йг | а4(г)ёг.

0 10 Тогда граничное условие у стенки пузырька запишется в следующем виде:

и'

~д=0 = НР,°{р,)т - с~ + НР,0Р -(Р,}ті

р, =

_ Р, _

0

т а

1 т0 а о

и задача с осциллирующей частью граничного условия примет следующий вид:

дс_ д

дт дд

Ґ \4/Ь~

і+4] др

а у1

дд

со следующими граничными и начальными условиями:

^Сі д=0 = НР,

-3у

■(р,-*

С~ОЭС |т=0 СОЭС |д=„ 0.

Задача с постоянной частью граничного условия будет иметь следующий вид:

дс_ д

дт дд

і + ^_V" дрт а3) дд

с граничными и начальными условиями

эт д=0

-37

С

с = с

эт |х=0 эт |д=ет

= 0.

Сумма решений будет решением для полной задачи, т.е.

т) = ^т)+^т).

В дальнейшем в данной работе будет рассматриваться решение только задачи с осциллирующей частью граничного условия.

Численная реализация

Численное исследование уравнения Келлера-Миксиса для одиночного пузырька (1) осуществлялось с помощью метода Дормана-Принца восьмого порядка точности [5]. Численный метод верифицирован с использованием встроенной процедуры МаЙаЪ оёе45 для решения обыкновенных дифференциальных уравнений. Вычисления усложнились после включения в задачу динамики одиночного пузырька расчета массопереноса через стенку пузырька.

Решение диффузионной задачи экспоненциально стремится к нулю при приближении к бесконечной границе, когда д ^ . Таким образом,

граница может быть установлена для некоторой конечной величины д = дс .

Г = а

г=а

Т

После всех преобразований, задача диффузии (3) сводится к виду

ди д

дт дд

. ди

к(д, т)—

дд_

и 0 = 0,

т=0

'д=дс

= 0,

(4)

и n= ua,

1д=0 й

и решается по схеме

,4/3

где к(д,т) = 11 + ^3

Кранка-Николсон [6].

Дискретизация уравнения (4) по пространству и времени заключается в следующем уравнении для определения решения на п-м временном шаге методом итераций:

( \ ( 1. \

хп+1 =

К

/ + \ в,

v 2hg у

Xn +

(5)

7 - 2^ в-

ч 2Пд ;

+ }Ь^(¥п+1 + Рп )

где Ьт и и - шаги по времени и пространству; 7 -т д

единичная матрица; Вп ~ симметричная трехдиагональная матрица размера (М - 1)х(М - 1), состоящая из коэффициентов к; М - размер сетки по пространству; X - вектор решения и ¥ - силовой вектор, зависящий от граничного условия. Численные расчеты и обсуждение результатов Численный метод, предложенный в данной статье, был сопоставлен с автомодельным аналитическим решением [7], где зависимость радиуса от времени имеет вид а = ^/2^г + а1 (в берется пропорционально Б). При замене П = г / а (г) в уравнении (2) аналитическое решение для с примет следующий вид:

J 1

U (n) = -C J-exp

n n

u| = u ,

rn=1 a

J 1

U c = - C J-eXP

1 n

P_

D

^n2 1 v — + -v 2 n у

C = - U,

J -exp

1 n

' р_ f n2 +- у

D V, 2 n )_

Было получено хорошее совпадение численных результатов с аналитическими. После некоторого переходного периода при установлении устойчивого профиля концентрации относительная погрешность решения не превышала 10-2.

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

При расчетах рассматривался пузырек воздуха в воде при температуре 20°С и следующих значениях параметров: у = 1.4, а = 0.0725 Н/м,

Р0 = 105 Па, С1 = 1500 м/с, Ш = 2п-20 кГц, В = 2-109 м2/с и с0 = 2.5-10-5.

На рис. 1 показано сравнение с результатами, полученными в работе [4] в рамках асимптотической теории с колебаниями малой амплитуды. Безразмерная концентрация газа с* = с / сх, зависящая от безразмерной эйлеровой координаты г* = г / а0, показана в различные моменты времени для первого периода колебания пузырька для числа Пекле Ре = а^Ю / В = 103 и радиуса пузырька, рассчитанного по формуле а(г) = 1 + 0Л8т(2п£). Поток массы через границу пузырька не учитывается. Пунктирной линией обозначены результаты расчетов, проведенных с использованием схемы Кранка-Николсон, сплошной линией показаны данные, полученные в работе [4]. Искривленная и вертикальная пунктирная линии показывают положения стенки пузырька, полученные в [4] и в данной работе, соответственно. Как можно заметить, результаты хорошо согласуются. Различия могут быть объяснены разной точностью используемых численных методик.

(а) (Ь) (с) (і)

Рис. 1. Безразмерная концентрация газа с* в зависимости от безразмерной эйлеровой координаты г* для Ре = 103, а(г) = 1 + 0.Ып(2яг) в моменты времени (а) г = Г/4; (Ь) г = Г/2; (с) г = 3Г/4 и (і) г = Г Результаты получены в работе [4]

(пунктирная линия) и в настоящей работе (сплошная линия).

(a)

Ф)

(с)

(d)

t/T

(e)

Рис. 2. Сравнение асимптотических профилей концентрации для моментов времени (а) г = Т/4; (Ь) г = Т/2; (с) г = 3Т/4 и (й) г = Т, а также радиусов пузырька (е) для первого периода акустического поля. Сплошная линия соответствует случаю а(г) = 1 + 0.Ып(2лг) и Ре = 104; пунктирная и штриховая линии - а(г), рассчитанному по уравнению (1) без учета эффекта изменения массы и с учетом соответственно для а0 = 1 мкм и АР = 0.65 атм.

На рис. 2 показаны профили концентрации газа у стенки пузырька на первом периоде колебания для радиуса пузырька, заданного формулой a(t) = 1 + 0.1sin(2nt) (сплошная линия); рассчитанного из нелинейного уравнения Келлера-Миксиса (1) без изменения массы в функции давления газа pg (пунктирная линия) и с изменением массы (штриховая линия). Расчеты для линейных колебаний пузырька были проведены при Pe = 104, что соответствует начальному радиусу пузырька a0 = 1 мкм и амплитуде давления AP = 0.65 атм и наиболее приближено к функции a(t) = 1 + 0.1sin(2nt) для первого полупериода. Однако гармоническое приближение отклоняется от вычисленной кривой a(t) на втором полупериоде (рис. 2d), что ведет к видимым различиям профилей концентрации для

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

Далее было проведено сравнение с приближенной теорией колебаний пузырька большой амплитуды для больших чисел Пекле [8]. На рис. 3 представлен нормированный радиус одиночного газового пузырька атах/а0 и осредненная концентрация газа около стенки пузырька < с >т для различных амплитуд давления акустического поля АР. Рис. 3(Ь) позволяет оценить решение, полученное представленным в данной работе численным методом с приближенным решением, использованным в статье [8] (см. рис. 3(а)).

(а) (Ь)

Рис. 3. (а, Ь) Нормированный радиус пузырька атах/а0; (с, й) осредненная концентрация газа около стенки пузырька <с>т в зависимости от начального радиуса а0 для различных амплитуд давления др = 1.1 ^1.5 Па. Результаты получены (а, с) в работе [8] и (Ь, й) в настоящей работе. Сплошные линии - для кривых с учетом изменения массы газа на 20-м периоде колебания пузырька; штриховые линии - для кривых без учета изменения массы газа.

(а) (Ь) (с)

Рис. 4. (а) Максимальный радиус пузырька атах и (Ь) максимальная масса ттах в зависимости от периода колебаний Т; (с) функция нормированного радиуса пузырька а/а0 в промежутке от 4676 до 4680 периодов.

Видно, что характер кривых для осредненной концентрации газа, представленных в данной работе, отличается от кривых в [8], особенно для малых начальных радиусов. Действительно, на рис. Зі средняя концентрация увеличивается при малых а0, в то время как на рис. Зс она растет для тех же кривых. Кроме того, значения средней концентрации в первом и втором случаях различаются на один порядок. Изменение массы начинает влиять на систему при амплитудах акустического поля ДР = 1.4 атм и выше (см. рис. 3(і)).

На рис. 4 показаны максимальные радиус пузырька (рис. 4а) и его масса (рис. 4Ь) за 6600 периодов при ю = 20 кГц, АР = 1.5 атм для насыщенной система вода-воздух при 1 атм и 293 К. Из рис. 4с видно, что режим колебаний меняется на некотором переходном периоде Г = 4679, когда текущий радиус достигает примерно 1.186 мкм. Связано это с тем, что при учете в уравнении для давления газа изменения массы радиус пузырька будет расти вследствие направленной диффузии. При достижении радиусом некоторого порогового значения (см. рис. 3а,Ь при а0 ~ 1.2 мкм), его амплитуда колебания резко увеличивается, поскольку внешнее давление преодолевает давление поверхностного натяжения Лапласа. Этот эффект называется «Гигантский отклик» [8]. При этом скорость направленной диффузии сильно увеличивается.

Выводы

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

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

Авторы выражают благодарность за помощь при постановке задачи и ценные замечания И. Ш. Ахатову и Н. А. Г умерову. Работа выполнена при финансовой поддержке грантов Министерства образования и науки РФ (11.G34.31.0040), соглашения №8861 и РФФИ (проект №11-08-00823-а).

ЛИТЕРАТУРА

1. Hsieh D. Y., Plesset M. S. Theory of rectified diffusion of mass into gas bubbles // J. Acoust. Soc. Am. 1961. V. 33. №2. P. 206-215.

2. Eller A., Flynn H. G. Rectified diffusion during nonlinear pulsations of cavitation bubbles // J. Acoust. Soc. Am. 1965. V. 37. №3. P. 493-503.

3. Eller A. Growth of bubbles by rectified diffusion // J. Acoust. Soc. Am. 1969. V. 46. №5. P. 1246-1250.

4. Fyrillas M. M., Szeri A. J. Dissolution or growth of soluble spherical oscillating bubbles // J. Fluid Mech. 1994. V. 277. P. 381-407.

5. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных

дифференциальных уравнений. Нежесткие задачи.

М.: Мир, 1990. 512 с.

6. Crank J., Nicolson P. A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type // Proc. Camb. Phil. Soc. 1947. V. 43 №1. P. 50-67.

7. Scriven L. E. On the dynamics of phase growth // Chem. Eng. Sci. 1959. V. 10. P. 1-13.

8. Akhatov I., Gumerov N., Ohl C.-D., Parlitz U., Lauterborn W. The role of surface tension in stable single-bubble sonolumi-nescence // Ph. Rev. Lett. 1997. V. 78. №2. P. 227-230.

Поступила в редакцию 07.09.2012 г.