УДК 532.529
ВОЗДЕЙСТВИЕ ударных волн на пузырьковые и пенные
СТРУКТУРЫ В ДВУМЕРНЫХ ОСЕСИММЕТРИЧНЫХ ОБЪЕМАХ
© У. О. Агишева
Башкирский государственный университет Россия, Республика Башкортостан, 450074 г. Уфа, ул. Заки Валиди, 32.
Тел./факс: +7 (347) 229 9616.
E-mail: [email protected]
На основе двумерной с цилиндрической симметрией двухфазной модели газожидкостной смеси проведено численное исследование процессов распространения ударных волн в пузырьковой жидкости и водной пене.
Ключевые слова: ударная волна, пузырьковая и пенная структуры, двухфазная смесь, численное моделирование.
Введение
Двухфазные среды с компонентами различной сжимаемости используются как рабочие жидкости во многих областях техники и промышленности (паровые котлы, химическая и нефтегазодобывающая промышленность). Важность исследования структуры УВ в водной пене обусловлена ее практическим использованием в качестве мобильного барьера, позволяющего рассеивать и поглощать энергию технологических взрывов и применяемого для защиты производственных объектов и транспортных средств. Способность пузырьковых и пенных структур уменьшать все основные характеристики УВ показана в работах Р. И. Нигматулина [1], В. К. Кедринского [2]. Для исследования динамического воздействия на такие среды в условиях, приближенных к натурным экспериментам, становится важным численное моделирование пространственных задач на основе модели двухфазной газожидкостной смеси с использованием реалистических уравнений состояния фаз для получения достоверных значений термодинамических параметров, оказывающих определяющее влияние на формирование волновой картины течения.
В одномерной постановке задача динамики ударных волн (УВ) в пузырьковых жидкостях с небольшим газосодержанием (а20 < 0.06) исследовалась в работе [3]. Для верификации полученных результатов проводилось сравнение с экспериментальными осциллограммами давления [4]. В работе [5] изучено прохождение и отражение плоской нестационарной ударной волны в жидкости с парогазовыми пузырьками в одномерной постановке в диапазоне давлений 1-12 бар и газосодержанием а20 = 0.005-0.01. В [6] численно исследована динамика поля давления, сформированного ударной волной интенсивности 10-120 бар в цилиндре пузырьковой жидкости (а20 = 0.005-0.05), окруженном чистой жидкостью, в двумерной осесиммет-ричной постановке.
Выявлению роли плотности пены для ослабления акустических и ударных волн посвящены экспериментальные работы [7, 8]. Демпфирующая способность пены определяется скоростью звука Су
в среде в зависимости от газосодержания. Она уменьшается с 330 м/с в чистом газе (а: = 0) до ~50 м/с в водной пене (а: = 0.5) [9]. С ростом водо-содержания Су постепенно возрастает до 1500 м/с (а1 = 1). Согласно экспериментальным данным [7] степень защиты зависит также от расстояния до заряда и его массы. В работе [10] исследуется практическое применение пенной оболочки для увеличения энергоэффективности мишени термоядерного синтеза за счет сдерживания УВ, образующейся в дейтерий-тритиевом ядре. В работе [11] проводится численное моделирование процесса распространения ударной волны в пене на основе модели газокапельной среды с использованием метода Годунова в одномерной постановке. В [12] на основе модели двухфазной смеси в односкоростном одно-температурном приближении в одномерной постановке решены задачи, связанные с ударным воздействием на газожидкостную пену с объемным содержанием воды а!0 = 0.05-0.3.
В настоящей работе теоретически исследуются нестационарные процессы взаимодействия ударно-волновых импульсов с пенными и пузырьковыми структурами при давлениях порядка 20-1000 бар с использованием двумерной с цилиндрической симметрией двухфазной модели газожидкостной смеси в однодавленческом, односкоро-стном, двухтемпературном приближениях, базирующейся на положениях механики многофазных сред [1]. Использование уравнения состояния воды [13] позволяет учесть влияние нелинейной сжимаемости среды при расчетах больших давлений в УВ.
Основные уравнения
Система уравнений предлагаемой модели в цилиндрических координатах (где г - ось симметрии, г - радиус) включает в себя [14]:
Уравнения движения и неразрывности смеси:
Г + Э(р + д) = 0, рг +Э(р±^) = 0, дг дг
р
Р
дг дг
— + —
дг дг
+ - = 0,
(1)
r
где p - давление, q - квадратичная искусственная вязкость.
Для описания термодинамических свойств жидкой фазы используется уравнение состояния воды в форме Ми - Грюнайзена [13]:
Р11 = р1 p ](Pii) + Pi(1 W Til), eii = el p '(rii) + eli '(Til),
где тепловая и потенциальная составляющие давления имеют вид:
Pii)(Wii,Tii) = G(piiWiPiiTi, eii) = CviTii, (2)
piP )(Pii ) = A
/ \ 1-ß / / \ / \
P11 1P10 i exp b 1 - (P11 1P0 i - K IP111 1P i
pi
1(1 '(Г11)
Г121
dP 11
здесь А, К, Ь, в, - постоянные.
Для газа принимается уравнение состояния совершенного газа
Р21 = РпЯТ21 = Рп (7- 1)сК2Т21 . (3)
Законы сохранения массы каждой фазы и всей смеси имеют вид:
р10 0 _ а11 р0 р20а20 _ а20 р0
Pii
P
P21
P
p10a10 + p20a20
P p
(4)
Р11 Р21 В уравнениях 2, г - проекции скоростей, р = 1/V - средняя плотность смеси, р10, р20, р11, р21, а10, а20, а11, а21- начальные и текущие плотности и объемного содержания жидкой и газовой фаз, Т11, Т21 - текущие температуры жидкости и газа, р 0 -плотность жидкости при нормальных условиях, ск1, ск2 - теплоемкости жидкости и газа при постоянном объеме, Я - универсальная газовая постоянная, у -показатель адиабаты газа.
В предположении, что внутренняя энергии смеси в случае отсутствия фазовых переходов является аддитивной функцией по массе [1], воспользуемся утверждением [9], что при ударном сжатии смеси каждая фаза сжимается по индивидуальной ударной адиабате:
- ударная адиабата жидкости
/ V
Pii =
2
Р(р) p11
G(pn)
•2(e¡p
"e10)P11 + P10
Pii Pió
-1
1 + -
2
Pii Pió
Г(Рц) - ударная адиабата газа
Р21 = Р21 -(Г+!) + Р20 \у-1) Р20 Р20 \г+1) + Р21 -(Г-1) . Из условия равенства давлений = Р11 = Р21 определяется температура:
(5)
(6)
фаз
Tii =
p - P1(p)
P11 r(P11)^F 1 T2
_1L = (7)
20
P _Pl0_ p20 P21
Метод расчета
При численном моделировании на основе метода сквозного счета Уилкинса [14] расчетная область делится на четырехугольники сеткой / - к, которая движется вместе со средой. Обозначим координаты центров четырехугольников следующим образом: 1 = (/+1/2, к+1/2); 2 = (/-1/2, к+1/2); 3 = (/-1/2, к-1/2); 4 = (/+1/2, к-1/2). Для вершин четырехугольников примем: 1 = (/, к); 2 = (/+1, к); 3 = (/+1, к+1); 4 = (/, к+1).
к+1 к к-1 2 1
0
/ i r севдо нчейк i / /
\ \
ve 4 ÍÍI 3 " V
) IX V b \ .II/' )
\ \ } í 3 2 4
-O" 'VV'
\ \ / /
2' ./-/ j j+1
m
Рис. 1 Схема расчетной области
Для определения граничных условий вводятся псевдоячейки. Масса псевдоячеек у жесткой стенки задается зеркальным отображением через границу. В случае свободной границы плотности и давления в псевдоячейках соответствуют газу при нормальных условиях.
В начальный момент времени в центрах ячеек задается распределение давления, вычисляются масса и объем смеси и каждой фазы. Из условия сохранения массы ячейки вычисляется относительный объем, соответствующий четырехугольнику 1:
V"
Уз (%
k" + r" + r" )• A"
+
где A" и Ab - площади треугольников:
< A")" = 1
r21+
(A")"=2
z" (гз" - Г"" ) + z" (г" - r" ) ■"" (2 - Гз")
(гз" - Г"" ) + z" (r" - r")
2\'з _'4
4" (г" - Гз
( a)1 =( A")" +(Ab)".
ró +
b
После вычисления объемов находится масса всей ячейки и каждой фазы:
Ап
т =Р0\ -^п , т11 = Р^ХХУХ , т21 = Р210.1У\ •
Компоненты скорости вычисляются в узле сетки ;, к на половинном шаге по времени ¿п+1/2:
.п+12 .п-12 , , г; к = г; к +1
V/2$к )•
[рП (гп - гш)+ рп к, - п"у)+ рп (г1п - гп)], г;+/2 = г;/2+(&п/2.ф"к )•
[п I п п I п п I п п
Р1 II - г III ) + р 2 1г III - г IV ) + Р 3 1г I - г II Л
Ч>/ к =-
Рс Ап
Уп
у 1
Р0 Ап
Р0 Ап
Уп
Уп
у 2
Р Ап
Уп
где координаты узлов I = /, к-1); II = (/+1, к); III = /, к+1); IV = (/-1, к).
Далее определяются новые координаты узлов сетки:
п +1 п п + . п + %
/ = г;,к + г, /2 & /2,
п +1 п п + % А п + %
г;;:1 = ;+г^ &
Для нахождения термодинамических параметров смеси воспользуемся методом, предложенным В. Ф. Куропатенко для одномерных задач [15], и обобщим его на двумерный случай. Будем отличать ячейки, в которых находится разрыв решения, от ячеек с гладким решением, определяя знак разности скоростей в соседних ячейках. Положим, что в ячейке ;, к находится ударная волна, если выполнен° условие г; к+1 - г; к < 0 либ° г/+!,к - г; к < ^ и в ячейке находится волна разгрузки, если
'; к+1
■Г; к > 0 либо ^ +и
Л,к >0.
Расчет новых значений давления и плотностей фаз на шаге 1п+1 производится методом поиска корня из уравнения сохранения массы смеси (4), где Ри, Р21 берутся из уравнений ударной адиабаты (5)-(6) или адиабаты Пуассона, а значения плотности р10, р20 с предыдущего шага по времени Л Далее определяется температура фаз (7). Искусственная вязкость (1) имеет вид [14]:
Я." =(с 02Р0 ; ,к (Р; к/Р; )2 Ап,к Р]к к ) 1 На конечном этапе вычислительного цикла находятся объемные содержания фаз:
ап;,к = (Рп,к -Р\;к)/(Р2п;,к -Рр;,к), к =1 - а;, к
Результаты и их обсуждение
а) Ударно-волновое воздействие на пузырьковую жидкость
Тестовые расчеты для проверки построенной модели проводились в условиях экспериментов
работы [4]. В камере высокого давления в ударной трубе (рис. 2) генерируются ударные волны. Камера низкого давления (Ь0 = 3.55 м, Ь\ = 0 м) заполнена пузырьковой жидкостью с содержанием газа от 0.005 до 0.06. С помощью датчиков давления на стенках трубы измеряется амплитуда падающих и отраженных ударных волн.
Рис. 2 Схема ударной трубы.
9 10 11 12 13 мс Рис. 3 Сравнение с экспериментальными данными [4].
На рис. 3 представлено сравнение экспериментальной осциллограммы давления на расстоянии 0.535 м от торца трубы и расчетной кривой, полученной по двумерной модели в квазиодномерной постановке для пузырьковой жидкости с содержанием газа а20 = 0.04. Расчеты показали, что амплитуда падающей и отраженной волны, а также время прихода отраженной волны удовлетворительно согласуются с экспериментом. Низкочастотные осцилляции в экспериментальных данных при отражении, являющиеся следствием конструктивных особенностей ударной трубы, не учитывались при численном моделировании.
Для изучения влияния двумерных эффектов на динамику УВ в жидкости моделировалась следующая задача: в цилиндре с жесткими стенками (рис. 4, г = 0.5 м, г = 1 м) в начальный момент времени в зоне 1 размером 6г = йг = 0.2 м задается импульс давления р0тах = 1000 бар (рис. 5), аппроксимируемый гладкой экспоненциальной функцией вида
р(г, г) = е
(10)
+
+
1
+
+
4 д 4
правая граница
граница симметрии 3
Рис. 4 Расчетная схема задачи
в начальный момент времени.
p, бар 1000
500
0.5 r, м
0.25
z, м
10
Рис. 5 Начальный импульс давления.
Вся область заполнена либо чистой, либо пузырьковой жидкостью с содержанием газа а20 = 0.04. На рис. 6 показаны расчетные распределения давления в различные моменты времени в цилиндре, заполненном водой (а2 = 0). При воздействии первоначального импульса давления с течением времени формируется радиально распространяющаяся УВ, направленная к границам области. При этом в области 1 возникает зона кавитации, в которой содержание газа достигает а2 = 0.6. После частичного отражения импульса от боковой границы УВ движется к оси симметрии, обходя зону кавитации 1 и
р, бар ^ ^ 400 " ' '
300
200
100
0
0.5
r, м
догоняя волну, распространяющуюся в направлении правой границы (см. рис. 6 при t = 0.4 мс). При достижении оси симметрии наблюдается усиление амплитуды давления до 400 бар (см. рис. 6 при t = 0.6 мс), и волна продолжает движение к правой границе с уменьшением амплитуды давления (см. рис. 6 при t = 0.7 мс, 0.8 мс). У правой границы амплитуда УВ возрастает до 426 бар ближе к оси симметрии г (см. рис. 6 при t = 0.85 мс). Аналогичная задача в случае симметричной расчетной области размером г = г = 0.5 м рассмотрена в [16], где, в отличие от настоящей работы, наблюдалась кумуляция импульса давления на пересечении оси симметрии и правой границы.
На следующем этапе моделировалось распространение УВ в цилиндре, заполненном пузырьковой жидкостью (а20 = 0.04). Первоначальный импульс формирует радиальную УВ затухающей периодической структуры (см. рис. 7 при t = 1, 3 мс), но значительно меньшей амплитуды, чем при движении в чистой жидкости. При отражении от жестких границ амплитуда УВ достигает 7.2 бар. Далее отраженная УВ обходит возникшую в зоне 1 газовую полость (а2 = 1) и движется по сжатой среде (а2 = 0.025) к оси симметрии г, где ее интенсивность увеличивается до 9.5 бар (см. рис. 7 при t = 4, 5 мс). По мере приближения к правой границе волновой импульс приобретает форму, близкую к одномерной (см. рис. 7 при t = 7 мс). Далее происходит отражение сформировавшейся УВ от правой границы области, и отраженная УВ амплитуды 13.2 бар принимает одномерный вид (см. рис. 7 при t = 9 мс). В случае симметричной расчетной области (г = г = 0.5 м) в [16] наблюдалась кумуляция УВ на оси симметрии возле правой стенки с амплитудой давления ~ 30 бар, и, в отличие от настоящих расчетов, квазиодномерная УВ не формировалась.
0.4 мс
Шшшш
0.5
0.25
0.5
r, м
0.5
r, м
1 0
0.5
r, м
p, бар 400
300
200
100
0
0.8 мс
0.5
1 0
0.5
r, м
p, бар 400-'
300-'
200-'
100-''
0-U
0.85 мс |
10
0.5
r, м
10
Рис. 6 Динамика волнового импульса в сечении цилиндра, заполненного чистой жидкостью.
0
z. м
0.5
0.5
z, м
0.5
0.5
г, м
0.5
г, м
0.5
г, м
1 0
1 0
1 0
0.5
г, м
0.5
1 0
0.5
г, м
1 0
Р, бар
15Т" 100 <
10
0.5
г, м
Рис. 7 Динамика волнового импульса в сечении цилиндра, заполненного пузырьковой жидкостью
б) Взаимодействие ударной волны с барьером из водной пены
Рис. 8 Сравнение с экспериментальными данными [8]
При решении пространственных задач о взаимодействии УВ с пенным слоем, как и в расчетах в одномерной постановке [12], тестирование двумерной модели проводилось по данным экспериментальной работы [8]. Исследование проводилось в ударной трубе (рис. 2), в которой камера низкого давления состояла из двух секций: первая заполнена воздухом (Ь0 = 3 м), вторая - пеной (Ь1 = 0.363 м). Начальное содержание воды в пене а10 = 0.2. На рис. 8 приведено сравнение расчетной кривой и осциллограммы давления, полученной на датчике, расположенном в газе на расстоянии 0.093 м от границы пены. Расчетная амплитуда проходящей волны и время прихода отраженной
Р, бар 0^
2 мс
0.5
Р, бар
4т-"
1
г, м
УВ удовлетворительно согласуются с экспериментом. Интенсивность отраженной УВ имеет завышение, что объясняется неучетом в модели вязкости пены. После завершения релаксационных процессов к моменту времени t = 2 мс амплитуда волны в расчете и эксперименте выравнивается.
Для исследования эффективности пенной защиты решается следующая задача: в цилиндрической области (рис. 4, г = £ = 1 м) в начальный момент времени в зоне 1 размером 62 = 6г = 0.2 м задается импульс давления с амплитудой Р0тах = 20 бар, аппроксимируемый функцией вида (10). Область 2 заполнена чистым газом (азот). В слое 3 толщиной 0.2 м располагается слой пены (а10 = 0.1). Боковая и правая граница полагаются свободными, оси 2 и г - границы симметрии.
В результате распада разрыва на границе скачка давления возникает УВ, имеющая затухающую периодическую структуру, которая распространяется к границам области (см. рис. 9 при t = 2 мс). При t = 3 мс волновой импульс достигает границы пены. Происходит отражение волны сжатия от более плотной по сравнению с газом поверхности пены, что приводит к усилению амплитуды давления и к существенному уменьшению скорости распространения УВ в слое пены. Взаимодействие отраженной от пенного слоя волны со встречными импульсами приводит к усложнению волновой картины течения (см. рис. 9 при t = 4, 5 мс): пенный слой блокирует распространение УВ в направлении оси 2.
0.5
3 мс
0.5
Р, бар
4 г "
" 1
г, м
0.5
5 мс
0.5
1
г, м
10 1 0 ^ Рис. 9 Динамика волнового импульса в области, содержащей слой пены
1 0
м
м
2
0
0
м
Выводы
На основе двумерной с цилиндрической симметрией модели газожидкостной среды проведены численные исследования процесса формирования и распространения УВ в воде, в пузырьковых и пенных структурах. Получено удовлетворительное согласование тестовых расчетов с экспериментальными данными. Показано, что в чистой жидкости формирование УВ сопровождается образованием зоны кавитации и после частичного отражения импульса от внешней жесткой границы формируются два фронта УВ, взаимодействие которых приводит к кумуляции УВ на оси симметрии у правой стенки сосуда. При распространении волнового импульса в пузырьковой жидкости формируется существенно неодномерная картина течения в центральной области, которая в дальнейшем трансформируется в квазиодномерную при приближении УВ к правой границе области. Получены оценки эффективности пространственного взаимодействия волнового импульса с пенной преградой в условиях, приближенных к реальным задачам защиты промышленных объектов от технологических взрывов. Расчеты по двумерной модели показали более высокую эффективность демпфирующих свойств пены по сравнению с одномерным приближением.
Автор работы выражает благодарность научному руководителю д.ф.-м.н. Раисе Хакимовне Бо-лотновой за постановку задачи и полезные обсуждения в ходе решения.
Работа выполнена при финансовой поддержке грантов РФФИ (11-01-97004_р_поволжье и 11-0100171 -а), и Программы фундаментальных исследований ОЭММПУ РАН (ОЭ-13).
ЛИТЕРАТУРА
1. Нигматулин Р. И. Динамика многофазных сред. Ч. 1. М.: Наука, 1987. 464 с.
2. Кедринский В. К. Гидродинамика взрыва. Эксперимент и модели. Новосибирск: Изд-во СО РАН. 2000. 435 с.
3. Болотнова Р. Х., Галимзянов М. Н., Агишева У. О. Моделирование процессов взаимодействия сильных ударных волн в газожидкостных смесях // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. 2011. №2. С. 3-14.
4. Сычев А. И. Сильные ударные волны в пузырьковых средах // Журнал технической физики. 2010. Т. 80. №6. С. 31-35.
5. Ганиев О. Р., Хабеев Н. С. Ударные волны в жидкости с пузырями, содержащими испаряющиеся капли сжиженного газа // Изв. РАН. МЖГ. 2002. №3. С. 98-107.
6. Кедринский В. К., Вшивков В. А., Лазарева Г. Г. Формирование и усиление ударных волн в пузырьковом «шнуре» // ПМТФ. 2005. Т. 46. №5. С. 46-52.
7. Raspet R., Griffiths S. K. The reduction of blast noise with aqueous foam // Journal of Acoustical Society of America.
8.
1983. №74. P. 1757-1763.
Britan A., Ben-Dor G., Shapiro H., Liverts M., Shreiber I. Drainage effects on shock wave propagating through aqueous foams // Colloids and Surfaces A: Physicochem. Eng. Aspects. 2007. № 309. P. 137-150.
9. Агишева У. О., Болотнова Р. Х., Бузина В. А., Галимзя-нов М. Н. Параметрический анализ режимов ударно-волнового воздействия на газожидкостные среды // Известия Российской академии наук. Механика жидкости и газа. 2013. №2. C. 15-28.
10. Norimatsu T., Nagai K., Takaki T., Yamanaka T. Issues in capsule fabrication and injection into a wet-walled IFE reactor // Fusion Engineering and Design. 2001. Vol. 55. Is. 4. P. 387-396.
11. Васильев Е. И., Митичкин С. Ю., Тестов В. Г., Ху Х. Динамика давления при ударном нагружении газожидкостных пен // Журнал технической физики. 1998. Т. 68. №7. С. 19-23.
12. Болотнова Р. Х., Агишева У. О. Особенности распространения ударных волн в водных пенах с неоднородной плотностью // Сборник трудов Института механики УНЦ РАН. Выпуск 9. 2012. С. 41-46.
13. Нигматулин Р. И., Болотнова Р. Х. Широкодиапазонное уравнение состояния воды и пара. Упрощенная форма // Теплофизика высоких температур. 2011. Т. 49. № 2. С. 310-313.
14. Олдер Б., Фернбах С., Ротенберг М. Вычислительные методы в гидродинамике. М.: Мир. 1967. 384 с.
15. Куропатенко В. Ф., Мустафин В. К. Методика расчета нестационарных течений в многокомпонентных неравновесных смесях веществ // Вестник Челябинского университета. 1997. №1. С. 97-102.
16. Bolotnova R. Kh., Agisheva U. O. Numerical investigation of two-dimensional axisymmetric bubbly liquid shock wave flow in bounded volume // Fluxes and Structures in Fluids: international conference, 25-28 June, 2013. Saint Petersburg. Submitted for publication.
Поступила в редакцию 28.05.2013 г.