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

Механика

Вестник Нижегородского университета им. Н.И. Лобачевского, 2011, № 6 (1), с. 163-168

УДК 532.593

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПРОСТРАНСТВЕННОГО ВЗАИМОДЕЙСТВИЯ ВОЗДУШНОЙ УДАРНОЙ ВОЛНЫ С ПРОНИЦАЕМЫМ ГРАНУЛИРОВАННЫМ СЛОЕМ

© 2011 г. Е.Г. Глазова, А.В. Кочетков, И.А. Ходыкина

НИИ механики Нижегородского госуниверситета им. Н.И. Лобачевского

kochetkov@mech.unn.ru

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

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

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

Насыпные гранулированные слои из различных твердых материалов, являясь проницаемыми для газа, используются для уменьшения амплитуды проходящих ударных волн [1]. Предложены упрощенные модели, описывающие эти процессы [2]. Эти модели включают в себя ряд параметров, сильно влияющих на характеристики отраженных и проходящих ударных волн, значения которых определяются либо из экспериментов, либо из весьма общих рассуждений. Реальные физические процессы, происходящие при взаимодействии ударной волны с подобными преградами, отличаются высокой степенью сложности [3,4] и требуют тщательного изучения. В данной работе с помощью вычислительного комплекса STAR-CCM+ моделируется пространственное газодинамическое взаимодействие плоской воздушной ударной волны с насыпным слоем металлических шариков.

1. Постановка задачи.

Основные уравнения

Вводится декартова система координат Oxyz, ось Oz которой направлена перпендикулярно насыпному слою толщиной Н. Начало координат располагается на свободной поверхности насыпного слоя. Тогда z = -H соответствует твердой поверхности, на которую опирается слой.

Выделяется элементарный канал квадратного сечения г х г вдоль оси Oz , боковые границы которого являются плоскостями симметрии,

г - радиус шариков. Область определения задачи включает в себя объем канала, ограниченный плоскостями симметрии x = 0 , x = г , у = 0, у = г , поверхностями шариков и плоскостями г = L, г = — H (рис.1). Параметры в набегающей плоской ударной волне (УВ) бесконечной длительности задаются в области 0 < г < L . Расстояние L выбирается из таких соображений, чтобы в процессе моделирования отраженная от слоя ударная волна не исказила условия на искусственно поставленной границе г = L .

Рис. 1

Движение газа (воздуха) описывается системой уравнений [5], включающей в себя уравнение неразрывности, уравнения Навье-Стокса и уравнение закона сохранения энергии:

— + V(pV) = О дt

дu др д

P---- =--------1---

дt дх дх

Dv др д

P-------- =---------1----

дt ду дх

Dw др д

P-------- =---------1-----

дt Dz дх

и21J-f*V

ду

Du Dv

Н| д + ~дХ

