CC BY
143
Вычислительные технологии
Том 8, № 5, 2003
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ УДАРНЫХ ВОЛН В ПУЗЫРЬКОВЫХ СИСТЕМАХ
В. А. Вшивков, Г. Г. Лазарева Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: vsh@ict.nsc.ru
A numerical model of pressure and velocities field dynamics in the process of interaction of a plain shock wave with a “free” bubble system is suggested in the axisymmetric case.
A comparison of the explicit and the splitting finite-difference schemes is made. Boundary conditions and the method of definition of error calculation of the exact value of maximum pressure amplitude in the focusing region are considered in details on the basis of estimation of the energy variation over the emanation wave.
Введение
Среди существующих теоретических моделей сазера (акустического аналога импульсных лазерных систем) известна схема, аналогичная схеме лазера на свободных электронах. Она предполагает в качестве активной среды жидкость и в качестве дисперсных частиц газовые пузырьки [1]. Так как для экспериментальных исследований в этом направлении наиболее принципиальными могут оказаться “свободные” системы, особенно интересны исследования поведения ударных волн в системе распределенных в жидкости пузырьковых кластеров и механизмов их усиления. Существует много методов решения пространственных задач нестационарной газовой динамики и уравнений химической кинетики. Однако при исследовании волновых процессов в бинарных средах возникает трудность разработки алгоритмов для описания быстропротекающих процессов в пузырьковых средах с физикохимическими превращениями. Это обусловлено не только широким диапазоном временных и пространственных масштабов течений в такого рода системах, но и особенностями, связанными с большими амплитудами ударных волн, генерируемых в кластере (порядка десятков мегапаскалей). Численное моделирование акустически активных систем, способных к генерации мощного излучения, исследование динамики их состояния и волновых процессов позволяют находить значения параметров, которые трудно, а порой и невозможно измерить в ходе физического эксперимента. В связи с этим были разработаны адекватные изучаемым явлениям физико-математические модели, алгоритмы и методы их численной реализации, позволяющие исследовать взаимодействие волн с динамически меняющимися акустическими системами [2]. Вопрос о взаимодействии ударной волны с пузырьковым кластером изучался в рамках квазиодномерной постановки задачи в [3], где получен эффект цилиндрической фокусировки на ось в пассивной пузырьковой среде. Динамика волновых процессов в химически активных пузырьковых средах в лагранжевой
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2003.
системе координат для цилиндрической симметрии рассмотрена в [4-8]. Ниже эта задача рассматривается на примере решения двумерной аксиально-симметричной задачи о столкновении плоских ударных волн с пузырьковыми кластерами сферической и тороидальной формы.
1. Постановка задачи. Основные уравнения
На торце кругового цилиндра радиуса гтах, заполненного некоторой жидкостью, в момент
і = 0 задается скачок давления. На оси г цилиндра на расстоянии Ісі от этой границы расположен центр пузырькового кластера с объемной концентрацией газовой фазы ко. Пузырьки газа в кластере имеют один и тот же радиус До. При і > 0 ударная волна, распространяясь вдоль положительной оси г, сталкивается с облаком пузырьков, огибает его и в зоне контакта фронта преломляется в кластер. Состояние реальных жидкостей с микронеоднородностями носит явно двухфазный характер, что позволило использовать в качестве управляющей системы модифицированную систему Иорданского — Когарко — ван Вингардена (ИКВ-модель) [6]. Рассматриваемая модель базируется на известных законах сохранения массы и импульса для средних значений давления р, плотности р, скорости и и подсистеме, определяющей состояние среды.
Система двумерных уравнений газовой динамики:
др + аіу(ри) = 0, ди + и(Уи) = -1 Ур, (1)
ді ді р
где
Р = р(р) = 1 +
про
1к
к = ГНГрв3-
1 - ко
Система уравнений (1) записана в безразмерных переменных, где в качестве нормирующих величин взяты: ь0 = 10 м/с, р0 = 1 г/см3, Ь = 1 см, £0 = L/v0 = 10-3 с, р0 = рь% = 106 дин/см2 ^ 1 атм. Здесь р0 — плотность жидкости, с0 — скорость звука, р0 — давление в жидкости, к — удельная доля газовой фазы в кластере, в = К/К0 — относительный радиус пузырьков. В дальнейшем данная система уравнений будет рассматриваться в цилиндрической системе координат (г, г).
Заметим, что особенность модели состоит в том, что она позволяет игнорировать наличие пузырьков и считать среду однородной с особыми свойствами состояния, описываемого уравнением для к0. Движением пузырькового кластера можно пренебречь, так как сдвиг кластера за все время рассматриваемого процесса много меньше пространственного шага сетки.
Система уравнений состояния пузырьковой среды:
дв = „
31 = ^
^ = -^ "2 - С. - с" - р - в-3^, (2)
дС = __! с 2_ С-
ді 2в в2 2 в2
где
^ _ 2^ ^ _ 2^ ^
Сі = --, С2 = —--------, Сз
Щро’ Коро' Во\/роро ’
2
роСо
а — коэффициент поверхностного натяжения, ^ — коэффициент вязкости. Расчеты проводились в диапазонах изменения параметров п =7.15, к0 = 0.001... 0.1, Я0 = 0.01... 0.4 см. Пределы изменения параметров аналогичны рассмотренным в работе [8].
Прямым следствием системы уравнений газовой динамики (1) является уравнение
Здесь первое слагаемое — плотность энергии, состоящая, в свою очередь, из кинетической и потенциальной частей, второе и третье — потоки энергии через границы. Отсюда видно, что при нулевых компонентах вектора скорости на границах энергия внутри области сохраняется. Таким образом, уравнение (3) позволяет согласовывать граничные условия и проводить контроль правильности решения.
Область решения в цилиндрических координатах представляет собой прямоугольник 0 < г < гтах, 0 < г < гтах. Граничные условия (г = 0) описывают постоянную ударную волну с амплитудой р = Р3^ путем задания осевой компоненты скорости в предположении равенства нулю радиальной компоненты. Расчеты проводились в диапазонах изменения параметра Р3^ =10 ... 120 атм. На границе г = 0 задаются граничные условия симметрии. На границе г = гтах необходимо поставить такое граничное условие, чтобы возмущения уходили из области, не возвращаясь. Выход волны из области при г = гтах определен равенством нулю вторых производных всех функций по осевому направлению. Подробнее эти граничные условия будут описаны в разд. 4.
2. Разностные схемы
Рассмотрим два разностных метода, используемых для решения газодинамической системы уравнений (1): явную схему с направленными разностями и схему расщепления, описанную в [9] и адаптированную к данной задаче. В качестве схемы для реализации первого этапа может быть взята и другая схема [10].
2.1. Явная схема с направленными разностями
Введем новые переменные
(3)
пги . . _. Л
------ (р + р0(с0 - 1)м = °.
п - 1
Я = гр, Р (Я) = (р(р)), и = пг, V = п.
В новых переменных система уравнений газовой динамики (1) примет вид
зи гг3и т,3и 1 3Р
—+ и--------------+ V— =--------,
З£ Зг Зг Я Зг
зу изу &у = _^ЗР
З£ + Зг + Зг Я Зг.
В двумерной области решения введена равномерная прямоугольная сетка с узлами
На введенной сетке газодинамическая часть (4) исходной системы уравнений аппроксимируется явной схемой с направленными разностями [12]. При этом система уравнений разбивается на две подсистемы, связь между которыми осуществляется только через вычисление давления. Это облегчает контроль над работой схемы. Таким образом, схема решения задачи на одном шаге разбивается на два этапа. На первом (газодинамическом) этапе вычисляются давление и плотность, затем определяются осевая и радиальная компоненты скорости. На втором (кинетическом) этапе решается система обыкновенных дифференциальных уравнений. Выбор явной схемы с направленными разностями обусловлен ее надежностью, простотой и удобством контроля численных результатов.
Для реализации этого численного метода необходимо перейти от функций с непрерывными аргументами к дискретным наборам чисел, их заменяющих. Для этого введем сеточные функции
Т'г — (І 1/2)^г, І — !)•••) ^тах + 1)
0,
где
/ — любая из функций и, V, Д.
2.2. Схема расщепления
Систему уравнений газовой динамики (1) с использованием новых переменных
Д — гр, Р (Д) — (р(р)), и — григ, V — гри,
запишем в дивергентном виде:
д д д
-Д + —и + — V — 0,
ді дг дг
д д д
діи + дГ (Чи + Рг) + дг(Чи*) — Р (5)
д д д
+ дг (^) + ^ (Vй, + Рг) —
На введенной выше сетке рассмотрим аппроксимацию газодинамической части системы уравнений (5) схемой расщепления. Схема состоит из предиктора с расщеплением по координатам г и г и корректора. Схема расщепления интересна абсолютной устойчивостью (при 0.5 < а < 1), что позволяет сократить время счета.
Описанным выше образом перейдем от функций с непрерывными аргументами к дискретным наборам чисел, их заменяющих. В качестве предиктора используется неявная схема расщепления по направлениям. На первом этапе предиктора учитываются производные только по направлению г:
и
п+1/4
і,к
та
+
- Чі . Дг(и"к+1/4ь?к/ДіПк)
+г-
Д * рп+1/4 Дг Рі,к
кг
Д“+1/4 - Д?к + Дг и”к+1/4
та
Кг
(6)
п+1/4 іт д (\Гп+1/4Т Тп / \
і,к *і,к + Дг( Vі,к иі,к/Ді,к)
та К
где Дг определена указанным выше образом,
Д*/іПк
/п /п
*/ і,к */ і— 1
/і
к,
п /п
і+1,к о і,к,
и™к < о, ип > 0,
f — любая из функций [/, V, Д.
Получим неизвестное значение рп+1/4 с помощью разложения в ряд Тейлора, используя значение функции давления на предыдущем временном слое п:
Р
п+1/4
дР
і,к
др
і,к
др
ді
РП +
2 п+1/4
таа2 гіРі,к
гі Р'ік
і,к
та
0
0
0
П
г
где а2 = др/дР на временном слое п.
Выразив др/дР из уравнений (6), получим
рП+1/4 ____ рп
рг,к р г,к
2 Л ип+1/4 таа2 Лгц,-к
к,
На втором этапе предиктора учитываются производные только по направлению г:
Кк+1/2 - Кк+1/“ . л.(К'к+1/2у-пк/дпк)
та
к,
Л * рп+1/2 . Л. р г,к п
+Г--------------- = 0,
к,
т-)п+1'2 Оп+1/4 Д Т/п+1/2
Дг,к — Дг,к + Л ^г,к = о
та
к.
(7)
иП+1/2 - и"+1/4 . л,(и;‘+1/2у-пк^)
г,к
+
та
к,
где Л, определена указанным выше образом,
л;/Пк
п п п
/г,к /г,к-1, Ч,к < 0,
/-- - /П, V: > о,
</ г,к+1
/ — любая из функций Ц, V, Д.
Аналогичным образом получим неизвестное значение Р™+1/2 с помощью разложения в ряд Тейлора, используя значение функции давления на слое п и выражение для др/д£ из уравнений (6), (7):
п+1/2 = РПк + ^ (Лг иПк+1/4/к, + Л, У,П+1/2Л;).
Г г
г,к
В качестве корректора используется явная схема, восстанавливающая консервативность схемы расщепления:
д:к‘ - д^ + л, иу1/2 + л, к,к
п+1/2
т
к,
к,
- шпк+л, (и::+1/2и::+1/2/д;:+1/2)
т
к,
Л*(гг^1/2) , Л,(иТ+1/2К:к+1/2/ДП+1/2) „п+1/2
к,
+
к,
Р
г,к
V?1 - У,к + л, (С^Рц^/Да11*)
т
л,«:кн/2к;+1/2/ДП+1/2). д;рп+1/2
к,
кг
+ Г
к,
0.
Г
0
0
2.3. Решение обыкновенных дифференциальных уравнений
Для расчета подсистемы уравнений состояния пузырьковой среды применялась неявная схема Рунге — Кутты — Мерсона 4-го порядка [12]. Так как график решения может иметь асимптоты, близкие к вертикальным, шаг т может измельчаться вплоть до нуля. В уравнениях (2) перейдем от аргумента £ к длине дуги функции £(£) по формуле
Тогда вместо системы (2) получаем следующую систему уравнений:
= (1 + с2)-1/2, Щ = й (1 + С2)-1/2, д® = С(1 + с2)-1/2,
где О — правая часть второго уравнения состояния пузырьковой среды, переписанного в виде д®/д£ = О. Переход к аргументу £ позволяет использовать любой из методов решения обыкновенных дифференциальных уравнений с постоянным изменением шага, ограниченного размером длины дуги (порядка т) [13].
3. Постановка граничных условий
Реализация граничных условий для входа и выхода ударной волны из области в случае ее взаимодействия с пассивным сферическим кластером имеет несколько характерных особенностей.
3.1. Граничные условия на входе
При постановке граничных условий для входа ударной волны в расчетную область используется адиабата Гюгонио
= (Р - 1)(Р - Ро)-
Непосредственное выполнение граничного условия из адиабаты Гюгонио не дает удовлетворительных результатов в процессе решения поставленной задачи. В ходе методических расчетов выяснилось, что когда возмущение, инициированное пузырьковым кластером, возвращается ко входу в область (г = 0), возникают численные колебания значений функций. В связи с этим была произведена модификация граничных условий следующим образом: если градиенты давления и осевой компоненты скорости имеют разные знаки на входе в область, то задается условие равенства нулю производной давления по осевому направлению. Осевая компонента скорости задается, как и раньше. При такой постановке граничных условий возмущения, инициированные пузырьковой средой, возвращаются ко входу в область и проходят сквозь границу (г = 0). При этом не возникает колебаний значения функций. Когда излученные кластером ударные волны не достигают входной границы, работает граничное условие из адиабаты Гюгонио.
3.2. Граничные условия на выходе
Традиционные граничные условия на выходе плоской ударной волны из области (г = гтах) плохо работают в силу существенно нелинейной зависимости р(р). Последовательный перебор всевозможных граничных условий для выхода ударных волн показал, что наиболее
приемлемые результаты демонстрирует самое простое граничное условие
^ = 0, ЁР = 0, ^ = 0.
дг , дг , дг
3.3. Исследование граничных условий в одномерной постановке
Исследуем зависимость решения от граничных условий в одномерном случае. Система одномерных уравнений газовой динамики имеет вид
др др^ д^ д^ 1 др д£ + дг 0, д£ + ^дг рдг-
Рассмотрим прохождение ударной волны в чистой жидкости в области (0, гтах). Скорость ударной волны, согласно адиабате Гюгонио, равна 0.15 см/мкс. Следовательно, время пробега волны расчетной области длиной гтах = 3 см равно 20 мкс. Решение имеет следующий вид. Ударная волна (рис. 1, а) выходит из области через правую границу (рис. 1, б) , при этом небольшое возмущение отражается от границы выхода (рис. 1, в) и движется к левой границе расчетной области со скоростью около 0.2 см/мкс. Возникшее таким образом несущественное колебание, не влияющее на характер решения, достигнув левой границы области, начинает взаимодействовать с граничными условиями на входе (рис. 1, г). Такое взаимодействие приводит к сильному искажению решения (рис. 1, д). Большое время счета
Рис. 1. Взаимодействие граничных условий в одномерной постановке. Распределение давления (на всех графиках Р — относительное давление, з — расстояние, см) для различных моментов времени: а — £ =10 мкс, б — £ = 18 мкс, в — £ = 22 мкс, г — £ = 35 мкс, д — £ = 50 мкс, е — £ = 70 мкс.
показывает неограниченное падение амплитуды давления волны (рис. 1, е). Таким образом, реализация граничных условий для входа и выхода ударной волны из области требует их согласования. Постановка независимых граничных условий приводит к неустойчивости и, следовательно, нефизичности решения. Постановка граничного условия на входе волны в область, описанная в п. 3.1, позволяет избежать нежелательных эффектов взаимодействия условий для входа и выхода на границах области.
4. Метод нахождения значения максимальной амплитуды давления в области фокусировки
Свойства используемой явной схемы предполагают хорошее качественное соответствие результатов расчета характеру истинного физического процесса. Использование схемы именно первого порядка дает уверенность в том, что большие градиенты решения, получаемые в результате расчетов, не являются счетными эффектами и не усилены ими. Таким образом, выбор схемы обусловлен ее простотой, надежностью, удобством и хорошими результатами расчетов, полученными по этой схеме. Побочным эффектом такого выбора схемы является сглаживание острых пиков решения, например максимума амплитуды давления. Следовательно при больших шагах пространственной сетки максимальная величина давления, полученная в результате расчетов, существенно меньше ее действительного значения. Рост максимума давления Р/ с измельчением шага по пространству к и времени т демонстрирует табл. 1. В ней приведены результаты расчета при параметрах Лсг = 5 см, гсг = 10 см, к = 0.01, Д0 = 0.2 см, Р<^ = 30 амп. Таким образом, имея результат расчета, можно с уверенностью рассматривать его как некоторую оценку снизу для максимальной амплитуды давления в момент фокусировки. К сожалению, точного значения максимальной амплитуды давления в момент фокусировки схема дать не может вследствие вышеизложенных причин. Но есть возможность определить искомые величины с известной погрешностью. Для этого необходимо получить оценки пика решения как снизу, так и сверху. Ниже описан метод, построенный на основе оценки величины изменения энергии волны переизлучения, для определения погрешности нахождения точного значения максимума давления в области фокусировки.
Обратимся к одномерной аналогии описываемого метода. Предлагаемый метод основан на следующих свойствах функции давления Р(г) в окрестности максимума. Для определенности будем считать, что окрестность максимума амплитуды давления включает в себя все значения Р(г) > Р1. Так как схемы первого порядка “размазывают” острые пики решения на несколько ячеек, каждому разбиению по пространству соответствует своя функция Р(г). Расчеты показывают, что в окрестности максимума существует достаточно широкая область значений г, в которой функция Р(г) выпуклая при любых параметрах сетки. Свойством используемой разностной схемы является слабая зависимость от параметров
Таблица 1
н т Р/
0.1 0.2 505.83
0.05 0.2 584.86
0.05 0.1 603.56
0.025 0.1 609.68
0.025 0.05 715.64
сетки значения соотношения
Е (р)/У =
и >р
Р^г
Л
р >Р1
Можно показать, что на заданном промежутке наибольшим максимумом среди всех выпуклых функций Р(г), площадь подграфика которых Е(р) равна заданной величине, обладает функция, имеющая график в виде треугольника АВС (рис. 2, а). Таким образом, зная площадь подграфика, можно однозначно определить высоту треугольника, которая равна максимальной амплитуде давления. Определенная таким образом высота и является искомой оценкой сверху. Выбор линейного распределения, или треугольника в одномерном случае, обусловлен не только простотой его геометрии и большим количеством точек, лежащих в его основании. Как видно из рис. 2, б, при измельчении пространственной сетки пик решения в пятне фокусировки становится круче и острее. Следовательно, гипотеза о линейности функции давления вполне обоснована.
Распространим эти рассуждения на пространственный случай. Рассмотрим трехмерную область, где амплитуда давления принимает самые высокие значения Р(г, г, <^) > Р1
а
Рис. 2. Схема, иллюстрирующая одномерную аналогию метода определения погрешности нахождения значения максимальной амплитуды давления в области фокусировки (а) и изолинии давления в пятне фокусировки для расчетных сеток с шагами Н = 0.1 см, т = 0.2 мкс (слева) и Н = 0.025 см, т = 0.05 мкс (справа) при параметрах задачи к = 0.01, До = 0.2 см, Д4 = 2.25 см, ДС1 = 1.25 см (б).
Рис. 3. Распределение давления в момент фокусировки: во всей расчетной области (а) и в пятне фокусировки (б); £ = 180 мкс, к = 0.01, До = 0.2 см, ДС1 = 5 см.
(рис. 3). Изучение динамики максимума амплитуды давления при измельчении шагов расчетной сетки, т.е. при приближении расчетного решения исходной задачи к точному, позволяет найти требуемую оценку сверху. Рассматриваемая область решения представляет собой осесимметричный объем, в котором расчетное распределение давления является нелинейным. Предположим, что при устремлении пространственного и временного шагов к нулю давление распределено по рассматриваемой области линейно от ее границ к центру. Назовем график функции линейного распределения четырехмерным конусом. Тогда основанием конуса является область в трехмерном пространстве, аналогом которой выступает в одномерном случае основание трегольника АВС. Образующие конуса аналогичны в одномерном случае сторонам треугольника.
Таким образом, построен конус, отражающий предполагаемое в предельном случае линейное распределение давления. Высота этого конуса, а следовательно, и значение давления в точке, характеризующейся максимальной амплитудой давления, однозначно определяются через значение потенциала объемного действия сил давления Е (р) и объем основания конуса V.
Запишем их отношение в цилиндрических координатах для сеточных функций:
2п ^ Рг,к Г;Л,Г ^
Р1,к>Р 1
2п ^ Г;Л,Г ^ '
Рг,к >Р1
Произведем оценки величины изменения энергии волны переизлучения. В табл. 2 приведены результаты расчета при параметрах = 1.25 см, гсг = 2.25 см, к = 0.01, Д0 = 0.2 см,
= 30 атм. Как видно из таблицы, выражающее количество энергии в единице объема соотношение Е (р)^ остается практически постоянным вне зависимости от шага сетки по времени и пространству. Следовательно, описываемый метод позволяет достаточно точно оценить сверху максимальное значение амплитуды давления.
Целесообразным критерием определения окрестности максимума явлется достаточное для определения искомых энергий число точек. Для определения Р1 предлагается выбрать значение амплитуды давления в точке перегиба распределения давления вдоль оси 0^ от левой границы области фокусировки (со стороны падающей на кластер волны) до характеризующейся расчетным максимумом амплитуды давления точки этой области. Связь критерия выбора значения Р1 с решением позволяет получать достаточно большое количество точек, лежащих в пятне фокусировки, при любых вычислительных параметрах. Метод сохраняет физическую интерпретацию, так как определенная по этому критерию окрестность максимума обладает всеми свойствами области фокусировки. Критерий нагляден и прост в реализации, но может быть заменен на любой другой. Выбор точки перегиба именно на такой кривой обусловлен такими особенностями структуры ударных волн
Таблица 2
н т Р/ Р1 Е/2п У/2п Е/У
0.1 0.4 107.63 50.0 67.92 1.13 60.09
70.0 15.47 0.198 78.12
0.05 0.2 124.31 50.0 79.44 1.29 61.61
70.0 22.59 0.267 84.78
0.025 0.1 136.24 50.0 83.61 1.35 62.12
70.0 22.52 0.264 85.40
в пузырьковых средах, как существование предвестников. Так как используемая схема сохраняет качественную картину рассматриваемого процесса при любых допустимых разбиениях по пространству и времени, при устремлении шагов сетки к нулю точка перегиба не будет изменять свое месторасположение. С измельчением сетки значение амплитуды давления Р1 в точке перегиба будет расти. Для удобства численной реализации метода Р1 можно определять как некоторое среднее значение амплитуды давления в точке перегиба при различных шагах сетки.
Определив Р1 и вычислив Е(p)/V на сетке с наименьшими шагами сетки по пространству и времени, можно определить высоту конуса того же объема. Высота такого конуса Рт и будет требуемой оценкой сверху для значения максимальной амплитуды давления. Для определения Рт используем выражения для объема конуса в цилиндрических координатах, подставив его в интегральное выражение для Е (р):
Следовательно Рт = 4Е(р)^ — 3Р1.
Разработанный метод позволяет находить значение максимального давления в области фокусировки с точностью 2-5 % при достаточно мелкой сетке. Шаги по времени и пространству для сетки, использованной при расчете для получения оценки снизу, не должны превышать 0.1 мкс и 0.025 см соответственно. Таким образом, предлагаемый метод позволяет определить с высокой точностью величину максимального давления, знание которой имеет важное значение при проведении экспериментов.
5. Результаты численных расчетов
5.1. Постановка тестовых задач
Исследование численного алгоритма проводилось на модельной задаче о прохождении плоской стационарной ударной волны по чистой жидкости. Расчетная область этого течения представляет собой цилиндрическую область [0,гтах] х [0,гтах]. На левой входной границе задается осевая компонента скорости в предположении равенства нулю радиальной компоненты, при этом генерируется постоянная ударная волна с амплитудой .
Прохождение плоской стационарной ударной волны хорошо описывается приведенной выше явной схемой с направленными разностями. Тестовые расчеты показывают, что ширина фронта пропорциональна длине шага сетки. Таким образом, ширина фронта определяется количеством точек, и, измельчая разбиение области, можно увеличивать крутизну фронта. Расчеты тестовой задачи с использованием схемы расщепления показали, что она хорошо описывает прохождение плоской ударной волны. С уменьшением параметра а ширина фронта уменьшается, но начиная с некоторого значения а, близкого к 0.5, возникают осциляции решения (рис. 4). Измельчение шага сетки приводит к усилению крутизны фронта падающей волны. Недостатком схемы расщепления в сравнении с явной схемой является большее время счета (в четыре раза) при единых параметрах расчета. Этот недостаток частично компенсируется увеличением временного шага в два раза с сохранением качества решения по сравнению с используемым шагом при счете по явной
я
я
о
о
Рис. 4. Распределение давления на фронте ударной волны, полученное с помощью схемы расщепления, при различных значениях параметра а.
Таблица 3
н ||Р1 - Рн11 ||Р2 - Рн11
0.1000 2.27 1.50
0.0750 1.56 1.40
0.0500 1.27 0.86
0.0250 0.67 0.49
0.0100 0.31 0.20
0.0050 0.18 0.12
0.0025 0.10 0.06
схеме с тем же разбиением по пространству. Отметим, что структура решения сохраняется и при увеличении временного шага в десять раз. Ширина фронта, получаемая из расчетов с использованием этих разностных схем при единых параметрах, практически одинакова. Полная энергия системы в расчетах по явной схеме и схеме расщепления сохраняется одинаково хорошо. Из табл. 3, где приведены квадраты отклонения давления от точного решения в норме Ь2, следует, что с измельчением шага пространственной сетки значения сеточных функций стремятся к соответствующим точным значениям со скоростью О (К) для явной схемы и О (К2) для схемы расщепления. Как в случае явной схемы (р1), так и в случае схемы расщепления (р2) правомерно говорить о близости решений разностной и дифференциальной задач.
5.2. Взаимодействие ударных волн в активной пузырьковой системе
В качестве иллюстрации разработанной численной модели приведем результаты расчетов задач о взаимодействии плоской ударной волны со сферическим и тороидальным облаком пузырьков [14, 15]. Прохождение ударной волны с амплитудой РзН = 3 МПа через сферический и тороидальный пузырьковый кластеры для двумерного случая показан на рис. 5 и 6 в виде распределения относительного (в единицах гидростатического давления р0 = 0.1 МПа) давления р в пространстве (г, г). Радиус ударной трубки гтах = 15 см, ее длина гтах = 40 см. Известно, что по мере распространения волны по пузырьковой среде происходит ее затухание из-за потерь на сообщение жидкому компоненту кинетиче-
ской энергии и увеличение внутренней энергии газовой фазы в пузырьках при их сжатии. Первая фаза процесса состоит из поглощения пузырьковым кластером ударной волны, при этом волна затухает и ее амплитуда падает (см. рис. 5, б, 6, б). Поглощение ударной волны объясняет сохранение невозмущенной области в центральной части сферического кластера, уже охваченного ее действием (рис. 5, в). Следующая фаза состоит из генерации пузырьковым кластером новой ударной волны и дальнейшего распространения этой ударной волны в пространстве. При взаимодействии с ударной волной пузырьковый кластер генерирует мощный импульс давления в виде волны пузырьковой детонации в случае сферического кластера, с амплитудой, в десятки раз превышающей амплитуду падающей волны (рис. 5, г). Этот эффект является следствием выбора реальной геометрии кластера и различия в скоростях распространения волн в кластере и окружающей его жидкости.
Амплитуда давления в фокусе зависит от ряда параметров — объемной концентрации газовой фазы к0, радиуса кластера КС1 и пузырьков К0 [6]. В дальнейшем аналогичные процессы порождают пульсации в окрестности пятна фокусировки (рис. 5, д, 6, д). Вторая и прочие фокусировки имеют гораздо меньшую амплитуду и быстрее затухают, чем первая, что обусловлено уменьшением пятна фокусировки. Ударная волна, излученная тороидальным облаком пузырьков, имеет ряд особенностей. В процессе схождения к центру симметрии излученной пузырьковым кластером волны (рис. 6, в) имеет место кумулятив-
Рис. 5. Изолинии давления для случая сферического пузырькового кластера. Левый обрез на рисунках соответствует левому торцу кругового цилиндра, нижний — оси симметрии OZ. Параметры течения: ко = 0.01, Яо = 0.1 см, 1С; = 10 см, ЛС; = 5 см. Схема, иллюстрирующая расположение сферического кластера (сечение заштриховано) (а), £ = 60 мкс (б), £ = 120 мкс (в), £ = 180 мкс (г), £ = 250 мкс (б).
Рис. 6. Изолинии давления для случая пузырькового кластера в форме тора. Левый обрез на рисунках соответствует левому торцу кругового цилиндра, нижний — оси симметрии OZ. Параметры течения: ко = 0.001, Ко = 0.1 см, Ісі = 10 см, Л^ог = 6 см, Ксігс = 1 см. Схема, иллюстрирующая расположение тороидального пузырькового кластера (меридиональное сечение тора заштриховано) (а), і = 80 мкс (б), і = 110 мкс (в), і = 130 мкс (г), і = 180 мкс (б).
ный эффект, причем взаимодействие волн, сходящихся к оси симметрии, носит характер нерегулярного маховского отражения (рис. 6, г). Особенности строения ударной волны, переизлученной пузырьковой средой, сохраняющиеся при ее прохождении через чистую жидкость, влияют на форму и структуру имеющего маховскую конфигурацию течения: такое течение располагается в осесимметричной объемной области сложной формы с ярко выраженным ядром высокого давления в отличие от плоского маховского диска в случае ударной волны в гомогенной среде. Интенсивность кумулятивного эффекта находится в зависимости от удельной доли газовой фазы, радиуса пузырьков, амплитуды падающей волны и геометрических параметров тороидального облака пузырьков.
Таким образом, выбор тороидальной формы пузырькового кластера позволяет наблюдать и исследовать такие тонкие волновые эффекты, как маховское отражение ударной волны. Следовательно, результаты численных расчетов показывают, что рассматриваемая модель вполне адекватно описывает динамику волнового поля в кластере, параметры и распределение его излучения. Устойчивость процесса фокусировки волны в кластере проверялась введением специального возмущения в виде жидкой сферы радиусом rdr = 0.5... 1 см, которая размещалась в сферическом кластере на оси z на расстоянии ldr < lei — rdr. Расчет показал, что возмущение, вносимое этой сферой в поле давления, незначительно и быстро исчезает, не оказывая влияния на конечный результат.
Выводы
В данной работе создана и реализована численная модель динамики поля давления и скоростей для аксиально-симметричной задачи о взаимодействии плоской ударной волны со “свободной” пузырьковой системой. Создание численной модели состоит в адаптации явной схемы с направленными разностями и схемы расщепления к решаемой задаче. Показано, что использование предложенных разностных схем для реализации численной модели позволяет вполне адекватно описывать динамику волнового поля в кластере, параметры и распределение его излучения. Приведены сравнительные расчеты по этим разностным схемам, из которых следует, что явная схема более производительна. Показана взаимозависимость граничных условий для входа и выхода плоской ударной волны в расчетной области. Предложена модификация граничного условия для входа ударной волны на основе адиабаты Гюгонио. Создан метод для определения погрешности нахождения значения максимума амплитуды давления в области фокусировки на основе оценки величины изменения энергии волны переизлучения. Приведенные результаты расчетов показывают эффективность предложенных схем для расчета динамики течений с пассивным пузырьковым кластером.
Авторы выражают признательность В.К. Кедринскому за постановку задачи, поддержку и плодотворные дискуссии, В.М. Ковене за помощь в адаптации схемы расщепления к рассматриваемой задаче и обсуждении полученных результатов, Г.И. Дудниковой за постоянное внимание и помощь в работе.
Список литературы
[1] ЗАВТРАК С.Т., Волков И.В. Лазер (Sound Amplification by Stimulated Emission of Radiation) // ЖТФ. 1997. Т. 67, № 4. C. 92-98.
[2] КЕдринский В.К., Вшивков В.А., Дудникова Г.И., Шокин Ю.И. Усиление ударных волн при столкновении и фокусировке в пузырьковых средах // Докл. РАН. 1998. Т. 361, № 1. C. 41-44.
[3] Кедринский В.К., Вшивков В.А., ДудниковА Г.И., Шокин Ю.И. Роль кавитационных эффектов в механизмах разрушения и в крупномасштабных взрывных процессах // Вычисл. технологии. 1997. Т. 2, № 2. C. 63-77.
[4] КЕдринский В.К., Вшивков В.А., ДудниковА Г.И., Шокин Ю.И. Взаимодействие волн в химически активных пузырьковых средах // Докл. РАН. 1996. Т. 349, № 2. C. 185-188.
[5] Кедринский В.К., Вшивков В.А., ДудниковА Г.И., Шокин Ю.И. Взаимодействие и усиление волн пузырьковыми системами // Лаврентьевские чтения по математике, механике и физике: Тез. докл. V Междунар. конф. 18-22 сент. 2000. 63 с.
[6] КЕдринский В.К. Гидродинамика взрыва: эксперимент и модели. Новосибирск: Изд-во СО РАН, 2000.
[7] Kedrinski V.K., Shokin Y.I., Vshivkov V.A., Dudnikova G.I. Shock amplification by bubbly systems with energy release (SABSER) // Proc. 6th Japan — Russian Joint Symp. on Computational Fluid Dynamics. Nagoya, Sept. 21-23, 1998. P. 58-61.
[8] Замараев Ф.Н., КЕдринский В.К., Мейдер Ч. Волны в химически активной пузырьковой среде // ПМТФ. 1990. № 2. C. 20-26.
[9] КовЕня В.М., Лебедев А.С. Модификации метода расщепления для построения экономичных разностных схем // Журн. вычисл. математики и мат. физики. 1994. Т. 34, № 6. C. 886-897.
[10] КовЕня В.М. Некоторые тенденции развития математического моделирования // Вычисл. технологии. 1992. Т. 7, № 2. C 59-71.
[11] Самарский А., Попов Ю. Разностные схемы в газовой динамике. М.: Наука, 1975.
[12] Мудров А.Е. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. Томск: МП “Раско”, 1992.
[13] Березин Ю.А., Вшивков В.А. О критических параметрах ударной волны в плазме // ПМТФ. 1976. № 2. C. 100-106.
[14] КЕдринский В.К., Шокин Ю.И., Вшивков В.А. и др. Генерация ударных волн в жидкости сферическими пузырьковыми кластерами // Докл. РАН. 2001. Т. 381, № 6. C. 773-776.
[15] Вшивков В.А., Дудникова Г.И., Лазарева Г.Г., Шокин Ю.И. Усиление ударных волн пассивными и активными пузырьковыми системами // Восьмой Всерос. съезд по теорет. и прикл. механике: Аннотации докл. Екатеринбург, 2001. C. 168.
Поступила в редакцию 13 мая 2003 г., в переработанном виде —19 июня 2003 г.
