CC BY
38
Мезоскопические течения в неоднородном газе
Д.А. Медведев, А.П. Ершов, Ю.В. Янилкин1, Е.С. Гаврилова1
Институт гидродинамики им. М.А. Лаврентьева СО РАН, Новосибирск, 630090, Россия 1 РФЯЦ Всероссийский научно-исследовательский институт экспериментальной физики, Саров, 607190, Россия
Рассматривается динамика газовой среды, в которой распределены включения другого газа. При импульсном воздействии на систему, например при прохождении ударной волны, возникают мезоскопические течения, приводящие к искажению формы включений и перемешиванию. Результаты применяются к анализу процесса детонации гетерогенных взрывчатых веществ.
1. Введение
Во многих задачах механики дисперсных смесей существенно гидродинамическое взаимодействие между фазами. Ряд работ посвящен поведению газовых включений в несущей среде, представляющей собой более тяжелый или более легкий газ. При прохождении ударной волны граница включения деформируется и генерируется неоднородное поле скорости. Дальнейшая деформация провоцируется обтеканием неоднородности. При этом могут проявляться гидродинамические неустойчивости, резко ускоряющие смешение. Поэтому такие задачи относятся к достаточно сложным проблемам гидродинамики, на которых тестировались многие гидродинамические коды. Как правило, рассматривалось уединенное включение — сфера (например в виде мыльного пузыря [1]) или цилиндр (например струя [2]). В данной работе рассматривается более сложный случай «дисперсного газа», в котором объемная доля обоих газов в системе не мала и существенно взаимодействие включений не только с несущим газом, но и между собой.
Практический интерес связан, в частности, с реагирующими системами. Из малых объемов горючего газа,
диспергированных в воздухе, за ударной волной может достаточно быстро сформироваться квазигомогенная смесь, способная к детонации. И напротив, при смешении с инертной или содержащей ингибиторы добавкой газ может стать невосприимчивым к детонации.
Другой пример, который непосредственно привел к постановке данной работы — взрывчатые вещества. В литом составе тротил/гексоген частицы гексогена распределены в твердой матрице тротила. Эти вещества реагируют за фронтом ударной волны в течение ~0.1 мкс в основном раздельно [3, 4]. Далее «пузыри» газообразных продуктов детонации гексогена взаимодействуют с продуктами детонации тротила в течение более долгого времени (микросекунд), когда может быть существенным перемешивание. Именно составы тротил/гексоген представляют собой лучшее «сырье» для детонационного синтеза алмаза [5, 6]. Мезопроцес-сы, сопровождающие взаимодействие включений с окружением, наблюдались при электрической диагностике детонации [7].
В данной работе мезоскопические течения, вызванные присутствием в газе неоднородностей, исследовались численно. Часть предварительных результатов
© Медведев Д.А., Ершов А.П., Янилкин Ю.В., Гаврилова Е.С., 2004
изложена в докладах [8, 9], но в основном ниже излагаются новые данные, с более детальными расчетами и более сложными моделями. Хотя некоторые оценки «привязаны» к продуктам детонации состава тротил/ гексоген, результаты относятся практически к любой системе, в которой мезообъемы одного газа диспергированы в другом газе. Требуется только, чтобы число Рейнольдса Re при обтекании включений было хотя бы умеренно большим (например порядка 1000). Выполнение этого условия вполне типично. Например, для «пузыря» размером 1 см в воздухе достаточна скорость обтекания 1.5 м/с.
Дополнительной целью работы была оценка возможностей развитых в последнее время решеточных методов газодинамики и сравнение их с классическими разностными методами.
2. Предварительные оценки
Взаимодействию включений и несущего газа и перемешиванию могут способствовать гидродинамические неустойчивости: Рэлея-Тейлора, Рихтмайера-Мешкова и Кельвина-Гельмгольца. К настоящему времени эти неустойчивости изучены в сравнительно простых ситуациях.
Неустойчивости Рэлея-Тейлора посвящено большое число как экспериментальных, так и расчетных работ (см., например, [10, 11]). Анализ литературных данных и сравнение их с результатами трехмерных расчетов можно найти в [12]. Для сравнительно небольших чисел Атвуда А = (р2 - рОДр2 + Рх) толщина слоя смешения при развитой неустойчивости Рэлея-Тейлора h ~ ~ 0.06А^2, где g — ускорение силы тяжести. Согласно [2, 13] характерная скорость V, до которой разгоняется включение на первом этапе (за время 1-2 пробегов звуковой волны по включению диаметра й), порядка (1 - А)У, где V — массовая скорость газа за ударной волной. Тогда характерная величина g ~ (1 - А)Ус/й, а характерное смещение границы gt2 /2 можно оценить как (1 - А)Уй/с. Имеем Н/й ~ 0.12А(1 - А)У/с <
< 0.03У/с. Для слабых ударных волн и тем более для малых либо близких к единице чисел Атвуда перемешивание Рэлея-Тейлора на этапе начального разгона несущественно.
Далее «пузырь» разгоняется обтекающим его потоком при относительной скорости и = АУ (скорость проскальзывания). Считая, что характерное ускорение на этом этапе g ~ и 2/й, за характерное время обтекания tу ~ й/и имеем Н/й ~ 0.06А, то есть этот эффект также слаб, особенно при малых А (для пары тротил/ гексоген А = 0.055). Этот вывод сохраняется и при использовании эмпирических зависимостей, применяемых для случая переменного ускорения [11]), в том числе и для больших времен, так как ускорение со временем убывает.
Неустойчивость Рихтмайера-Мешкова возникает при резком ускорении границ раздела в ударной волне. Здесь закон роста возмущений не экспоненциальный, а линейный, и поэтому начальные условия полностью не «забываются». При длине волны X ее амплитуда а растет согласно ёа/= (2лАу/Х)а0, где а0 — начальная амплитуда, V ~ (1 - А) V — скачок скорости границы в ударной волне. Закон роста можно записать и в виде ёа/ёг1 = 2п(1 -А)иа0/X. За время обтекания
= ^и возмущение вырастет до амплитуды а = а0(1 + 2п(1 - А^/X). Если положить X/2п ~ ^2, а0 ~ ^2, то за время форма заметно исказится. Однако не вполне корректно считать «начальным возмущением» исходную форму включения: это явно выходит за границы применимости понятия возмущения. Начальная фаза взаимодействия ударной волны с включением отнюдь не сводится к возбуждению неустойчивости Рихтмайера-Мешкова. Более того, из расчетов, приведенных ниже, следует ее сравнительно малая важность. Правильнее рассматривать этот этап как «генератор» начальных условий для последующего течения, не выделяя произвольно отдельных эффектов из всей их совокупности.
Если же рассматривать коротковолновые возмущения, X/2п << d, но достаточно гладкую форму поверхности а0 «X/2п, то скорость роста ёа/<< и, что за время обтекания не приведет к заметному росту. Таким образом, мелкомасштабная неустойчивость Рихтмайе-ра-Мешкова несущественна.
Развитие неустойчивости Кельвина-Гельмгольца подробно исследовано в ряде работ, начиная с [14]. Показано, что нелинейная стадия роста возмущений вначале формирует устойчивые вихревые структуры. Вихри захватывают объемы вещества с противоположной стороны слоя и растут, сливаясь друг с другом. Толщина вихревых слоев 8 на порядок меньше длины пути, т.е. 8 — 0.1d. После развития вторичных трехмерных структур происходит резкий переход смешения. Характерное число Рейнольдса Re 8 (вычисленное по толщине слоя), при котором возникают развитая турбулентность и молекулярное перемешивание, не менее 104 [15]. Число Рейнольдса Re, соответствующее диаметру частицы, на порядок больше. Следовательно, при умеренно больших числах Рейнольдса ^е меньше или порядка 104) сдвиговая неустойчивость не успевает развить стадию турбулентного перемешивания. Но она важна как фактор, усиливающий диффузию за счет появления развитой поверхности и диспергирования. В частности, в составе тротил/гексоген скорость проскальзывания оценена В.М. Титовым и В.В. Митрофановым [16] в 200 м/с. При вязкости V = 0.04 см2/с [17] и размере зерна 0.2 мм получаем Re = 104. С ростом числа Рейнольдса до 105 и выше роль неустойчивости Кельвина-Гельмгольца возрастает.
Таким образом, представляет интерес моделирование гидродинамических эффектов в дисперсной системе двух газов. Исходя из предыдущих оценок, в данной работе основное внимание уделяется проскальзыванию, приводящему к наиболее важной неустойчивости Кельвина-Гельмгольца.
3. Ячеечная модель скольжения
Гидродинамическое взаимодействие газов исследовалось в основном для одиночного включения [1, 2], т.е. разреженной системы. Взаимное влияние близко расположенных включений наблюдалось, например, в [18] для системы неоднородностей («газовой завесы», которую можно представить как набор цилиндрических пузырей).
В интересующем нас случае объемная доля обоих газов в системе не мала. Поэтому естественно рассматривать стесненное «ячеечное» обтекание, с периодическими граничными условиями на границе ячейки — расчетной области, содержащей неоднородность. С точки зрения расчета постановка получается более простой, чем для одиночного включения, когда граничные условия нетривиальны.
Относительная скорость проскальзывания и вырабатывается при прохождении фронта ударной волны, но при рассмотрении последующих процессов разумно задать ее как начальное условие. Поскольку мы интересуемся случаем, когда плотности веществ близки, а скорость проскальзывания задается, далее в первом приближении рассмотрим обтекание газового пузыря газом той же плотности.
Таким образом, простейший случай взаимодействия включений между собой и с несущим газом можно охарактеризовать как ячеечную модель с импульсным возбуждением проскальзывания. Основная часть расчетов проведена в двумерной постановке, с использованием уравнения состояния идеального газа.
Следует отдельно обсудить применимость описанной постановки к случаю детонации гетерогенных взрывчатых веществ. Будем рассматривать движение в системе отсчета, движущейся со средней скоростью продуктов детонации 2 км/с). Скорость проскальзывания 200 м/с) мала по сравнению со скоростью звука 6 км/с), и течение глубоко дозвуковое. Поэтому на стадии развития мезоскопических течений сжимаемость веществ несущественна, и при рассмотрении обтекание неоднородностей продуктов детонации можно описывать простейшим уравнением состояния идеального газа, причем не важны и такие детали, как конкретное значение показателя адиабаты. С другой стороны, реальные свойства взрывчатых веществ и продуктов детонации сказываются на этапе формирования поля скоростей, что может повлиять и на последующее течение. Некоторые проявления этого влияния обсуждаются в
Рис. 1. Геометрия решетки и возможные состояния для метода LBE
разделе 6. В данной работе мы в основном ограничимся общими для любой газовой системы закономерностями развития мезоскопических течений, не углубляясь в специфику взрывчатых веществ (которая при таком подходе, в принципе, может быть учтена заданием начальных условий для расчетов).
4. Численные подходы
В данной работе течение вещества рассчитывалось с использованием двух различных подходов: решеточного и конечно-разностного.
Решеточные методы возникли как альтернатива традиционным расчетным газодинамическим подходам. В методе решеточного уравнения Больцмана (lattice Boltzmann equation, LBE [19]) сплошная среда рассматривается как ансамбль одинаковых частиц, которые могут двигаться по ребрам регулярной пространственной сетки из одного узла в другой, так что существует небольшое число возможных скоростей частиц c k. Ищется функция распределения частиц по возможным состояниям, а макроскопические параметры (такие как плотность и скорость течения) определяются уже из этой функции. Оказывается, что уже при небольшом числе состояний удается достаточно точно описать гидродинамику.
В данной работе для двумерных расчетов использована 9-скоростная система (рис. 1). Частицы единичной массы населяют узлы квадратной сетки и движутся по сторонам квадрата (скорость 1, энергия 1/2) либо диагоналям (скорость V2, энергия 1). Кроме того, имеются покоящиеся частицы. За шаг по времени At = 1 все движущиеся частицы переносятся в соседние узлы, где происходят столкновения, изменяющие заселенности состояний. В качестве переменных используются компоненты функции распределения Nk, k = 0, ..., 8, уравнения эволюции которых имеют вид
Nk (x + c k At, t + At) = Nk (x, t) + Q k (N (x, t)),
где c k — векторы возможных скоростей; Q k — оператор столкновений. Фактически, решается кинетическое уравнение Больцмана для модельной системы частиц. Показано, что метод решеточного уравнения Больцмана приводит к уравнениям типа Навье-Стокса [20].
В настоящее время, в основном, используется BGK-форма оператора столкновений, подобная предложенной Bhatnagar, Gross, Krook для классического уравнения Больцмана [21]. Она представляет собой релаксацию к локальному равновесному состоянию
Qk (N) = - (Nk - Ne4)/т.
Равновесные функции распределения N|4 зависят от макроскопических параметров — плотности р =
= 2 Nk и вектора скорости вещества в узле u =
k
= 2 Nk c kl P и обычно выбираются в виде, близком к
k
максвелловскому. Поскольку метод используется для малых скоростей течения, N|4 представляются в виде ряда с точностью до членов порядка и2:
..................2 ^
Nkq = pwk
^ + (ck U) + (ck U)2
2T 2
u
2T
Здесь Т — безразмерная температура. Заселенности Wk зависят только от модуля скорости данного состояния, но не от ее направления, так что имеется четыре неизвестных: w(0) = w0 (заселенность для покоящихся частиц), w(1) = ™1+4 (для частиц, движущихся по сторонам квадрата), w(л/2) = и^8 (для «диагональных» частиц) и Т. В двумерном случае необходимо выполнение следующих равенств:
Е N1* =р,
к
Е ^Чс к = ри
к
Е N^1 = р(2Т + и 2).
к
Поскольку возможны три уровня энергии частиц: 0, 1/2 и 1, добавляется еще одно «больцмановское» уравнение w(0)/w(1) = и'(Г)/^^2)- Решение дает Т = 1/3, w(0) = 4/9, w(1) = 1/9, w(V2) = 1/36 [22]. Такие модели описывают течения газа фиксированной температуры, с изотермической скоростью звука с8 = 4Т = 1/,/3 (причем должно быть и << с8).
Время релаксации т определяет коэффициенты вязкости и диффузии: V = D = (т-1/2)/3. Выбором т в пределах ^2 < т < ^ можно регулировать свойства среды. При т < ^2 расчет теряет устойчивость.
BGK-вид оператора столкновений обеспечивает галилеевскую инвариантность и легко обобщается на трехмерный случай. Нами для трехмерных расчетов ис-
пользовалась 15-скоростная модель (модули скорости 0, 1, -/3, заселенности w(0) = 2/9, w(1) = 1/9, w(V?) = = 1/72) [22].
Для отслеживания распределений концентрации веществ и диффузии в решеточных подходах вводились два сорта частиц — «красные» и «синие», вначале пространственно разделенные.
«Расщепление» плотности на несколько составляющих Nk является некоторым усложнением описания. Но это более чем компенсируется резким упрощением расчета эволюции системы. По самой природе алгоритма масса и импульс течения сохраняются автоматически. Не возникает проблем с нелинейностью уравнений гидродинамики, с деформацией счетных ячеек, как в лагранжевых методах, или с аппроксимацией потоков через границы, как в эйлеровых. Примененные в данной работе изотермические варианты LBE являются простейшими и пригодны для расчетов в дозвуковом, почти несжимаемом случае. Однако недавно развиты гибридные методы [23], в которых явно вводится конечно-разностное уравнение энергии, а динамика «управляется» решением уравнения Больцмана. Такие подходы могут применяться и для течений с переменной температурой, сжимаемостью, ударными волнами.
Для сравнения с решеточными расчетами и для прямого моделирования ситуаций, в которых решеточные методы мало пригодны, использовался комплекс газодинамических программ ЭГАК [24]. Данный комплекс предназначен для численного моделирования двумерных (плоских и осесимметричных) нестационарных течений многокомпонентной среды, характерной особенностью которых является наличие больших деформаций, включая турбулентные течения. В комплексе реализованы методики, позволяющие проводить расчеты широкого класса задач, в том числе газо- и гидродинамических течений с учетом горения и детонации взрывчатых веществ, а также турбулентного перемешивания. Для индивидуализации различных веществ используется полный набор термодинамических парам ет-ров, в случае газовой динамики это удельные внутренние энергии, плотности и объемные концентрации. Могут рассматриваться разные вещества со своими уравнениями состояния; разные фазовые состояния одного вещества, такие как взрывчатые вещества - продукты детонации, вода - пар и др.; вакуум; абсолютно твердые тела. Количество веществ как в целом по счетной области, так и в отдельных смешанных ячейках не ограничено.
В комплексе используется континуальный подход для всех веществ и непрерывное представление потоков через стороны счетных ячеек. Для предотвращения размывания контактных границ в процессе счета используется метод концентраций [25]. При расчете потоков из смешанных ячеек анализируется поле концентраций
Рис. 2. Развитие неустойчивости и перемешивание при обтекании цилиндра: расчеты методом LBE (а); прямое газодинамическое моделирование без турбулентности (б); газодинамический расчет с учетом турбулентности (в). Для всех вариантов сверху вниз показаны моменты безразмерного времени t|tf: 0, 0.5, 1, 1.5, 3, 4.5. В варианте (а) зона существенного смешения выделена черным, в вариантах (б) и (в) — светлым
и восстанавливается положение контактной границы, а затем на этой основе определяются потоки веществ. Этот метод позволяет локализовать положение контактных границ с точностью до одной счетной ячейки. В расчетах могут быть использованы и смеси (например, при турбулентном перемешивании), в этом случае восстановление контактных границ не производится и потоки определяются для смесей.
Моделирование турбулентного перемешивания с помощью комплекса ЭГАК может производиться как напрямую, без привлечения каких-либо моделей турбулентности, так и с использованием полуэмпирических моделей. В настоящей работе использовалась к-е-мо-дель [26, 27].
Аппроксимация уравнений производится в лагран-жево-эйлеровых переменных. Конечно-разностные схе-
мы, использующиеся в комплексе, строятся интегро-дифференциальным способом, имеют первый порядок аппроксимации по времени и первый или второй порядок — по пространственным переменным. Более подробное описание методики расчета можно найти в работе [28].
5. Результаты расчетов
В качестве стандартного случая выбрано двумерное обтекание первоначально круглой «капли», или пузыря, газа, движущейся влево, потоком газа, движущимся вправо (рис. 2). Объемные доли обоих веществ одинаковы, начальные скорости в каждом веществе однородны и противоположны по знаку (их =± и/2, vy = 0). Граничные условия периодические по обеим координатам, что отражает стесненность обтекания. Реально
рассчитывалась половина квадратной ячейки, а на рисунках для наглядности достроены прилегающие области с учетом симметрии и периодичности. Приведены результаты трех параллельных расчетов: по методу LBE при количестве узлов в ячейке 1600 х 1600 и числе Рейнольдса Re = ий/V ^ 12000 (рис. 2, а), по методу ЭГАК (рис. 2, б) и по методу ЭГАК с учетом турбулентности (рис. 2, в). В двух последних расчетах использовалось уравнение состояния идеального газа с показателем адиабаты у = 3, при размере расчетной сетки 200x200; материальная вязкость не учитывалась. Схемная вязкость в чисто сдвиговом течении согласно теоретическим оценкам также должна быть равна нулю. Действительно, в специальном расчете сдвигового течения (первой задачи Стокса) по методу ЭГАК вязкость зафиксировать не удалось. Для решеточного же расчета вязкость, оцененная в тестовой задаче, хорошо совпадала с теоретическим значением. Картины течения даны для одинаковых моментов безразмерного времени , где
= й/и — характерное время обтекания.
Начальные стадии расчетов хорошо согласуются. На момент 1.5/у- крупномасштабные деформации пузырей практически одинаковы. Граница сильно растягивается и искривляется с образованием вихрей. Газодинамический расчет дает почти прямые углы на наветренной стороне пузыря и резкий «пик» на оси симметрии. Эти детали размазаны в решеточном варианте вязкостью и диффузией. Длины волн неустойчивостей поверхности и амплитуды волн в решеточном расчете приблизительно те же, что и в газодинамическом, хотя мелкомасштабные детали различаются, что связано с разными числами Рейнольдса (по существу, с различием в вязкости среды).
На больших временах газодинамический расчет дает заметное диспергирование. В областях, затронутых диспергированием, проскальзывание быстро прекращается. В результате меняются и условия обтекания «оставшейся» капли. На момент 4.5 /у- можно говорить лишь
о качественном согласии расчетов по форме включений.
Газодинамические расчеты с учетом турбулентности дали результаты, очень близкие к расчетам без турбулентности, с той разницей, что мелкомасштабные детали размывались (рис. 2, в).
Таким образом, различие в результатах решеточного и классического расчета связано в основном с разными свойствами моделируемых сред. Можно также считать, что описанный решеточный расчет лучше приспособлен к описанию поведения малых включений, а газодинамический — для сравнительно больших. При увеличении числа узлов сетки (чему препятствовали только технические причины) следует ожидать дальнейшего сближения результатов.
Представляет интерес сравнение расчетов стесненного обтекания со случаем уединенных включений. В
Рис. 3. Влияние числа Рейнольдса на обтекание
экспериментах [2] демонстрируется обтекание сравнительно тяжелых «газовых цилиндров» за слабой ударной волной. Для небольших времен результаты сходные, с естественным отличием, что более плотные цилиндры из SF6 деформируются медленнее: стадии наших расчетов t|tf = 0.5, когда форма включения близка к полукругу, соответствует в [2] заметно большее время = 1.5. Однако далее деформации качественно отличаются. В [2] из полукруга развивается сначала тонкий серп, значительно превосходящий исходный цилиндр по поперечным размерам, а затем с острых концов серпа стекают четко выраженные вихри, в окрестности которых локализуется неперемешанный материал включения. В наших же расчетах деформация в основном продольная: из полукруга развивается «наконечник стрелы» (/Ду = 1.5), а затем — «человеческая фигура» (/Ду = 3) и «ракета» (/Ду = 4.5). Это различие связано с блокированием поперечного расширения при стесненном обтекании.
В эксперименте [2] не наблюдалось проявлений сдвиговой неустойчивости. Это объясняется размытием границы раздела цилиндрической струи и окружающего воздуха, хорошо видным на иллюстрациях (исходная толщина диффузионного слоя была около половины радиуса цилиндра). Поэтому коротковолновые возмущения подавлялись, а длинноволновые, если и развивались, то медленнее основной моды деформации включения. Отметим здесь, что своеобразие неоднородных продуктов детонации состоит, в том числе, и в осуществ-
лении достаточно резких границ раздела. В более «мягких» атмосферных условиях относительно столь же резкие границы требуют существенного увеличения масштаба. В опытах [1] наблюдалась хорошо выраженная неустойчивость границ пузырей. Здесь резкость границы достигалась использованием разделяющей газы разрушаемой пленки, которая, однако, могла вносить искажения в процесс взаимодействия веществ.
При обтекании более мелких включений зона диффузии относительно расширяется. Это приводит к подавлению волн на границе раздела. На рис. 3 показаны несколько конфигураций, полученных методом LBE при одном и том же безразмерном времени г/г^ = 1.5, но при различных числах Рейнольдса: 12000, 5 600, 2500 и 1600 (сверху вниз). При уменьшении Re до 5600 неустойчивость Кельвина-Гельмгольца на границе уже проявляется слабо, а при еще меньших масштабах возмущения не развиваются; виден диффузионный слой (область, в которой разница концентраций веществ менее 30 % от начальной), выделенный черным цветом. При числах Re порядка сотен диффузионный слой становится довольно толстым, что приводит даже к перекрытию языков перемешанного вещества. При этом основная деформация включения все еще хорошо воспроизводится.
Поэтому для продвижения в диапазон больших размеров допустимо отказаться от полного моделирования по вязкости, ограничиваясь условием Re >> 1. Например, условия на рис. 2, а соответствуют размеру пузыря d = 10 мм, скорости обтекания 18 м/с (/у = 0.55 мс) и вязкости V = 0.15 см2/с. Но можно принять и d = 100 мм, г^ = 5.5 мс. При сохранении скорости обтекания и числа Рейнольдса вязкость будет преувеличена в 10 раз, но вязкие силы останутся малыми. Это аналогично искусственной или схемной вязкости в традиционных газодинамических расчетах. Если расчет претендует только на воспроизведение основных мод искажения поверхности, такое частичное моделирование имеет смысл; в ряде случаев оно используется ниже. Меньшие числа Рейнольдса соответствуют меньшему количеству узлов решетки и, значит, меньшему объему вычислений.
В качестве интегральной меры степени взаимодействия веществ использовалась эффективная электрическая проводимость смеси. Эта характеристика определяется искажениями начальной геометрии и перемешиванием и может быть сопоставлена с экспериментальными результатами [7]. Считалось, что проводимость вещества «капель» несущественна по сравнению с проводимостью несущего газа, как это имеет место для пары тротил/гексоген. Проводимость вычислялась по полученным гидродинамическим конфигурациям. Задавались значения электрического потенциала: ф = 0 на одной из граней ячейки, ф = ф 0 — на противоположной, и для каждого ребра сетки, соединяющего узлы I и
1 1 1 1 1 > Ж' 1 2 3 4 5 6
;-7,—'Ч1
1111
0 1 2 3 4 ^
Рис. 4. Зависимость электрической проводимости от безразмерного времени при обтекании круглых капель: 1-5 — расчет методом LBE, Яе - 400 (1), 800 (2), 1 600 (3), 5600 (4), 12000 (5); 6 — расчет с
использованием комплекса ЭГАК
назначалось значение проводимости Gij = а0hлJninj . Здесь а0 — постоянная величина; п{ , nj■ — «эффективные концентрации» проводящей фазы в конечных узлах ребра (п = тах(р1 - р2, 0)); h — шаг сетки. Формула для Gij — модельная, но она отражает как «разбавление» хорошо проводящих продуктов детонации тринитротолуола (уменьшение р1), так и «отравление» проводимости при смешении (что может быть результатом вторичных химических реакций между продуктами детонации тротила и гексогена, приводящих, например, к выжиганию свободного углерода). Методом установления рассчитывалось распределение потенциала. После этого находился ток через ячейку I = ^ Дф i,
i
суммирование производилось по слою, прилегающему к одной из граней с постоянным потенциалом. Для квадратной ячейки проводимость G = I/ ф 0. Проводимости ячейки в горизонтальном направлении Gx и в вертикальном Gy, как правило, не одинаковы. Далее на графиках показано среднее значение (р) = (Gx + Gy )/2. На рис. 4 представлена зависимость проводимости от безразмерного времени г/г^ для различных чисел Рейнольдса. Графики близки, т.е. взаимодействие практически автомодельно. Характерное время спада проводимости га - (1 - 2) гу. Поскольку ~ d, время т растет
с увеличением размера включения, что и наблюдалось в [7]. Для зерен гексогена размером 200 мкм га ~ 12 мкс, что не противоречит экспериментам [7]. Для миллиметровых частиц га ~ 5-10 мкс — заметно больше экспериментального. Однако необходимо учесть, что перемешивание — не единственный процесс, приводящий к спаду электропроводности: из-за общего расширения спад наблюдался и в чистом тротиле. С учетом этого замечания можно говорить об удовлетворительном объяснении экспериментальных данных.
Представляет интерес влияние начальных особенностей формы включений на течение. На рис. 5 показаны различные стадии развития неустойчивости при обтекании первоначально квадратных «капель», при раз-
Рис. 5. Развитие обтекания квадратных включений: «прямой» квадрат (а); наклонный (б). Расчет с использованием комплекса ЭГАК без учета турбулентности. Время для обоих вариантов (сверху вниз): г/гу: 0, 1.5, 3
личных ориентациях относительно потока (газодинамическое моделирование). Процесс в целом сходен с рис. 3, но перемешивание более быстрое, по-видимому, из-за наличия углов, на которых неустойчивости развиваются быстрее. Об этом же говорит рис. 6, проводимость спадает быстрее, чем для круглых капель, хотя качественно графики близки.
Интересно исследовать, как меняется картина перемешивания при изменении начального расположения включений. В следующей серии расчетов методом LBE начальное расположение капель было диагональным (квадратная решетка, повернутая на угол 45°), т.е. координаты центра одной из круглых капель (0,0) (из-за периодических граничных условий эта капля выглядит как четыре области в углах расчетной ячейки). Другая капля находилась в центре ячейки. Проведены также расчеты, в которых одна подсистема цилиндров смещена относительно другой в горизонтальном либо вертикальном направлении (это достигалось смеще-
1 —
2
з
- \
1 1 I 4
0 1 2 3 4 ^
Рис. 6. Зависимость проводимости от безразмерного времени для капель различной формы (газодинамический расчет): цилиндр (1); квадрат (2); наклонный квадрат (3)
нием центрального цилиндра на половину радиуса (рис. 7, б, в) или на четверть радиуса по обоим направлениям (рис. 7, г). В этих расчетах размер ячейки L = = 400, скорость обтекания и = 0.2, кинематическая вязкость V = 1/12, диаметр капель d = - 225, харак-
терное время обтекания гу - 560, число Рейнольдса Яе - 1 080. Развитие неустойчивостей при различных начальных конфигурациях представлено на рис. 7. Для всех вариантов верхний кадр соответствует г/гу = 0, второй — г/г^ = 1, третий — г/г^ = 3, четвертый —
фу =5
Картины при «шахматном» начальном расположении заметно отличаются от «квадратного» случая. Течение довольно быстро перестраивается от в основном горизонтального к в основном вертикальному с системой вихрей. Происходит значительная деформация капель с фрагментацией и объединением на поздних стадиях. При этом в некоторых случаях происходит перекрытие ячейки в одном из направлений (например, рис. 7, б-г). На рис. 8 представлена зависимость проводимости ячейки от времени для шахматного расположения. Спад проводимости для «сдвинутых» конфигураций более быстрый.
6. Более общие постановки
Все предыдущие расчеты были выполнены в двумерной постановке. В реальном трехмерном течении развитие неустойчивостей происходит более интенсивно (неоднородность по третьей оси можно представить себе как дополнительное возмущение). С другой стороны, связность каналов в трех измерениях менее подвержена нарушению, чем в двумерном случае. Указанные факторы противоположным образом влияют на эффективную электрическую проводимость среды.
Рис. 7. Течение при обтекании «диагональной» системы пузырей. Средний цилиндр: в центре ячейки (а); сдвинут влево (б); сдвинут вниз (в); сдвинут вниз и влево (г)
Было проведено трехмерное моделирование в импульсной ячеечной постановке методом LBE. На рис. 9 показаны «трехмерные» графики проводимости для
Рис. 8. Зависимость электрической проводимости от безразмерного времени для шахматного расположения капель. Кривая 1 соответствует рис 7, а; кривая 2 — рис 7, б; кривая 3 — рис 7, в; кривая 4 — рис 7, г
трех чисел Рейнольдса: 500, 1 000 и 1400. Видно, что существенного отличия от двумерных данных нет, т.е. трехмерность, по крайней мере интегрально, не сказывается на результатах.
Рис. 9. Зависимость электрической проводимости от безразмерного времени для трехмерного расчета: Яе - 500 (1); 1000 (2); 1400 (3)
а
0 5 10 15 20 25 30 35 40 45 50
1 = 6
15 20 25 30 35 Т 45 50
25 30 35 40 45 50 55 60 65 70
Рис. 10. Развитие проскальзывания в цепочке, содержащей 20 пузырей. Время указано в микросекундах, координата — в миллиметрах. Для ґ = 6 мкс стрелкой показано положение фронта ударной волны. Последний кадр (ґ = 18 мкс) сдвинут влево относительно двух первых
Формирование проскальзывания моделировалось методом ЭГАК в двумерной постановке. По покоящейся неоднородной среде пропускалась ударная или детонационная волна. Расчетная область представляла собой «полоску» с 20 круглыми включениями радиуса 1 мм, по существу — 20 последовательных ячеек. Начальные плотности в пузырях и несущей среде относились как Р2/Р1 = 1.82/1.63 (А = 0.055), уравнения состояния — идеального газа с у = 3 (включения) и у = 3.2 (несущий газ), что соответствует свойствам гексогена и тротила. Начальное давление однородно. Левой границе сообщалось движение со скоростью 1.25 мм/мкс, приводящее к сжатию «усредненного» вещества в 1.23 раза (число Маха М = 1.28). На рис. 10 видно, как постепенно развивается течение за ударной волной. Если не брать самые крайние пузыри, на которые влияют граничные условия, то деформации включений практически те же, которые получены в импульсной постановке. Пузыри в ударной волне на ранних стадиях сплющиваются, но это не сказывается сколько-нибудь заметно на их обтекании для больших интервалов времени. Не видно и проявлений неустойчивостей Рэлея-Тейлора и Рихт-майера-Мешкова. Полоска дает пространственную развертку временного хода течения в отдельной ячейке. Характерное время обтекания примерно соответствует последней стадии расчета. Отметим, что такой расчет
мало реален для метода LBE, более приспособленного для расчета почти несжимаемого течения.
Наконец, интерес представляет случай гетерогенной смеси двух взрывчатых веществ, когда за фронтом ударной волны выделяется энергия. В той же полоске, при тех же начальных условиях задавались теплоты реакции Q: 5.4167 кДж/г (гексоген) и 4.1667 кДж/г (тротил). Реакция считалась мгновенной, что соответствует достаточно большим размерам включений. Граничные условия — жесткие стенки, волна идет слева направо. На рис. 11 показаны распределения концентраций, давления и плотности незадолго до момента выхода фронта волны из расчетной области. Давление, как и следовало ожидать, выравнивается быстро (несколько периодов решетки), а концентрация на этих временах существенно не «размешивается». Несколько медленнее, чем давление, выравнивается плотность (из-за большего энерговыделения пузыри, вначале более плотные, расширяются). Формы пузырей несколько другие, чем в инертной среде, что также связано с энерговыделением: наряду с деформацией в потоке играют роль пульсации включений.
Представляет некоторый интерес более полное моделирование процессов в гетерогенном взрывчатом веществе, начиная с исходного состояния веществ. При этом учет реального уравнения состояния не представ-
б
Рис. 11. Детонация в газе с неоднородностями. Размеры указаны в мм, время с момента инициирования 4 мкс; объемные концентрации (а); растровая картина давления (б); растровая картина плотности (в)
ляет проблемы. Но, к сожалению, на данный момент трудно учесть различия в кинетике реакций. Применяемые сейчас кинетики — макроскопические, и их пригодность к описанию процессов на уровне зерен взрывчатых веществ далеко не очевидна. Поэтому принятый в данной работе простейший «импульсный» подход дает, возможно, оптимальное (хотя и качественное) описание реального сложного процесса.
7. Заключение
Моделирование течения газовой среды с включениями другого газа демонстрирует существенное гидродинамическое взаимодействие включений и несущего газа, приводящее к сильной деформации исходной структуры. Характерный временной масштаб задается временем обтекания включения. В диспергировании и перемешивании преобладающую роль играет сдвиговая неустойчивость Кельвина-Гельмгольца. В концентрированной системе стесненность обтекания порождает формы, качественно отличные от тех, которые получаются из одиночных неоднородностей. Всякая начальная неидеальность (негладкость границ, несимметрия расположения включений) ускоряет взаимодействие.
Двумерная ячеечная постановка с импульсным возбуждением скольжения оказалась вполне адекватной, несмотря на ее простоту. Заметные отличия получены только для реагирующей системы, хотя качественные закономерности сохраняются. Более серьезное отличие возможно при существенно разных временах реакции веществ, что может вызвать пульсационные движения с заметным ускорением перемешивания.
Решеточные расчеты показали хорошее совпадение с результатами апробированных газодинамических кодов. Для обоих подходов существуют «ниши», где их использование предпочтительно. Сравнение результатов способствовало лучшему пониманию процесса.
В применении к гетерогенным взрывчатым веществам расчеты согласуются с известными изотопными данными [3, 4] и измерениями электрической проводимости [7]. При микронных размерах неоднородностей за субмикросекундные времена существенно диффузионное смешение. При миллиметровых размерах зерен перемешивание слабое, а для штатной дисперсности гексогена около 200 мкм степень перемешивания в течение микросекунд можно оценить как заметную.
Работа выполнена при поддержке РФФИ (гранты №№ 99-03-32336, 02-03-32873).
Литература
1. Haas J.-F., Sturtevant B. Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities // Journal of Fluid Mechanics. - 1987. - V. 181. - No. 1. - P. 41-76.
2. Jacobs J.W. The dynamics of shock accelerated light and heavy gas cylinders // Physics of Fluids A. - 1993. - V. 5. - No. 9. - P. 22392247.
3. Aнuсuчкuн B-Ф., Дерендяев Б.Г., Kqптюг B.A., Мальков И.Ю., Ca-лахутдинов H-Ф., Титов B-М. Исследование процесса разложения в детонационной волне изотопным методом // Физика горения и взрыва. -1988. - Т. l4. - М 3. - С. 1l1-1ll.
4. Koзырев H.B., Брыляков П.М., Cy Ceн Чел, Штейн МЛ. Исследование процесса синтеза ультрадисперсных алмазов методом меченых атомов // Доклады AR СССР. - 199Q. - Т. 314. - М 4. - С. 889891.
5. Лямкин A-И., ПетровE.A., Ершов A-П., CaкoвuчLB., Cmaвeр A-М.,
Титов B-М. Получение алмазов из взрывчатых веществ // Доклады An СССР. - 1988. - Т. 3Ql. - М 3. - С. 611-613.
6. Титов B-М., Aнuсuчкuн B-Ф., Мальков И.Ю. Исследование процесса синтеза ультрадисперсного алмаза в детонационных волнах // Физика горения и взрыва. - 1989. - Т. l5. - М 3. - С. 117-1l6.
7. Ершов A-П., Camoнкuнa H-П., Дибиров O.A., Цыкин C.B., Янилкин Ю-B. Исследование взаимодействия компонентов гетерогенных взрывчатых веществ методом электропроводности // Физика горения и взрыва. - lQQQ. - Т. 36. - М 5. - С. 97-1Q8.
8. Медведев Д-A., Ершов A-П., Янилкин Ю-B. Моделирование переме-
шивания в двухкомпонентной системе // Экстремальные состояния вещества. Детонация. Ударные волны. Труды Международной конференции. III Харитоновские тематические научные чтения, Саров, lQQ1. - Саров: РФЯЦ - ВНИИЭФ, lQQl. - С. 247-253.
9. Медведев Д-A., Ершов A-П., Kyперштох A-Л. Численное исследование гидродинамических и электрогидродинамических неустойчивостей // Математические проблемы механики сплошных сред. Матер. 5-й Сибирской школы-семинара. - Динамика сплошной среды. - Вып. 1lQ. - Новосибирск: Ин-т гидродинамики, lQQl. -С. 93-1Q3.
1Q. Youngs D.L. Modelling turbulent mixing by Rayleigh-Taylor instability // Physica D. - 1989. - У 37. - No. 1-3. - P. l7Q-l87.
11. Dimonte G., Schneider M. Turbulent Rayleigh-Taylor instability experiments with variable accelerations // Physical Review E. - 1996. -У. 54. - No. 4. - P. 374Q-3743.
1l. Янилкин Ю-B., Cmaцeнкo B.П., Ребров C.B., Cuнькова О.Г., Cmaд-ник A-Л. Исследование гравитационного турбулентного перемешивания при больших разноплотностях с помощью прямого трехмерного численного моделирования // ВAНТ. Сер. ММФП. -lQQl. - Вып. l. - С. 3-9.
13. Rudinger G., Somers L.M. Behavior of small regions of different gases carried in accelerated gas flows // Journal of Fluid Mechanics. -196Q. - У. 7. - No. l. - P. 161-176.
14. Brown G.R., Roshko A. On density effects and large structure in turbulent mixing layers // Journal of Fluid Mechanics. - 1974. - У. 64. -No. 4. - P. 775-816.
15. SlessorM.D., Bond C.L., DimotakisP.E. Turbulent shear-layer mixing at high Reynolds numbers: effects of inflow conditions // Journal of Fluid Mechanics. - 1998. - У. 376. - P. 115-138.
16. Титов B-М., Митрофанов B.B., Ершов A-П., Kyперштох A-Л., Мальков И.Ю. Углерод в детонационных процессах (часть Б). Отчет, выполненный для Ливерморской лаб. - Новосибирск: Институт гидродинамики СО РAН, 1994. - 69 с.
17. Давыдова O.H., Kyзнецов H-М., Лавров B.B., Шведов K.K. О недосжатой детонации конденсированных взрывчатых веществ с инертными примесями // Химическая физика. - 1999. - Т. 18. -М 4. - С. 53-66.
18. Rightley PM., Vorobieff P., Martin R., Benjamin R.F. Experimental observations of the mixing transitions in a shock-accelerated gas curtain // Physics of Fluids. - 1999. - У. 11. - No. 1. - P. 186-lQQ.
19. McNamara G., Zanetti G. Use of the Boltzmann equation to simulate lattice-gas automata // Physical Review Letters. - 1988. - У. 61. -No. lQ. - P. 2332-2335.
lQ. Chen H., Chen S., Matthaeus W.H. Recovery of the Navier-Stokes equations using a lattice Boltzmann method // Physical Review A. -199l. - У. 45. - No. 8. - P. R5339-R534l.
l1. BhatnagarP., GrossE.P., KrookM.K. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-
component systems // Physical Review. - 1954. - V. 94. - No. 3. -P. 511-525.
22. Qian Y. H., d’Humieres D., Lallemand P. Lattice BGK models for Navier-Stokes equation // Europhysics Letters. - 1992. - V. 17. -No. 6. - P. 479-484.
23. ZhangR., ChenH. Lattice Boltzmann method for simulation of liquid-vapor thermal flows // Physical Review E. - 2003. - V. 67. - No. 6. -P. 066711-1-066711-6.
24. Янилкин Ю.В., Шанин А.А., Ковалев Н.П., Гаврилова Е.С., Губ-ков Е.В., Дарова Н.С., Дибиров О.А., Жарова Г.В., Калмано-вичА.И., Павлуша И.Н., Самигулин М.С., Симонов Г.П., Синько-ва О.Г., Сотникова М.Г., Тарасов В.И., Торопова Т.А. Комплекс программ ЭГАК для расчетов двумерных течений многокомпонентной среды // ВАНТ. Сер. ММФП. - 1993. - Вып. 4. - С. 69-75.
25. Бахрах С.М., Глаголева Ю.П., Самигулин М.С., Фролов В.Д., Яненко Н.Н., Янилкин Ю.В. Расчет газодинамических течений на
основе метода концентраций // Доклады АН СССР. - 1981. -Т. 257. - № 3. - С. 566-569.
26. Янилкин Ю.В., Никифоров В.В., Жарова Г.В. Модель с двумя уравнениями и методика расчета турбулентного перемешивания в двумерных сжимаемых течениях // ВАНТ. Сер. ММФП. - 1994. -Вып. 4. - С. 79-84.
27. Yanilkin Yu.V, Nikiforov V.V., Bondarenko Yu.A., GubkovE.V., Zharova G.V, Statsenko VP, Tarasov VI. Two-parameter model and method for computations of turbulent mixing in 2D compressible flows // Proc. 5 Int. Workshop on the Physics of Compressible Turbulent Mixing. Stony Brook, USA, 1995. - P. 105-110.
28. Янилкин Ю.В. Численное моделирование двухмерных течений многокомпонентной среды с учетом некоторых мелкомасштабных процессов // Физическая мезомеханика. - 1999. - Т. 2. - № 5. -С. 27-48.
Mesoscopic flows in a heterogeneous gas
D.A. Medvedev, A.P. Ershov, Yu.V. Yanilkin1, and E.S. Gavrilova1
M.A. Lavrentyev Institute of Hydrodynamics SB RAS, Novosibirsk, 630090, Russia 1 RFNC All-Russian Scientific Research Institute of Experimental Physics, Sarov, 607190, Russia
Consideration is given to the dynamics of a gaseous medium with the distributed inclusions of another gas in it. The system being under the impulse action, for example, during shock wave passage, one can observe the appearance of mesoscopic flows that give rise to shape distortion of the inclusions and mixing. The results are applied to analyze the detonation of heterogeneous explosives.