( Du DvI ду дх j

(Du Dw

aZ+~дх

_д_

ду

д

+------

ду

(i Dv і t 2 —

І дУ J

Dv Dw

-—divV З

Dz

д

Dz

Du Dw

н| +~дх

( Dv Dw^ Dz ду

д_

Dz

H2 ї У ^^divV

DE DE DE DE і

P| — + u — + v— + w— I + pdivV = div(k-VT) + Ф, Dt дх ду Dz J

д

д

+

+

+

где р - плотность, V = (и,V, w) - скорость, ц -коэффициент вязкости, Е - удельная полная энергия, р - давление, к - коэффициент теплопроводности, Т - температура,

» , | ди У (ду ^ (дwУ 1 ( ди ду^

Ф = 1их) +ЬуI +1& I+ 2 Ьу+¥ I+

1 ( Dv Dw +----I------1-----

2 iDz ду

1

+ — 2

Du Dw

---------1---------

Dz дх

- -3- H(divV )2

- диссипативная функция. Для воздуха при умеренной температуре и давлении термодинамические параметры связаны уравнением состояния идеального газа:

m

Р = 77 RPT,

M

где m - масса, M - молярная масса, R - газовая постоянная, R = 0.278 кДж/(кг-К).

Для описания турбулентных движений используется стандартная к-є модель. В качестве краевых условий на плоскостях симметрии, на поверхности z = - H и на поверхностях шариков за-дается условие непротекания, т.е. нормальная скорость газа равна нулю (vn = 0). На поверхностях неподвижных шариков по температурному режиму задаются адиабатические или изотермические условия. На искусственной границе z = L задаются параметры за фронтом набегающей ударной волны. При t = 0 газ внутри слоя шариков покоится при условиях нормальной атмосферы. В области z > 0 в качестве начальных условий задаются постоянные параметры за фронтом набегающей ударной волны.

Решение задачи проводится с помощью пакета STAR CCM+ версия 5.0б [б]. При задании физических моделей используется нестационарный решатель на основе неявного метода конечных объемов [5, 7], аппроксимирующий

полные уравнения Навье-Стокса для сжимаемого вязкого теплопроводного газа с первым и вторым порядком точности. Начальный временной шаг = 0.3 мкс. Для области Ь = 9Н, что гарантировало отсутствие возмущений с искусственной границы г = Ь в интересующем интервале времени, типовая пространственная сетка из многогранников содержала 11 395 ячеек, 51 845 внутренних граней, 51 128 вершин.

В начальный момент времени t = 0 параметры невозмущенного газа имеют следующие значения: р = р0 = 0.1 МПа, Т = Т0 = 300 К, V = 0. За фронтом плоской ударной волны: рув = = 0.1805 Мпа, иув = уув = 0, wув = -149.2 м/с, Туе = = 355 К. Динамическая вязкость воздуха ц = = 1.855 • 10—5 Па-с. Коэффициент теплопроводности к = 0.026Вт/(м-К), г = 0.52мм, Н = 17мм. Геометрические размеры и параметры ударной волны выбраны в соответствии с условиями проведения экспериментов [1], где размеры соответствуют объемной концентрации шариков а 2 = 0.61.

В поисках упаковки с заданной объемной концентрацией твердой фазы были рассмотрены несколько видов упаковок. Упаковка первого вида (рис. 2а), в которой центры соседних шариков образуют кубическую ячейку, оказалась недостаточно плотной, ее объемная концентрация твердой фазы составила а2 = 0.52. Упаковка второго вида (рис. 2б), когда центры соседних шариков образуют тетраэдр, позволяет получить максимальное значение а 2 = 0.74 . За счет увеличения расстояния между центрами шариков стало возможным смоделировать заявленное в эксперименте [1] значение а 2 = 0.61.

Ниже приводятся результаты расчетов для второго типа упаковки шариков.

2

2

б)

Рис. 2

Рис. 3

2. Результаты численного моделирования

При воздействии ударной волны на слой шариков происходит отражение ударной волны, распространяющейся в обратном направлении по отношению к падающей волне. Внутрь слоя по поровому газу распространяются волны сжатия, обтекающие каждый шарик как жесткое тело. Внутреннее течение является весьма сложным. Его параметры определяют процессы многократного отражения волн и взаимодействия с турбулентными образованиями за шариками. В результате этих взаимодействий происходит относительно медленное нарастание давления на жесткую стенку, на которой располагается слой. Важные значения имеют вязкость среды и теплообмен между поровым газом и металлическими шариками. С целью оценки значимости этих эффектов проведены расчеты с различными значениями вязкости и различными условиями теплообмена. Кроме того, интерес вызывает влияние точности применяемой численной схемы на параметры процессов. Численные решения сопоставляются с экспериментальными данными [1].

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

слоя ударной волны. Наблюдается некоторое повышение амплитуды отраженной ударной волны с течением времени, вызванное увеличением давления порового газа в преграде. Начало отсчета времени в эксперименте связано с моментом разрыва мембраны между камерами высокого и низкого давления. Здесь же красным цветом показано решение для а2 = 0.61, полученное по численной схеме второго порядка точности. Наблюдается соответствие численных и экспериментальных результатов.

На рис. 4а показаны зависимости от времени осредненного по площади сечения канала давления газа на тыльной стенке г = —Н за пористым слоем толщиной Н = 17 мм при а2 = 0.61. Синим цветом отмечены экспериментальные данные [1], зеленым и красным — численные решения с адиабатическими условиями на шариках, при динамической вязкости газа ц = = 1.855 Па-с, с первым и вторым порядком точности соответственно. На рис. 4б: красным и желтым цветом — численные решения с адиабатическими и изотермическими условиями на шариках соответственно, при динамической вязкости газа ц = 1.855 Па-с, со вторым порядком точности. Сиреневым цветом отмечено численное решение с увеличенной вязкостью газа, равной ц = 6.0 Па-с, с адиабатическими условиями на шариках.

а

//* * /

У/

У" Г

0 с 1 0 2 0 3 о 1, мс б Рис. 4 4 0 5 0 6 0

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

Далее приведены расчеты для слоя толщиной Н = 37 мм и объемной концентрацией твердой фазы, равной а2 = 0.61. На рис. 5а синим цветом отмечены экспериментальные данные [1], красным и зеленым — численные решения с изотермическими условиями на шариках, динамической вязкостью газа ц = 1.855 Па-с, с первым и вторым порядком точности соответственно. На рис. 5б фиолетовым и зеленым цветом отмечены численные решения с адиабатическими и изотермическими условиями на шариках соответственно, динамической вязкостью газа ц = 1.855 Па-с, со вторым порядком точности. Желтым и голубым цветом отмечено чис-

ленное решение с увеличенной вязкостью газа, равной ц = 6 Па-с, адиабатическими и изотермическими условиями на шариках соответственно.

Изотермическое условие на шариках, т.е. максимальный учет теплообмена между газом и металлом, а также увеличение вязкости приближает численное решение к эксперименту [1]. В этом случае численные кривые не только качественно, но и количественно соответствуют экспериментальным.

С целью выяснения возможных причин недостаточно точного соответствия расчетов и эксперимента для Н = 17 мм проведены исследования процессов для нескольких значений объемной концентрации шариков. На рис. 6 представлены расчеты с различной начальной плотностью упаковки. Синим цветом отмечен эксперимент [1], голубым - расчет с объемной концентрацией, равной а2 = 0.61, зеленым - а2 = = 0.67, желтым - а2 = 0.71, розовым - а2 = 0.74. Из приведенных результатов следует, что соответствие численных результатов экспериментальным улучшается при более плотной упаковке шариков, чем начальная. Этот эффект объясняется переукладкой самих шариков в

И Л |1 i1 1

/

// / // /

J /

У

О 0,2 0,4 0,6 0,8 1 1,2 1,4

t, мс

а

: - -——

—

/ / / / s'^

///

///

V \

О 0,2 0,4 0,6 0,8 1 1,2 1,4

t. мс

б

Рис. 5

t. мс

Рис. 6

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

3. Выводы

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

же переукладка самих частиц под действием давления. Соответствие численных и экспериментальных данных свидетельствует о корректности используемых математических и численных моделей.

Авторы выражают благодарность Николаевой А.С. за помощь в проведении расчетов на STAR CCM+ 5.06.

Работа выполнена при частичном финансировании РФФИ (грант № 09-08-00711а), Программы

поддержки ведущих научных школ России (грант № НШ-4807.2010.8) и Федеральной целевой программы «Научные и научно-педагогические кадры инновационной России на 2009—2013 годы» (ГК № 16.740.11.0087).

Список литературы

1. Ben-Dor G., Britan A., Elperin T., Igra O., Jiang J. P. Mechanism of compressive stress formation during weak shock waves impact with granular materials // Experiments in Fluids. 1997. V. 22. P. 507-518.

2. Губайдулин А.А., Дудко Д.Н., Урманчеев С.Ф. Моделирование взаимодействия воздушной ударной волны с пористым экраном // ФГВ. 2000. Т. 36. № 4. С. 87-96.

3. Нигматулин Р.И. Динамика многофазных сред. Т. 1, 2. М: Наука, 1987.

4. Ударно-волновые процессы в двухкомпонентных и двухфазных средах / Киселев С.П., Руев Г.А., Трунев А.П., Фомин В.М., Шавалиев М.Ш. Новосибирск: ВО «Наука», Сибирская издательская фирма, 1992.

5. Флетчер К. Вычислительные методы в динамике жидкостей. Т. 1, 2. М.: Мир, 1991.

6. Документация, сопровождающая вычислительный комплекс STAR CCM+ for Windows 4.02.011.

7. Пейре Р., Тейлор Т. Вычислительные методы в задачах механики жидкости Л: Гидрометеоиздат, 1986.

NUMERICAL SIMULATION OF THE SPATIAL INTERACTION OF AN AIR SHOCK WAVE WITH A PERMEABLE GRANULAR LAYER

E. G. Glazova, A. V. Kochetkov, I. A Khodykina

The processes of a shock wave passage through bulk gas-permeable layers of metal balls are studied by numerical simulation methods. The main parameters of transmitted and reflected shock waves (pressure, density, speed, temperature) are analyzed. We reveal the most significant factors determining these processes, including gas viscosity, heat transfer between gas and metal, and packing density.

Keywords: numerical simulation, shock wave, gas, gas-permeable barrier, viscosity.