CC BY
53
Вычислительные технологии
Том 15, № 1, 2010
Метод численного решения задач сильного сжатия несферического кавитационного пузырька*
A.A. Аганин, Т. Ф. Халитова, H.A. Хисматуллина Институт механики и машиностроения Казанского научного центра РАН, Россия
e-mail: [email protected]
В предлагаемом методе движение пара в пузырьке и окружающей жидкости описывается уравнениями динамики сжимаемой теплопроводной жидкости с учетом фазовых переходов. Используются подвижные сетки, UNO-модификация метода Годунова второго порядка, явно-неявная схема.
Ключевые слова: динамика сжимаемой теплопроводной жидкости, кавитацион-ный пузырек, ударная волна, UNO-схема, метод Годунова.
Введение
В настоящее время сильное сжатие парогазовых пузырьков в жидкости привлекает интерес в связи с открытием таких явлений как однопузырьковая и многопузырьковая сополюмипесцепция [1-3], нейтронная эмиссия при акустической кавитации дей-терированного ацетона [4, 5]. Теоретические и экспериментальные исследования [6-11] привели к формированию представления о том, что световое излучение и нейтронная эмиссия осуществляются за счет сильного высокоскоростного сжатия пузырька. В финальной стадии такого сжатия в пузырьке формируется ударная волна, сходящаяся к его центру [6, 7]. Со временем ее интенсивность возрастает так, что кратковременно в центре пузырька образуется сферическое ядро с очень высокими значениями температуры и плотности [9, 11]. Согласно численным исследованиям сонолюминесценции [9], в области, размеры которой составляют ~ 10-9 м, в течение ~ 10-11 — 10-10 с сохраняются температуры ~ 106 К и плотности ~ 10 г/ем3, В случае нейтронной эмиссии [11] в области размером ~ 10-7 м в течение долей пикосекунд температуры имеют значения ~ 108 К, а плотности ~ 10 г/см3.
Малость размеров пузырька и горячего ядра в нем, а также короткая продолжительность происходящих процессов затрудняют экспериментальное изучение финальной стадии сжатия пузырька. В связи с этим важную роль в понимании происходящих при этом процессов приобретает математическое моделирование. В частности, указанное выше представление о механизмах сонолюминесценции и ядерного излучения сформировалось в результате вычислительных экспериментов с применением сферически симметричных моделей, в которых для описания движения газа в полости пузырька используются уравнения газовой динамики [4, 6, 7, 10-13].
В данных экспериментах одномерные уравнения газовой динамики решаются методами первого порядка точности Лакса — Фридрихса [6], Годунова [4, 11, 12], комплексом
*Работа выполнена в рамках программы ОЭММПУ РАН и при финансовой поддержке РФФИ (грант № 08-01-00215 и 08-01-97029).
© ИВТ СО РАН, 2010.
программ 1)У.\Л21) [7, 9, 10]. При этом в [7, 9, 10] численное моделирование на основе уравнений газовой динамики в области газа и жидкости используется как в ходе сжатия пузырька, так и в ходе его расширения. Применение такого подхода требует больших затрат компьютерного времени. В связи с этим в работах [11, 14, 15] предложен более экономичный метод расчета, в котором при расширении и на начальной стадии сжатия пузырька, когда его поверхность движется с малыми числами Маха, газ полагается гомобарическим, а жидкость малоежимаемой. Полная гидродинамическая модель используется лишь на высокоскоростной стадии сжатия, когда число Маха движения стенки пузырька становится близким к единице.
Указанный выше ударно-волновой механизм достижения высоких значений температуры и плотности вещества в полости пузырька существенным образом опирается на гипотезу сферической симметрии сжатия. Однако в реальности пузырек при сжатии чисто сферическим не бывает. Отклонения от сферической формы могут привести к значительному снижению экстремальных значений плотности и температуры, что, например, и происходит при лазерном сжатии сферических мишеней [16-18]. Поэтому изучение влияния несферичности пузырька на основные параметры происходящих при сильном сжатии процессов является весьма актуальным.
В тех фазах сильного сжатия, где искажения сферической поверхности пузырька остаются малыми, изучение их влияния можно провести существующими методами, основанными на расщеплении движения газа в пузырьке и окружающей жидкости на сферическую составляющую и ее малое несферическое возмущение [19-31]. При этом поверхность пузырька представляется в виде суммы поверхностных сферических гармоник с коэффициентами, характеризующими амплитуду отклонения поверхности пузырька от сферической в виде соответствующих гармоник.
Если по искажениям сферичности задача ставится как линейная, что соответствует случаю очень малых искажений, то полученный при расчете сферической составляющей закон изменения радиуса пузырька служит в качестве входных данных для решения уравнений эволюции искажения его сферичности. При описании эволюции искажения сферичности пузырька жидкость обычно предполагается несжимаемой, а плотность содержимого пузырька однородной [19, 25, 32-39]. В результате этого отклонение от сферической формы описывается обыкновенным дифференциальным уравнением второго порядка. Методы исследования малых искажений сферической формы пузырька основаны на предположении о том, что обратным влиянием искажений сферичности на радиальную составляющую движения газа в пузырьке и окружающей жидкости можно пренебречь. Вместе с тем в финальной высокоскоростной стадии сильного сжатия сферическая и несферическая составляющие становятся взаимозависимыми как в силу больших скоростей сжатия, так и в результате роста искажений сферичности. Поэтому в финальной стадии сжатия указанные методы оказываются неприемлемыми.
Изучение динамики пузырька при значительных деформациях его поверхности обычно ведется методом граничных элементов. Этот метод с успехом применяется в задачах с большими формоизменениями пузырька, находящегося как вдали от разнообразных границ [40], так и вблизи жестких стенок [41-44], свободных поверхностей [45], других пузырьков [46]. Основным допущением метода граничных элементов является предположение о несжимаемости жидкости и потенциальности ее движения. На высокоскоростной стадии сжатия такие допущения несправедливы.
Представляется, что наиболее подходящим методом исследования влияния несферичности пузырька при сильном сжатии является прямое численное моделирование на
основе уравнений газовой динамики для описания движения как газа в пузырьке, так и окружающей жидкости. Таким методом можно изучать сжатие при малых и больших деформациях поверхности и скоростях сжатия. Успех применения прямого численного моделирования во многом зависит от экономичности используемого при этом метода решения. Обзоры методов расчета нестационарных многомерных газодинамических задач можно найти, например, в [47-50],
В настоящей работе предлагается метод численного решения задач сильного высокоскоростного сжатия осесимметричного кавитационного пузырька. Этот метод является по существу развитием предложенного ранее для расчета аналогичных одномерных [51] и двумерных [52] задач сжатия пузырька без учета теплопроводности газа и жидкости и испарения-конденсации на поверхности пузырька.
Для описания движения газа и жидкости применяются двумерные уравнения динамики газа и сжимаемой жидкости с учетом эффектов теплопроводности, испарения и конденсации на межфазной поверхности. Их решение находится в два этапа на основе принципа расщепления по физическим процессам. На первом этапе, подробно изложенном в [52], рассматриваются уравнения газовой динамики без учета теплопроводности газа и жидкости, испарения и конденсации на межфазной поверхности. Уравнения состояния принимаются в виде довольно общих зависимостей давления и внутренней энергии от плотности и температуры. Решение находится с помощью модификации метода Годунова на основе UNO-схемы (UNO — Uniformly NonOseillatorv, равномерно неосцпллнрующая) второго порядка точности по пространству и времени [53], На втором этапе, одномерный аналог которого применяется в [51], решается двумерное уравнение теплопроводности с граничными условиями на поверхности пузырька, учитывающими эффекты испарения и конденсации по формуле Герца — Кнудеена —Ленгмюра [11], При этом используются схемы со вторым порядком точности по пространству: неявная схема переменных направлений и явная схема.
Для иллюстрации работоспособности предлагаемого подхода рассматривается задача сильного сжатия кавитационного пузырька с малыми начальными отклонениями его формы от сферической. Для демонстрации влияния теплопроводности и фазовых переходов применяется решение аналогичной задачи без учета этих эффектов. Результаты тестирования работы предлагаемого метода и оценка его экономичности при решении двумерных задач динамики пузырька изложены в [52],
1. Постановка задачи
В данном разделе задача динамики пузырька сначала формулируется в компактной векторной форме. Затем она приводится в произвольных подвижных координатах относительно неподвижных сферических. По мнению авторов, при таком выборе эйлеровой системы отсчета расчет задач с относительно небольшой несферичностью пузырька будет наиболее эффективным,
1.1. Основные уравнения и граничные условия
Рассматривается динамика осесимметричного пузырька в слое жидкости. Форма поверхностей пузырька и внешней границы жидкости в сферических координатах (r, в) однозначно описывается зависимостями r = rs(e,t) и r = г/(e,t), Здееь r — радиальная
в
Движение содержимого пузырька (неконденеируемого газа или пара) и окружающей жидкости описывается системой уравнений
¿р „ ¿и „ ¿Е „ . „ , . .
-Р+рУ- и = 0, р— + Ур = 0, р— + У- ри-кУТ =0 1)
оЬ оЬ оЬ
с замыкающими соотношениями (уравнениями состояния) вида р = р(р, Т), е = е(р, Т). Здесь Ь — время, р — плотность, и — скорость, р — давление, Е = е + и2/2 — удельная полная энергия, е — удельная внутренняя энергия, Т — температура, к — коэффициент теплопроводности,
Граничные условия на межфазной поверхности г = г3 имеют вид
р+ (Б - и+) ■ п = р- (Б - и-) ■ п = р+ = р- + 2Ио, (2)
Т+ = Т-, (кУТ ■ п)+ - (кУТ ■ п)- = Ц, (3)
где Б = О ■ п — скорость смещения элемента поверхности пузырька; п — внешняя единичная нормаль; о — коэффициент поверхностного натяжения; 2И — средняя кривизна поверхности; — интенсивность фазовых преобразований; I — теплота парообразования, Знак плюс в верхнем индексе означает отношение к стороне жидкости, минус — к стороне пара.
Для пузырька, заполненного неконденсируемым газом, в (2), (3) нужно положить ] = 0. Если пузырек паровой, то интенсивность фазовых преобразований ] определяется по формуле Герца — Кнудсена — Ленгмюра [И]
у -Ш7ГЫГ ] (4)
где
Х = е"п2 [ е~х2с1х \ , П= у/Е^Т-.
V к ) V2р-
Здесь р$ — давление насыщения, аас — коэффициент аккомодации, Ки — газовая постоянная для пара. Посредством параметра х учитывается подвижность поверхности пузырька. Зависимости р^ (Т), к±(Т), 1(Т), о(Т) представляют собой аппроксимации экспериментальных данных.
На внешней поверхности жидкости г = Гf обычно задаются давление и температура
р = Pf (Ь),Т = Tf (Ь). (5)
Вместо динамического условия в (5) можно использовать кинематическое или более общее импеданеное соотношение [54].
1.2. Формулировка задачи в подвижных координатах
В реальных условиях при сильном сжатии пузырька форма его поверхности обычно отклоняется от сферической. В окрестности межфазной границы в жидкости возможны большие градиенты давления, в жидкости и паре — тонкие тепловые пограничные слои. Градиенты давления влияют на скорость движения поверхности пузырька, а от тепловых пограничных слоев зависит масса пара в пузырьке, а значит и изменение радиуса пузырька в финальной стадии сжатия. К тому же при сильном сжатии в полости
пузырька могут возникать очень интенсивные радиально сходящиеся волны сжатия, в том числе ударные. Они определяют экстремальные значения газодинамических параметров в окрестности центра пузырька. Методы, применяемые для расчета таких задач, должны корректно описывать все перечисленные особенности процесса сильного сжатия, С этой целью в настоящей работе применяется смешанная эйлерово-лагранжева система (СЭЛ) координат, связанная с поверхностью пузырька. При этом предполагается, что поверхность пузырька при сжатии может значительно отклоняться от сферической, но без нарушения однозначности вдоль радиальной координаты (т, е, без нарушения однозначности функции r = rs(9,t)), СЭЛ-координаты не только позволяют учесть перечисленные особенности задачи, но и дают возможность плавно переходить от сетки с хорошим разрешением пограничных слоев к сетке с хорошим разрешением ударной волны в газе в финальной стадии процесса.
Соотношения между эйлеровыми координатами r, 9 и связанным с ними временем t и СЭЛ-координатами п и связанным с ними временем т имеют следующий вид:
r = r(£,n,T), 9 = 9(П^ t = т.
Предполагается, что в координатах п линии, соответствующие поверхности пузырька и внешней границе жидкости, всегда остаются координатными линиями £ = const,
В СЭЛ-координатах уравнения динамики газа и сжимаемой жидкости можно записать как
где
q
Q
( р \
ри pv \рЕ J
q,
Qt + F + Gn = S,
F = Vhi, g = Vhg, s = Vhs, (
(6)
p (U - Uw) pu (U - Uw) + pCr pv (U - Uw)+ pr-1Ce \ pE (U - Uw) + p U - K(TrCr + Ter-2Ce)
g
( pV
puV pvV + pr-l,qO
\
\ pEV + pV - кТвr-2ne )
Здесь U = u£r + vr-1&, V = vr-l,qe; Uw = uw^r\ uw = rT; h = J2g g = r4 sin2 6; J = r^6n; u, v — компоненты вектора скорости частицы среды; uw — радиальная компонента вектора скорости подвижных координат.
Контактные условия (2), (3) в СЭЛ-координатах принимают следующий вид:
( 0 \
(2p + pv2) r-1
-pr-1uv + pr-1 ctg 9 0
D = Uo+ + j/p+ = Uo - + j/p-, p+ = p~ + 2Ha, T+ = T-, [к (TrCO + Te$r-2)]+ - [к (TrCO + TeCOOr-2)}- = jl
а условия (5) — следующий:
p = pf (t), T = Tf (t),
(7)
(8)
(9)
где
/° = //VeT+^v32
2. Метод численного решения
При решении уравнений (6) е граничными условиями (7)-(9) на каждом временном шаге сначала строится расчетная сетка, а затем в два этапа выполняются вычисления, построенные по принципу расщепления по физическим процессам. На первом этапе игнорируются теплопроводности пара и жидкости, фазовые переходы на поверхности пузырька. Уравнения (6) решаются без учета тепловых потоков с соответствующими граничными условиями. Влияние теплопроводности и фазовых переходов учитывается на втором этапе. Для этого поле полной энергии, вычисленное на первом этапе, корректируется с учетом тепловых потоков.
2.1. Построение расчетной сетки
Расчетная область в координатах г, в представляет собой криволинейный четырехугольник, состоящий из двух зон: газа и жидкости, В случае, когда зона жидкости является неограниченной, в качестве внешней границы расчетной области принимается достаточно далеко удаленная искусственная граница г = га(в,Ь). Ее начальное положение зависит от конкретной задачи, а ее перемещение находится из граничных условий, которые ставятся так, чтобы удовлетворить условиям на бесконечности. Если движение жидкости и газа симметрично относительно плоскости в = п/2 (ниже для простоты рассматривается только этот случай), то расчетная область представляет собой следующее множество точек:
{г, в; 0 < г < гf, 0 < в < п/2} .
В СЭЛ-координатах £, п эта область имеет форму прямоугольника. Она покрывается сеткой из (М9 + N1) х N0 ячеек, где Мд и N — количество ячеек в паправлении £ в зонах газа и жидкости соответственно, N0 — количество ячеек в направлении п-Сеточный аналог произвольной функции /обозначается как /, В дальнейшем также используются обозначения: нижние индексы г, к указывают на отношение к центру ячейки (г = 0,...,^ + N1 — нумерация ячеек в напр авлении £; к = 0, ...,N0 — нумерация в направлении п); г ± 1/2, к и г, к ± 1/2 — к гран ям £ = £^±1/2 и П = Пк±\/2 ячейки; г ± 1/2, к ± 1/2 — к вершинам ячейки; верхние индексы п и п +1/2 показывают отношение к моментам времени т = тп и т = тп+1/2 = (тп + тп+1) /2, Напри-
г, к + 1
т = тп будет обозначаться РПк+р Сферические координаты узлов ячейки гк имеют вид
г1±1/2,к±1/2 = г(£г±1/2, пk±1/2, ^ вк±1/2 = в(Пк+1/2)-
Для адекватного описания особенностей решения в начале сжатия (тонких тепловых пограничных слоев и градиентов давления) используется сгущение ячеек в радиальном направлении (по геометрической прогрессии) к межфазной поверхности как со стороны газа, так и со стороны жидкости, В направлении в размер ячеек одинаков и равен Дв, В СЭЛ-координатах сетка всегда остается равномерной с шагом △£ по ^направлению и Дп по п-направлению. Это значительно упрощает построение разностной схемы. Сгущение (или разрежение) сетки по геометрической прогрессии определяется ее начальным элементом и знаменателем, У прогрессий, используемых для неравномерного размещения узлов сетки на луче в = вк+1/2, начальными элементами являются длины прилежащих к поверхности пузырька отрезков в областях газа
г«(вк+1/2,Ь) , ,
Гмд+1/2,к+1/2 ~ Гмд-1/2,к+1/2 = РКЧ-Тт--(1и)
1Уд
и жидкости
гМд +3/2,к+1/2 — гМд+1/2,к+1/2 — т(^)
^№+1/2,^
N
(Н)
Знаменатели этих прогрессий определяются по известным значениям сумм Мд и N их элементов (первая сумма равна г3(0к+1/2, ¿), вторая — г/(0к+1/2,£) — Гц(6к+1/2,1)).
Функции в(^) и 7(¿) можно выбирать в зависимости от конкретной задачи. Если в(¿) < 1, то сетка в газе сгущается к поверхности пузырька, если в(¿) > 1 — в противоположную сторону. При в(^) — 1 сетка равномерная. Аналогично при 7(¿) < 7* сетка в газе сгущается к поверхности пузырька, при > 7* _ в противоположную сторону. При 7(¿) — 7* сетка равномерная. Здесь 7* — (г/ — г3)^/(гsN,1), При решении рассматриваемых ниже задач функции в(¿) и 7^) принимаются кусочно-линейными в виде
ф(1)
Фъ Ф1
ф2,
Ф2 — Ф
¿2 — *
1 (* " ¿1)
* < ¿1,
¿1 < * < ¿2,
¿2 < *,
(12)
где ¿1 — ¿а, ¿2 — ¿в для ф — в, и ¿1 — Ьс, ¿2 — ¿о для ф — 7,
Параметры вь в2> ¿а функции в(^) выбираются с учетом того, что в начале сжатия нужно добиться хорошего разрешения тонкого теплового пограничного слоя в окрестности пузырька, а в конце сжатия — сходящейся к центру пузырька ударной волны, формирующейся в его полости. Поэтому указанные параметры выбираются так, чтобы сначала сетка сгущалась к поверхности пузырька, а затем плавно переходила в равномерную (резкое изменение сетки может приводить к появлению нефизических осцилляций). Для хорошего разрешения ударной волны вблизи центра пузырька необходимо использовать сетку со сгущением к центру. Однако такой цели в настоящей работе не ставилось, а потому сетка со сгущением к центру в рассматриваемых ниже примерах не применяется. Параметры функции 7^) (7^ 72, ¿с, ¿ц) выбираются исходя из того, что в жидкости, как и в газе, вблизи поверхности пузырька возникает тонкий тепловой пограничный слой. Для его хорошего разрешения используется сетка со сгущением к поверхности пузырька. Кроме того, расчеты показали, что лучше выбирать 10 < в(¿)-
Описанный способ построения сетки позволяет эффективно изучать как эволюцию формы пузырька, так и динамику газа внутри пузырька в процессе сильного сжатия.
2.2. Алгоритм расчета первого этапа
На первом этапе влияние теплопроводности не учитывается. Решаемые при этом уравнения можно было бы получить непосредственно из уравнений (6)-(9) при к — 0, Однако расчеты оказываются более экономичными, если уравнения представить в виде
+ е + — Б,
(13)
где Ц — Зц, Е — Л, Ъ — 3%, Б — Зб,
( Р (и — и*)
ри (и — и,ш) + р£г
рь (и — и*,) + рг-1£ рЕ (и — и*)+ ри
рУ
риУ рьУ + рг-1'чв \ рЕУ + рУ )
/ —р (2u + v ctg 0) r-1 \ _ —р (2u2 + uv ctg в — v2) r-1 —р (3uv + v2 ctg в) r-1 V — (P + PE) (2u + v ctg в) r-1 /
с граничными условиями (7) и p = pf (т) из (9),
На каждом временном шаге система (13) с указанными граничными условиями решается численно с применением модификации метода Годунова на основе UNO-ехемы второго порядка точности по пространству и времени [53] (эту схему называют также схемой предиктор-корректор Ханкока [50]), Алгоритм решения здесь не приводится, так как он подробно описан в [52],
После выполнения первого этапа становятся известными поля плотности и вектора количества движения в момент времени т = тп+1, а также промежуточное поле полной энергии. На втором этапе последнее корректируется с учетом теплопроводности.
2.3. Алгоритм расчета второго этапа
На втором этапе для учета влияния теплопроводности решается уравнение энергии, С этой целью в предлагаемом методе применяется явно-неявная схема, в которой решение уравнения энергии находится по явной схеме с использованием промежуточного поля полной энергии, полученного на первом этапе, и вспомогательного поля температуры, рассчитанного по явно-неявной схеме из уравнения теплопроводности. Последнее получается из уравнения энергии и имеет вид
Тт + T (U — Uw)+ Tn V +
Thdi
nVh (Tss2 + Tarier 2
Vhdr]
+ Tn Цв) Цв r
А
Vh
(14)
где к\ = — 1/ (р£т), к2 = (р — р2£р) / (рвт), з = л/£г2 + £в2г~2, Уравнение (14) решается с контактными условиями (8), условием на внешней поверхности жидкости Т = Tf (т) из (9) и дополнительным условием Т = 0 при г = 0,
Вспомогательное поле температуры Т^к определяется численно из уравнения (14) с указанными граничными условиями по явно-неявной схеме. Сначала применяется неявная схема переменных направлений, в которой по каждому из направлений £ и п решается система разностных уравнений: одна система по £-направлению для вычисления Т,(к
гр rpn
J-ik J-ik Атп
— \ (U — Uw)ТCr +
+
ki д
K(T)Vh (t^s2 + Т^вщг-2
VhcX
вторая — по ц-паправлепию для расчета Tik
(15)
Тгк — Т
ik
Ат n
— \ vr-1 (ТСв + TnЦв) +
+
k1 д
Т?Св + Т?цв) Цвr
-2
+
(16)
Затем по явной схеме вычисляется Т^к:
= {(Й - й*Шг + УГ-1 (Т& + 7» + [к(Г)>/Ь (т?52 + Т^ет-2)
к1 9 (Т& + пет-2] + 4 + ^
(17)
у/Кдч ■ ] ■ ^
В этих схемах координаты ячеек сетки относятся к моменту времени т — тп+1. Все остальные параметры, входящие в правую часть уравнений (15)—(17), кроме помеченных значками ~ и —, относятся к моменту времени т — тп.
При каждом значении к схема (15) содержит пять видов разностных уравнений: вблизи центра пузырька, внутри зон (газа и жидкости), вблизи межфазной границы со стороны газа (граница справа), вблизи межфазной границы со стороны жидкости (граница слева) и вблизи внешней границы, которые различаются способом аппроксимации пространственных производных по ^-направлению,
В ячейках, прилегающих к центру пузырька и к внешней границе расчетной области, а также внутри зон для аппроксимации пространственных производных используются центральные разности второго порядка. При этом уравнения, относящиеся к приграничным ячейкам, содержат по два неизвестных значения (Т1к.; Т2к) и (ГГмд+м{-1,к, Тмд+м1,к), а в каждое уравнение внутри зон входят неизвестные значения Тг-1к.; Т^к, Тг+1к в трех последовательных по ^-направлению ячейках сетки,
В разностных уравнениях, относящихся к ячейкам, примыкающим к межфазной границе слева и справа, используются односторонние аппроксимации второго порядка, В эти аппроксимации входят неизвестные значения температуры Т*к в точках межфазной границы. Они выражаются через значения температуры во внутренних узлах расчетной сетки следующим образом:
%к 8(с+ - с") ' (18)
где с — кз, й — к(£вщг-2)/з, / — ]I — ^¿+Т+ — . Соотношение (18) получается
из краевых условий (8), второе из которых используется в виде
с+Т+ — с-Т~ — /.
Разностные уравнения в ячейках, примыкающих к межфазной границе слева и справа, получаются подстановкой выражения (18) температуры Т*к в соответствующие односторонние аппроксимации производных. Они содержат по четыре неизвестных значения температуры: по два с каждой стороны межфазной границы.
Таким образом, для каждого к схема (15) приводит к системе Ng + N1 линейных алгебраических уравнений относительно Тгк с пятидиагональной матрицей,
В направлении п нет внутренних границ, а на внешних границах ставятся условия симметрии. Поэтому для каждого г схема (16) приводит к системе Nв уравнений относительно Тгк не с пяти-, а с трехдиагональной матрицей, которые решаются методом матричной прогонки.
Применение неявной схемы повышает быстродействие алгоритма за счет возможности использования более крупного шага по времени, чем в случае чисто явной схемы, а применение явной схемы уменьшает несферические возмущения численной природы.
После того как расчет вспомогательного поля температуры по схемам (15)—(17) закончен, вычисляются соответствующие значения температуры межфазной границы по (18),
Далее определяются значения интенсивности фазовых преобразований ] по формуле (4), Полученное вспомогательное поле температуры Т^к используется для корректировки поля полной энергии, рассчитанной на первом этапе. Второй этап завершается вычислением давления РП+1 и температуры ТП+1 то найденным в момент времени т = тп+1 значениям плотности, вектора количества движения и полной энергии с применением уравнений состояния жидкости и газа,
3. Иллюстрация применения метода
Для демонстрации работоспособности предлагаемого подхода рассматриваются задачи о сильном сжатии кавитационного пузырька в дейтерированном ацетоне в сферически симметричной и осесимметричной постановках. Во втором случае пузырек в начале сжатия является слабо несферическим, В процессе сжатия несферичность межфазной поверхности, а в результате и несферичность движения жидкости и пара значительно возрастают,
В этих задачах существенным образом проявляются сжимаемость жидкости, теплопроводность газа и жидкости, несовершенство пара, конденсация на поверхности пузырька, наличие ударной волны, поэтому для их решения необходимо использовать большинство из возможностей рассматриваемого подхода. Результаты верификации и оценки эффективности алгоритма расчета изложены в [51, 52],
3.1. Сильное сжатие сферического пузырька
До начального момента времени кавитационный (паровой) сферический пузырек расширяется в неограниченном объеме жидкости (жидкого дейтерированного ацетона). Сначала расширение происходит под действием давления пара в пузырьке рь, превышающего давление окружающей жидкости = р0, а затем, при рь < под действием сил инерции. При £ = 0 пузырек достигает максимальных размеров, его радиус при
этом становится равным Д°, и расширение прекращается, В этот момент температура
Т0
р^ при температуре ТНа бесконечном удалении от пузырька давление и температура жидкости в процессе сжатия не изменяются, оставаясь равными (рте > р^) и Т0 соответственно. Под действием перепада давления в жид кости пузырек при £ > 0 сжимается. В конце сжатия в полости пузырька формируется ударная волна.
В процессе сжатия можно выделить две стадии: начало и финал. В начальной стадии пузырек близок к гомобаричеекому, основные градиенты давления расположены в жидкости в окрестности межфазной поверхности. В финале сжатия распределение параметров и в паре, и в жидкости является сильно неоднородным, в паре формируется ударная волна.
Расчет распределения давления и скорости при £ = 0 осуществляется в предположении несжимаемости жидкости, так как ее сжимаемость в начале несущественна. Кроме того, полагается ] = 0. Распределение плотности жидкости определяется из уравнения состояния. При описании состояния жидкого и парообразного дейтерированного ацетона применяются реалистичные уравнения в форме Ми — Грюнайзена [11].
Решение ищется в расчетной области, ограниченной искусственной границей т = та(0, ¿). В начальный момент времени она имеет форму сферы радиуса К°а. Ее перемещение в дальнейшем определяется из условия р(та,£) = р°а-, гДе Р°а соответствует давлению при £ = 0 и т = К°а. Путем численного эксперимента установлено, что для адекватного описания процессов, происходящих внутри пузырька па рассматриваемом отрезке времени, в данной задаче достаточно принять К°а = 9В°.
В начале сжатия и в газе, и в жидкости используется сетка со сгущением к меж-фазпой границе, В дальнейшем, когда в газе появляются значительные неоднородности ноля давления и начинает формироваться ударная волна, осуществляется постепенный переход па равномерную сетку по формуле (12), При сравнении результатов расчетов на разных сетках удобнее в качестве аргумента функций в и 7 принять те само время ¿, а его однозначную функцию Я(£). Иначе в силу быстрого изменения радиуса пузырька в финальной стадии сжатия степень сгущения сеток с разным числом ячеек при фиксированных значениях радиуса пузырька может сильно различаться. В данной задаче параметры 7! и 72 функции 7(Л) выбираются равными. Для сопоставления результатов расчетов, полученных на разных сетках, используется относительное время £ — ¿*, где ¿* — момент времени, когда в текущем расчете температура газа в четвертой (от центра
т
Задача решается при следующих входных данных: В° = 500 мкм, Т0 = 273.15 К, р°8 = 0.089 бар, р0 = 15.2 бар, аас = 1, ЯА = 35.05 мкм, ЯБ = 56.31 мкм, 71 = 72 = 0.01, в! = 1, в2 = 0.08.
На рис. 1-3 приводятся результаты расчетов дня пяти последовательных моментов времени: ¿1-3 (см. рис, 1) — го начальной стадии сжатия, ¿4, ¿5 (см, рис, 2, 3) — из его финальной стадии, В начале сжатия около межфазпой поверхности с обеих ее сторон возникают топкие тепловые слои (из-за наличия тепломассообмена), появляются значительные неоднородности давления в жидкости. Ко времени ¿4 передний фронт волны сжатия, сформировавшейся в результате высокой скорости движения стенки пузырька, превратился в ударную волну. В последующем интенсивность этой ударной волны быстро возрастает, а профиль давления в следующей за пей волне сжатия усложняется.
На рис. 3 представлены результаты, полученные методом настоящей работы па сетке Мд х N1 = 600 +1540 (Мд — число ячеек в области пара, а N1 — число ячеек в области жидкости) и классическим методом Годунова |31|. Расчеты методом Годунова выполнены как на сетке 600 + 1540, так и на сетках в 2 и 4 раза грубее ее. Видно, что в гладкой
Рис. 1. Радиальное распределение давления и температуры в полости еферичеекхнх) пузырька и окружающей жидкости для трех последовательных моментов времени ¿1_з начальной стадии сжатия. Точками указаны положения межфазной границы
106
p, бар
10
10
0
10 г, мкм 20
10
т, К
Рис. 2. Радиальное распределение давления и температуры внутри пузырька в два последовательных момента времени ¿4 и ¿5 финальной стадии сжатия
р, бар
10
103
10' 10
2_
10
т, К
14 ш
---20
¡1 ¡1 ¡1 ------ 30
'1 п ¡1 ¡1 зи
3000
1000 т
г, мкм
20
10
г, мкм
20
7
10 1
10
3 10
10
р, бар
Т, К
г 5 /
{
//
/;
105
1041
шЧ
г, мкм
г, мкм
5
г
г
4
Рис. 3. Радиальные распределения давления и температуры газа в полости пузырька в моменты времени ¿4, ¿5, полученные по методу настоящей работы (сплошные линии) и классическим методом Годунова (пунктир и штрихпунктир). Буквы ( ! и и означают расчеты методами Годунова и настоящей работы соответственно. Цифрами указаны расчетные сетки: 1 — 150 + 385, 2 _ 300 + 770, 3 — 600 + 1540
области различие между этими приближениями меньше, чем вблизи фронта ударной волны. При этом и в гладкой области, и в области разрыва решение, полученное методом Годунова, при измельчении сетки стремится к решению предлагаемым авторами методом. В 1511 показано, что метод настоящей работы значительно экономичнее (более чем в 10 раз) применяемых ранее |11, 12, 14, 151.
3.2. Сильное сжатие осесимметричного пузырька
Постановка данной задачи отличается от предыдущей лишь тем, что здесь в начальный момент времени поверхность пузырька не является чисто сферической, а имеет форму близкого к сфере эллипсоида. Уравнение этого эллипсоида имеет вид г8 = Я°+а2Р2(со8 0), оде а° — амплитуда отклонения от сферы в виде Р2(со8 0), Р2 _ полином
Лежапдра второй степени. Эволюция формы пузырька в процессе сжатия определяется движением жидкости и пара.
Начальное распределение параметров и положение искусственной границы г = га(в,£) при £ > 0 находятся так же, как в предыдущей задаче. При этом распределения давления, скорости и плотности при £ = 0 вычисляются с учетом того, что несферичность пузырька является малой. Давление на внешней границе р°а задается без учета песфери чпости,
В начале сжатия и в газе, и в жидкости используется, как и в предыдущей задаче, подвижная сетка, сгущающаяся по геометрической прогрессии к межфазпой поверхности, Вместе с тем здесь для повышения экономичности вычислений показатель 7(£) в начальный момент времени больше. В ходе сжатия дня корректного описания пограничных тепловых слоев оп уменьшается согласно формулы (12), По мере формирования ударной волны осуществляется постепенный переход па равномерную сетку с помощью параметра в(£) по формуле (12),
Задача решается при следующих входных данных: Я° = 500 мкм, Я°а = 9Я°, Т0 = 273.15 К,р% = 0.089 бар,р° = 15.2 бар, аас = 1, ЬА = 11.52 мке, £в = 11.55 мке, £с = 10.82 мке, = 10.99 мке, 71 = 0.08, 72 = 0.04, в1 = 0.08, в2 = 1.
На рис, 4 представлены радиальные распределения давления и температуры внутри пузырька и в относительно небольшом слое окружающей жидкости в сечениях в = 0 и п/2 для двух моментов времени ¿1,2 (¿1 < ¿2). Видно, что в момент относящийся к начальной стадии сжатия, пузырек еще близок к гомобарическому. Основные градиенты давления находятся в жидкости, а у поверхности пузырька сформировались топкие тепловые слои, К моменту ¿2, относящемуся к финальной стадии сжатия, давление пара у поверхности пузырька становится сравнимым с давлением жидкости. Внутри пузырька у его поверхности формируется ударная волна. Ее наличие лучше видно по кривым температуры, так как максимум у температуры достигается непосредственно за фронтом ударной волны, в то время как у давления — па некотором удалении от пего. В начальной стадии сжатия форма пузырька еще мало отличается от сферической, поэтому профили давления и температуры при £ = £1 на луч ах в = 0 и п/2 различаются незначительно. По мере сжатия несферичность возрастает. В результате в момент ¿2 профили давления и температуры для в = 0 и п/2 уже заметно различны. В дальнейшем различие увеличивается.
Рис. 4. Радиальные распределения давления и температуры внутри полости пузырька и окружающей жидкости в двух сечениях 9 = 0и 9 = п/2 для моментов вре мени ¿1,2 (¿1 относится к начальной, ¿2 — к финальной стадии сжатия). Точками указаны положения межфазных границ
Развитие процесса сжатия пузырька в его финальной стадии в рассматриваемой задаче иллюстрирует рис, 5 выше горизонтальных прямых з1-3, По кривым температуры видно, что по мере распространения ударной волны к центру пузырька ее интенсив-
Рис. 5. Изолинии и радиальные зависимости давления и температуры пара в пузырьке для 9 = 0 п/2 в три последовательных момента времени ¿з_5 финальной стадии сжатия. Изолинии делят диапазон изменения параметров, указанный на соответствующих осях ординат, равномерно на 15 уровней (во избежание загромождения их числовые значения не приводятся). Жирная линия — поверхность пузырька. Выше горизонтальных прямых £1_з — результаты расчетов с учетом теплопроводности, испарения и конденсации, ниже без их учета
пость довольно быстро возрастает (наличие ударной волны в профилях давления будет наблюдаться ниже па рис, 6, где дня давления используется не линейная, а логарифмическая шкала). При этом в момент ¿5 максимальные температуры на лучах в = 0 и п/2 различаются уже более чем в 2 раза, а максимальные давления — в ~ 10 раз,
В процессе сжатия форма пузырька изменяется от слегка приплюснутой по оси симметрии сферы в начале сжатия до несколько вытянутой в его финале. Фронт ударной волны по мере приближения к центру пузырька деформируется значительно сильнее. Сначала его форма подобна форме поверхности пузырька и имеет вид вертикально вытянутого эллипсоида (момент ¿3). Затем сверху и снизу появляются вмятины (момент ¿4). В ходе дальнейшего сжатия изменение формы фронта ударной волны продолжается и в момент ¿5 она становится еще сложнее. При этом поверхность пузырька практически не изменяется, оставаясь вытянутым в направлении оси симметрии эллипсоидом.
Для оценки влияния теплопроводности, фазовых переходов па поверхности пузырька в рассматриваемой задаче па рис. 5 представлены также результаты ее численного решения без учета этих эффектов, т.е. в постановке работы |52|, Соответствие моментов времени устанавливается по совпадению радиусов Я^(£) сферической составляющей в представлении фронта ударной волны в виде суммы сферических поверхностных гармоник: Я^(£3) ~ 8.9 мкм, Я^(£4) ~ 4 мкм, Я^(£5) ~ 1.4 мкм, полученных с учетом и без учета указанных эффектов. Видно, что линейные размеры пузырьков в конце сжатия при учете и без учета теплопроводности и межфазпого массообмепа значительно (примерно в 1,5 раза) различаются. Несмотря па это, изменение формы ударной волны, изолиний давления и температуры в пузырьке в этих случаях оказывается весьма близким. При этом имеют место значительные количественные расхождения. Так, без учета теплопроводности и фазовых переходов по сравнению со случаем их учета давление па луче в = 0 в момент ¿5 в ~ 2 раза меньше.
Более детальная картина вблизи центра пузырька (г < 10 мкм) в момент ¿5 показана па рис. 6. Здесь дня радиальных распределений давления и температуры используется не линейная, а логарифмическая шкала, поскольку по рис. 5 можно ошибочно заключить, что разрыва в графике давления пет. На самом деле это не так. В логариф-
Рис. 6. Картина вблизи центра пузырька (см. рис. 5) в момент ¿5 для малой окрестности центра пузырька г < 10 мкм
мичеекой шкале (ем, рис, 6) хорошо видно, что разрыв, хотя и небольшой, но имеется и находится точно в том месте, что и на графиках температуры. Изолинии давления и температуры на рис, 6 показывают, что ко времени ¿5 в профиле фронта ударной волны появляется довольно острый угол, связанный с углублением в окрестности 0 = 0 вмятины, возникшей ранее в интервале £3 < £ < ¿4. При этом температура принимает большие значения не только непосредственно за фронтом ударной волны, но и в узкой области, уходящей в виде ленты вверх от угловой точки профиля волны. Область высоких давлений не совпадает с областью высоких температур. Она располагается в окрестности оси симметрии на довольно большом удалении от фронта волны.
Представленные на рисунках результаты получены на сетке (Л5 + N1) х Л = (300 + 700) х 50 (Лд и N1 — число ячеек по радиальной координате в паре и жидкости соответственно, Л — число ячеек по угловой координате как в паре, так и в жидкости). Расчеты па сходимость показывают, что вплоть до времени ¿5 такая сетка является приемлемой. Отметим также, что затраты времени счета для достижения аналогичной точности классическим методом Годунова первого порядка были бы многократно (примерно в 103 раз) больше.
Заключение
В работе предложен метод численного решения задач сильного сжатия кавитационного пузырька, В отличие от существующих подходов он позволяет изучать динамику пара в полости пузырька и окружающей жидкости в том случае, когда пузырек является несферическим (т.е. когда задачу нужно рассматривать как неодномерную), а процесс сжатия является высокоскоростным (т, е, когда в полости пузырька могут возникать радиально сходящиеся волны сжатия, в том числе и ударные, а сжимаемость жидкости становится существенной). Описанный метод, в частности, можно применять для изучения влияния деформаций поверхности пузырька на изменение максимальных значений давления, плотности и температуры, достигаемых при сжатии в полости пузырька, по сравнению со случаем чисто сферического сжатия,
В предлагаемом подходе движение пара в пузырьке и окружающей жидкости описывается уравнениями динамики сжимаемой жидкости с учетом эффектов поверхностного натяжения, теплопроводности пара и жидкости, тепломассообмена через межфазную границу совместно с реалистичными уравнениями состояния в виде зависимостей давления и внутренней энергии от плотности и температуры. Используются подвижные сетки с выделением поверхности пузырька и сгущением ячеек для улучшения разрешения тонких тепловых пограничных слоев у поверхности пузырька и радиально сходящихся ударных волн в полости пузырька.
Каждый шаг по времени выполняется в два этапа. На первом этапе влияние теплопроводности не учитывается. Решение находится с применением ГХО-.модификации метода Годунова второго порядка точности по пространству и времени. Полученное поле полной энергии корректируется на втором этапе путем учета тепловых потоков. При вычислении тепловых потоков используется решение уравнения теплопроводности, которое находится по явно-неявной схеме первого порядка точности по времени и второго порядка точности по пространству. Применение неявной схемы позволяет использовать единый шаг по времени для обоих этапов. Эффекты испарения и конденсации рассчитываются по формуле Герца — Кнудсена — Ленгмюра,
Работоспособность метода при описании происходящих в пузырьке в условиях его сильного сжатия сложных гидродинамических процессов продемонстрирована на задаче сильного сжатия несферического (осесимметричного) кавитационного пузырька в дейтерированном ацетоне. Сначала пузырек является слабо несферическим, В ходе сжатия несферичность межфазной поверхности, а в результате и несферичность движения жидкости и пара значительно возрастают, В этой задаче существенным образом проявляются сжимаемость жидкости, теплопроводность пара и жидкости, несовершенство пара, конденсация на поверхности пузырька. Кроме того, картину движения газа в пузырьке сильно усложняет наличие ударной волны.
Предложенный метод позволяет получить решение поставленной задачи за приемлемое время с приемлемой точностью с учетом всех перечисленных выше особенностей.
Список литературы
fl] Gaitan D.F., Crum L.A. Observation of sonoluminescence from a single, stable cavitation bubble in a water/glvcerine mixture // 12th Intern. Svmp. On Nonl. Acoustics. New York: Elsevier, 1990. P. 459-463.
[2] Brenner M.P. Single-bubble sonoluminescence // Rev. Modern Physics. 2002. Vol. 74. P. 425-484.
[3] Маргулис M.A. Сонолюминесценция // Успехи физ. наук. Обзоры актуальных проблем. 2000. Т. 170, № 3. С. 263-287.
[4] Taleyarkhan R.P., West C.D., Сно J.S. ет al. Evidence for nuclear emissions during acoustic cavitation // Science. 2002. Vol. 295. P. 1868-1873.
[51 Taleyarkhan R.P., West C.D., Lahey R.T. et al. Nuclear emissions during self-nucleated acoustic cavitation // Phvs. Revol. Lett. 2006. Vol. 96. 034301
[6] Wu C.C., Roberts P.H. Shock wave propagation in a sonoluminescencing gas bubble // Ibid. 1993. Vol. 70. P. 3424-3427.
[71 Moss W.C., Clarke D.B., White J.W., Young D.A. Hydrodvnamic simulations of bubble collapse and picosecond somoluminescence // Phvs. Fluids. 1994. Vol. 6, No. 9. P. 2979-2985.
[8] Barber B.P., Hiller R.A., Lopstedt R. et al. Defining the unknowns of sonoluminescence // Phvs. Rep. 1997. Vol. 281. P. 65-143.
[91 Moss W.C., Clarke D.B., Young D.A. Calculated pulse widths and spectra of a single sonoluminescencing bubble // Science. 1997. Vol. 276. P. 1398-1401.
[101 Moss W.C., Young D.A., Harte J.A. et al. Computed optical emissions from a sonoluminescing bubble // Phvs. Revol. E. 1999. Vol. 59, No. 3. P. 2986-2992.
[Ill Nigmatulin R.I., Akhatov I.Sh., Topolnikov A.S. et al. The theory of supercompression of vapor bubbles and nano-scale thermonuclear fusion // Phvs. Fluid. 2005. Vol. 17. P. 107105.
[121 Kondic L., Gersten J.I., Yuan C. Theoretical studies of sonoluminescence radiation: Radiative transfer and parametric dependence // Phvs. Revol. 1995. Vol. 52. P. 4976-4990.
[131 Storey B.D., Szeri A.J. Water vapour, sonoluminescence and sonochemistrv// Proc. R. Soc. bond. A. 2000. Vol. 456. P. 1685-1709.
[141 Akhatov I., Lindau O., Topolnikov A. et al. Collapse and rebound of a laser-induced cavitation bubble // Phvs. Fluids. 2001. Vol. 13, No. 10. P. 2805-2819.
[15] Аганин А.А., Ильгамов М.А. Динамика пузырька газа в центре сферического объема жидкости // Мат. моделирование. 2001. Т. 13, № 1. С. 26-40.
[16] Kull H.J. Excitation and amplification of cavity oscillations in spherically converging shells / / Phvs. Re vol. A. 1990. Vol. 41, No. 8. P. 4312-4325.
[17] Илькаев P.И., Гаранин С.Г. Исследование проблем термоядерного синтеза на мощных лазерных установках // Вестник РАН. 2006. Т. 76, № 6. С. 503-513.
[18] Лево И.Г., Тишкин В.Ф. Исследование гидродинамической неустойчивости в задачах лазерного термоядерного синтеза методами математического моделирования. М.: Физ-матлит, 2006. 304 с.
[19] Prosperetti A. Viscous effects on perturbed spherical flows // Quarterly Appl. Math. 1977. Vol. 34. P. 339-352.
[20] Keller J.В., Miksis M. Bubble oscillations of large amplitude //J. Acoust. Soc. Amer. 1980. Vol. 68, No. 2. P. 628-633.
[21] Brenner M.P., Lohse D., Dupont T.F. Bubble shape oscillations and the onset of sonoluminescense // Phvs. Revol. Lett. 1995. Vol. 75, No. 5. P. 954-957.
[22] Prosperetti A., Hao Y. Modeling of spherical gas bubble oscillations and sonolumines-cence // Phil. Trans. R. Soc. bond. A. 1999. Vol. 357. P. 203-223.
[23] Putterman S.J., Weninger K.P. Sonoluminescence: How bubbles turn sound into light // Annu. Revol. Fluid Mech. 2000. Vol. 32. P. 445-476.
[24] Lin H., Storey B.D., Szeri A.J. Inertiallv driven inhomogeneities in violently collapsing bubbles: The validity of the Ravleigh-Plesset equation //J. Fluid Mech. 2002. Vol. 452. P. 145-162.
[25] Lin H., Storey B.D., Szeri A.J. Ravleigh-Tavlor instability of violently collapsing bubbles // Phvs. Fluids. 2002. Vol. 14. P. 2925-2928.'
[26] Аганин А.А., Илвгамов M.A., Гусева T.C. Искажение сферической формы пузырька при больших расширениях-сжатиях из состояния покоя // Динамика газовых пузырьков и аэрозолей. Казань: Изд-во КГУ, 2003. С. 95-132.
[27] Аганин А.А., Илвгамов М.А., Топорков Д.Ю. Затухание начального искажения сферической формы пузырька // Там же. С. 66-94.
[28] Топорков Д.Ю. Динамика газового пузырька при периодическом изменении давления окружающей жидкости // Там же. С. 179-215.
[29] Аганин А.А., Гусева Т.С. Изменение малого искажения сферической формы газового пузырька при его сильном расширении-сжатии // ПМТФ. 2005. № 4. С. 17-22.
[30] Аганин А.А., Илвгамов М.А., Косолапова Л.А. и др. Эллипсоидальные колебания газового пузырька при периодическом изменении давления окружающей жидкости // Механика жидкости и газа. 2005. № 5. ('.15 52.
[31] Нигматулин Р.И., Аганин А.А., Илвгамов М.А., Топорков Д.Ю. Искажение сферичности парового пузырька в дейтерированном ацетоне // Докл. РАН. 2006. Т. 408, № 6. С. 767-771.
[32] Eller A.I., Crum L.A. Instability of the motion of a pulsating bubble in a sound field //J. Acoust. Soc. Amer. Suppl. 1970. Vol. 47, No. 3. P. 762-767.
[33] hllgenfeldt S., Lohse D., Brenner M. Phase diagrams for sonoluminescing bubbles // Phvs. Fluids. 1996. Vol. 8, No. 11. P. 2808-2826.
[34] Hilgenfeldt S., Brenner XL. Grossmann S., Lohse D. Analysis of Ravleigh-Plesset dynamics for sonoluminescing bubbles //J. Fluid Mech. 1998. Vol. 365. P. 171-204.
[35] Нао Y., Prosperetti A. The effect of viscosity on the spherical stability of oscillating gas bubbles // Phvs. Fluids. 1999. Vol. 11, No. 6. P. 1309-1317.
[36] Augsdorper U.H., Evans A.K., Oxley D.P. Thermal noise and the stability of single sonoluminescing bubbles // Phvs. Rev. E. 2000. Vol. 61. P. 5278-5286.
[37] Yuan L., Ho C.Y., Сни МЛ'.. Leung P.T. Role of gas density in the stability of single-bubble sonoluminescence // Ibid. 2001. Vol. 64. P. 016317.
[38] Kwak H.-Y., Karng S.W., Lee Y.P. Ravleigh-Taylor instability on a sonoluminescencing gas bubbles// J. Korean Phvs. Soc. 2005. Vol. 45, No. 4. P. 951-962.
[39] Ильгамов M.A. Качественный анализ развития отклонений от сферической формы при схлопывании полости в жидкости // Докл. РАН. 2005. Т. 401, № 1. С. 37-40.
[40] Calvisi M.L., Lindau О., Blake J.R. ет al. Shape stability and violent collapse of micro-bubbles in acoustic traveling waves // Phvs. Fluids. 2007. Vol. 19. P. 047101-047101-15.
[41] Plesset M.S., Chapman R.B. Collapse of initially spherical vapour cavity in the neighbourhood of a solid boundary // J. Fluid. Mech. 1971. Vol. 47, No. 2. P. 283-290.
[42] Blake J.R., Hooton M.C., Robinson P.В., Tong R.P. Collapsing cavities, toroidal bubbles and jet impact // Phil. Trans. R. Soc. Lond. A. 1997. Vol. 355. P. 537-550.
[43] Apanasiev K.E., Grigorieva I.V. Numerical investigation of three-dimensional bubble dynamics // J. Eng. Math. 2006. Vol. 55, No. 1-4. P. 65-80.
[44] Bremond N., Arora XL. Dammer S.M., Lohse D. Interaction of cavitation bubbles on a wall // Phvs. Fluids. 2006. Vol. 18. P. 121505-121505-10.
[45] Blake J.R., Gibson D.C. Growth and collapse of a vapour cavity near a free surface //J. Fluid Mech. 1981. Vol. 111. P. 123-140.
[46] Blake J.R., Keen G.S., Tong R.P., Wilson M. Acoustic cavitation: The fluid dynamics of non-spherical bubbles // Phvs. Trans. R. Soc. Lond. A. 1999. Vol. 357, No. 6. P. 251-267.
[47] Colella P., Puckett E.G. Modern Numerical Methods for Fluid Flow / Eds. U.C. Berkeley and U.C. Davis. Univ. California, 1994.
[48] Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. XI.: Физматлит, 2001. 608 с.
[49] Бондаренко Ю.А., Башуров В.В., Янилкин Ю.В. Математические модели и численные методы для решения задач нестационарной газовой динамики. Обзор зарубежной литературы. Саров: РФЯЦ-ВНИИЭФ, 2003. 53 с.
[50] Leer В. Review article: Upwind and high-resolution methods for compressible flow: From donor cell to residual-distribution schemes // Commun. Comput. Phvs. 2006. Vol. 1, No. 2. P. 192-206.
[51] Аганин А.А., Халитова Т.Ф., Хисматуллииа H.A. Расчет сильного сжатия сферического парогазового пузырька в жидкости // Вычисл. технологии. 2008. Т. 13, № 6. С. 17-27.
[52] Аганин А.А., Илвгамов М.А., Халитова Т.Ф. Моделирование сильного сжатия газовой полости в жидкости // Мат. моделирование. 2008. Т. 20, № 11. С. 89-103.
[53] Harten A., Engquist В., Osher S., Chakravarthy S.R. Uniformly high order accurate essentially non-oscillatorv schemes III //J. Сотр. Phvs. 1987. Vol. 71. P. 231-303.
[54] Илвгамов M.A. Колебания оболочек, содержащих жидкость и газ. XI.: Наука, 1969. 182 с.
Поступила в редакцию 30 апреля 2009 г.
