Мезопроцессы и мезоструктуры в гетерогенной детонации

Исследование детонационных процессов вначале естественно концентрировалось на простейших типах. Долгое время единственным важным структурным элементом детонационной волны считалась зона химической реакции в модели Зельдовича-Неймана-Деринга [1]. Но постепенно стало возникать понимание существенной роли более сложных структур. Спиновая и многофронтовая газовая детонация — наиболее известный пример самопроизвольного зарождения сложной структуры в весьма однородной среде из общих начальных условий. Выживание именно данной структуры связано с взаимодействием гидродинамических процессов и кинетики. Последняя задает характерное время химической реакции, которое и определяет замечательно регулярную геометрию фронта и последующего течения [2].

В конденсированных ВВ стабильные режимы более распространены, и им уделяется основное внимание из-за их практической значимоститем не менее существует интерес и к пульсирующим процессам в однородных [3] и неоднородных ВВ [4, 5]. Постепенно накапливается экспериментальный материал, свидетельствующий о важности начальной гетерогенности прессованных ВВ (в работе [6], например, наличие разрешаемого современными средствами химпика связывается с гетерогенностью вещества, обеспечивающей возникновение как горячих, так и холодных точек).

Парадоксально, что для двухфазной и пористой среды, где естественно ожидать сильно неоднородных проявлений, также преобладают представления об устойчивом, плоском фронте волны и практически одномерном течении за фронтом. Это можно объяснить тем, что в таких средах затруднено распространение газодинамических возмущений и, следовательно, сколько-нибудь регулярных структур. Однако ничему не противоречитСистема (1.1) выглядит как простое обобщение баротропной газодинамической модели (а при стремлении, а к нулю газодинамика получается как предельный случай). Однако между ними существует качественное различие. Практически сразу после формулировки модели выяснилось, что двухскоростная система негиперболична.

1.1.2 Коротковолновая неустойчивостьПри решении уравнений течения комплексные характеристики означают экспоненциальный рост возмущений. Действительно, на языке дисперсионного уравнения наклон характеристики соответствует скорости волны: и/к А. Инкремент неустойчивостит.е. скорость роста обратно пропорциональна длине волны возмущения. Хотя инкремент содержит малый параметр /3, скорость роста для коротких волн (больших к) формально не ограничена. Характерное время роста возмущения с волновым числом кБыстрый рост коротких волн приводит к некорректности задачи Коши. Физическая причина неустойчивости [10, 12] - повышение скорости фильтрации газа и соответственное падение давления в области, где возникла повышенная концентрация частиц. Появившийся градиент давления «загоняет» частицы в область сгущения, и возмущение экспоненциально растет. В первоначально плавном течении растущая волна движется вместе с частицами, которые из-за взаимодействия с протекающим газом собираются в плотные слои, или «блины».

Принципиален вопрос о природе этой неустойчивости: является она реальным физическим эффектом или формальной неустойчивостью, возникающей из-за неточностей описания двухфазной среды. Именно такой неточностью может объясняться сколь угодно быстрый рост коротковолновых возмущений. В пользу нефизичности этой картины говорит хотя бы возможность (и даже преимущественная) образования «блинов» толщиной меньше размера частиц. Выяснению этой ситуации и посвящен данный раздел. Im и = Ыт Л k/3w,(1.4)Обзор вариантов регуляризацииВ [12] отмечено различие в постановках задач, допускающих постоянное проскальзывание фаз (взвешенные слои, осаждение суспензий) и существенно нестационарных, типа прохождения ударной волны через газовзвесь.

Для первого случая вообще проблематична возможность «одномерного» поведения. Взаимодействие частиц через поток действительно приводит к образованию агрегатов. В псевдоожиженном слое типично как раз не однородное состояние среды, а чередование газовых «пузырей» и плотных «комков». В этом плане рассмотренная выше неустойчивость — это как бы одномерная карикатура на более сложное и богатое поведение реальной системы. Другая сторона взаимодействия частиц — развитие хаотических движений в дисперсной фазе, что также выходит за рамки простейшей двухскоростной модели. Хаотические движения сообщают псевдогазу частиц определенную упругость, и можно ожидать некоторого подавления концентрационной неустойчивости (во всяком случае, коротковолновой). Этому направлению посвящены работы [12, 13].

Напротив, если скорость относительного движения фаз стремится к нулю со временем (условно можно говорить о задачах с ограниченным временем проскальзывания) модель типа (1.1) может быть адекватной. Однако неограниченность роста коротких волн должна быть устранена. При численном решении таких задач это всегда подразумевается. Неустойчивость может гаситься либо неявно самой вычислительной процедурой (схемной вязкостью), либо специальными мерами (сглаживанием коротковолновых возмущений). Согласно [14], рост возмущений сглаженной (усредненной на размере, малом по сравнению с характерным размером течения) величины, а ограничен, и в такой норме задача Коши для системы (1.1) корректна. В [12] предложено ввести искусственное отталкивание частиц, физически понимаемое как результат хаотических движений в дисперсной фазе.

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

Оценка роста возмущенийОценим вначале характерное время выравнивания скоростей фаз. Пусть в начальный момент скорости фаз различны. Рассмотрим простейший случай однородного по пространству состояния, когда все величины не зависят от х. Тогда из (1.1) выпадают члены с градиентами, и остается чистое затухание относительной скорости: В практически важном случае закон сопротивления квадратичный: где, А — поперечное сечение частицы, Cd коэффициент сопротивления. Для «не слишком обтекаемых» частиц Cd 1. Отсюда получаем оценку времени проскальзыванияЗдесь R — размер частиц (для шаров — радиус).

Начальное возмущение заметно возрастет при tp < te. Сравнивая (1.4) и (1.5), находим, что при этом должно бытьМинимальное значение выражения в правой части равно 2, а при значительном контрасте массовых концентраций фаз оно гораздо больше. Следовательно, до опасного уровня развиваются возмущения с длиной волны, много меньшей R.(1.5)Коротковолновый случай требует специального анализа. Пригодность простейших осредненных уравнений (1.1) здесь сомнительна. Рассмотрим основные физические эффекты, которые могли бы «придать» модели необходимую устойчивость. Ограничимся случаем разреженной смеси (a 1).

Роль хаотических движенийХаотические движения частиц приводят к возникновению «давления» псевдогаза частиц и диффузионным эффектам [12, 13].

При обтекании частиц газом на каждую частицу со стороны окружающих действуют силы гидродинамического взаимодействия через поток. При хаотическом расположении равнодействующая будет порядка силы взаимодействия соседних частиц. По результатам [15] сила взаимодействия сфер радиуса R, удаленных на расстояние I R, в идеальной жидкости // 27гpw2R6/l4. Это выражение примем как оценку реальной силы. Считая, что // действует в течение времени выравнивания скоростей фаз te, оценим набранную хаотическую скорость: Согласно [12], давление Р2 обеспечивает устойчивость, если при малых, а оно растет как ат, причем т < 2. В данной же ситуации показатель степени при, а равен 8/3 > 2. Это значит, что псевдогаз частиц не имеет достаточной упругости.

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

Давление псевдогаза частиц Р2 Psvc Ps™2**8^3lc «vcte «{А+£fY2 < Rav^ «R.PS PS Vexp (R/lc)).

Роль флуктуаций теченияДвижение несущей фазы изменяется в окрестности обтекаемой частицы.

Исследуем влияние этих флуктуаций на устойчивость модели. Для этого потребуется воспроизвести классическую процедуру вывода двухфазных уравнений [11]. Рассмотрим контрольный объем толщины Ах, изображенный на рис. 1.1. Частицы занимают долю, а объема (заштрихованную), газ — долю (р = 1 — а. Масса твердой фазы на единицу сечения равна psaАх, масса газа — pip Ах. Здесь плотность р — средняя по объему. Отметим, что в макроскопически однородной сисРис. 1.1.теме средние по контрольной поверхности значения, а и ip совпадают со среднеобъемными, что следует из чисто геометрических соображений.

Через границы объема втекает (вытекает) поток массы газа pipuAx. Очевидно, что здесь плотность и скорость газа понимаются как средние по соответствующим поверхностям. Для замыкания уравнений существенно совпадение и этих средних со среднеобъемными величинами. С учетом этих замечаний баланс массы газа записывается в видеАхdpipAx dt= (pipu (x) — pipu (x 4- Ах)) что приводит к первому уравнению системы (1.1). Уравнения массы и импульса частиц получаются подобным же образом и также имеют классический вид.

Подробнее остановимся на уравнении импульса газа, в котором будем учитывать флуктуации. Удобнее рассмотреть вначале импульс смеси, чтобы не рассматривать межфазную силу внутри контрольного объема: др<�ри др (р (и2) fdav dav2 дdt+дх+ Ps+dt ' дх J ' дх+ — (pip — та) = 0.

Здесь {и2) — средний квадрат скорости газа. Среднее напряжение в твердой фазе обозначено г (для твердых частиц нельзя ввести микроскопическихнапряжений, но средние по сечениям каждой частицы вполне определены воздействием газа).

Введем флуктуацию скорости газа 5и2 = (и2) — и2 и выделим из напряжения в твердой фазе среднее по газу давление: т = — р — т'. Знак поправки т' выбран для удобства (сжатие означает положительное т'). С учетом предыдущих уравнений получаем (ди ди др д о / В результате преобразования получено уравнение импульса газа, в которое вошли флуктуационные поправки (последний член в левой части).

Флуктуацию 5и2, возникающую при обтекании частиц, для разреженной смеси и существенно дозвуковых скоростей можно записать в виде [11]5и2 = 6(a)w2 При не слишком малых скоростях дополнительные напряжения в частицахI 1 т порядка pw — характерного перепада давления в газе вокруг частицы. Таким образом, имеемр<�р5и2 + т’а = 6(а)р (и — v)2, 5(a) а.

Существенно, что 5(a) > 0, так как в среднем т' > 0 (частица «сдавлена» ускоряющим потоком).

1.1.3 Нелокальная модельПри исследовании коротковолновой устойчивости необходимы уравнения, пригодные для быстро меняющихся течений. Между тем обычно рассматриваются плавные течения, когда характерный размер L велик по сравнению с межчастичным расстоянием I [11]. Это позволяет выбрать представительный микрообъем с размером Ах I (как на рис. 1.1) и провести по нему осреднение уравнений. Возможность осредненного подхода к течениям с быстрыми изменениями параметров может показаться проблематичной.

Искомая система уравнений может быть получена, если последовательно применять менее распространенный способ осреднения — по площади. Как уже упоминалось, потоки через границы контрольного объема — среднеповерхностные величины, и для замыкания уравнений существенно используется эквивалентность объемных и поверхностных средних. Представительная площадка, размеры которой значительно превышают I, пересекает большое количество частиц, и по ней возможно статистически достоверное осреднение. Вид общих уравнений сохранения не изменяется при переходе к среднеповерхностным величинам [16].

Однако область применимости таких величин шире, чем среднеобъем-ных. В одномерном, наиболее наглядном, случае осреднение по площадкам с нормалью, направленной вдоль движения, позволяет рассматривать течеUо.оРис. 1.2.ния с резкими градиентами. Быстрые изменения величин поперек площадки не отразятся на процедуре осреднения в ее пределах. Таким образом, для вывода уравнений может быть применен «тонкий» контрольный объем (рис. 1.2). Толщина Ах меньше I и, возможно, меньше R. Частицы, имеющие общую часть с контрольным объемом, могут дважды пересекаться его границами. Процедура вывода не отличается от общепринятой. Приведем уравнения с необходимыми комментариями.

Для твердой фазы, поскольку пренебрегается хаотическими движениями частиц, удобно рассматривать их центры масс как материальные точки. Для концентрации центров п и скорости v имеемдп dnv п rrfdnv dnv2. + = 0, psV (— + —)=nfT (1.9)где V — объем частицыfr — полная сила, с которой газ действует на частицу, центр которой в момент t находится в точке х.

Величины, относящиеся к газу, осредняются по сечениям газовой фазы, причем скорость и осредняется с весом р. Как и ранее, а и (р — доли сечения двухфазной среды плоскостью х = const. Как ясно из рис. 1.3, частица, центр которой лежит на расстоянии? от контрольной плоскости, пересекается по площади определяемой формой частицы. Тогда, а можно записать как интеграл по поперечнику частицы: R х «a (x, t) = J А (?)п (х + рис 13-Rили более коротко в виде интеграла по объему частицы V: a = JndV, <р= 1-а. (1.10)vЕсли п — медленно меняющаяся функция, (1.10) сводится к, а = riV.

Различие полученной системы (1.14) и традиционной (1.1) состоит не в виде уравнений, а в задании взаимодействий фаз. Вместо простых алгебраических соотношений (например, а = nV) фигурируют интегральные выражения, как (1.10). Для плавных течений, когда подынтегральныефункции меняются медленно, интегрирование сводится к умножению на объем, и восстанавливаются классические зависимости. В общем же случае полученная модель является нелокальной. Например, архимедова сила в данной точке х определяется давлением в интервале (х — R, х + R), поскольку действует на частицу конечного размера.

Анализ устойчивостиЗамыкание модели требует еще задания силы трения /. Для плавных течений в пренебрежении силой присоединенных масс / обычно считают функцией средних величин.

Нет оснований распространять такую зависимость на быстроперемен-ные течения. Поэтому система (1.14) с традиционным выражением для силы / пригодна только для течений с характерным размером L R (но, возможно, L <1 — межчастичного расстояния). Вряд ли реально найти / для общего случая.

Далее, однако, будет рассмотрена более простая задача устойчивости плавного течения. Если имеется решение с наложенным на него малым коротковолновым возмущением, то разумно предположить, что / определяется в основном плавным течением, так как возмущения / из-за ее интегрального по поверхности частицы характера будут тем или иным способом сглажены. Поэтому при рассмотрении устойчивости примем, что / зависит от параметров газа, каким-то способом усредненных на размере частицы. Подробнее вид /• рассматривается ниже.

Покажем, что основная причина некорректности — использование предположения о плавности течения. Действительно, система (1.1) не описывает коротких волн, и поэтому оценки инкремента на основании этой системы неприменимы. Например, в то время как естественное выражение для архимедовой силы (1.13) — интеграл от градиента давления по объему частицы, в (1.1) архимедова сила равна /д = —Vdpfdx, т. е. оставлен первый член ряда Тейлора. Очевидно, (1.13) будет сглаживать коротковолновые пульсации. То же относится к интегральным выражениям для ои G.

Так как kL 1, произведения возмущений на градиенты невозмущенного решения (типа v’dn/dx) опущены. Однако учитывается взаимодействие возмущения и силы G, которая в затухающем течении отнюдь не мала (второй член в правой части последнего уравнения (1.15)).

Примем простейшую зависимость р' = с2р'. Согласно (1.10), возмущениеR1.0 F 0.60.2 О О4Рис. 1.4.

Скорость распространения этих волн близка к скорости частиц. Мнимая часть Л4 соответствует затуханию возмущения при fw > 0. Это неравенство может нарушаться только в узкой области кризиса сопротивления. В этой области инкремент возмущения —rfw/pV содержит малость г и не зависит от к. За конечное время, пока частица находится в области кризиса, рост возмущения будет ограничен, а с уменьшением w оно начнет затухать. Эта неустойчивость вполне реальна: быстрые частицы в области кризиса сопротивления тормозятся слабее, чем медленные, и начальная флуктуация скоростей частиц растет.

Для корней Ai^ квадратной скобки в левой части (1.21) имеемReAi^ и ± с, то есть они по-прежнему приблизительно соответствуют звуковым волнам в газе. Выпишем мнимую часть в первом порядке по 7:1тА12 = 1Шк±-1Em^JL.

2 2с 2сС уменьшением w, практически уже при w с, первое слагаемое преобладает (для не слишком коротких волн). Поэтому сказанное о Л4 относится и к Ai2, с той разницей, что малый множитель г заменяется на (также малый) jFH/2. В любом случае инкремент ограничен, а при больших к он спадает по крайней мере как 1 /к2.

Для коротких волн основным является третье слагаемое. Растет амплитуда волны, распространяющейся в направлении относительной скоростиw = и — v. Нетрудно показать, что этот «рост» кажущийся. Суммарная скорость и + и' затухает.

При больших kR из-за спада D преобладает первое слагаемое под корнемэтот случай уже разобран. При kR 1 или г < 7 может преобладать второе слагаемое. Тогда ImAi, 2 < w (ry)ll2(kR)-^2 и за время проскальзывания рост возмущения также будет конечным.

В околозвуковом случае (|ги| «с), когда близки три корня левой части (1.21), полученными оценками пользоваться нельзя. Наиболее существенно влияние D, когда три корня практически совпадают. Тогда А2, з, 4 = А2, з, 4 + {D/2с)1/3 и мнимая часть и, по крайней мере, ограничена. Из-за малости D ограничен и рост волны. В существенно нестационарной задаче |ги| уменьшается со временем и для каждой частицы околозвуковая ситуация возможна в течение малого по сравнению с te времени.

Отметим, что для удобства мы считали параметры г и 7 малыми. В действительности это требование необязательно. Для малых и средних к при ограниченном времени проскальзывания рост возмущения ограничен. Для больших же к правая часть (1.21) мала за счет множителя F, что позволяет провести практически те же рассуждения. Таким образом, нелокальная модель решает проблему коротковолновой неустойчивости в достаточно общей ситуации.

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

Интересно, что скорость волны Re (du/dk) больше скорости звука в газе, что приводит к кажущемуся нарушению причинности. Однако такие возмущения не переносят информациирост амплитуды из-за неустойчивости превышает рост за счет перемещения волнового пакета. Для движения со сверхзвуковой скоростью прежде надо создать «затравку» волны на всем рассматриваемом интервале. Наблюдатель будет видеть возникновение «из ничего» сверхскоростной волны, в линейной постановке растущей за счет внешней поддержки стационарного течения. Подобные явления возможны и в других активных системах, причем скорость света также может быть превзойдена. Некоторые виды таких «тахионов» рассмотрены в работах [17, 18].

Обсуждение результатовПрактически для нелокальной модели начальное коротковолновое возмущение может возрасти в несколько раз. Поэтому исходное плавное течение при наложении малого возмущения изменится мало и задача Коши будет корректной.

Предложенная система может быть применена и для расчета течений с разрывами. Их устойчивость требует отдельного анализа. Однако разрыв, размазанный на расстоянии порядка R, для данной модели может рассматриваться как непрерывное течениетолько при определении G' нельзя будет выносить п и fx за знак интеграла. Это приведет к несущественному изменению вида формфактора в (1.18), причем все выводы, касающиеся поведения коротких волн, останутся в силе. Следует отметить, что физическая реальность модели на размерах порядка R определяется достоверностью зависимости для силы /.

Вместо интегродифференциальной системы на первый взгляд было бы проще использовать следующее по сравнению с (1.1) дифференциальное приближение. Это можно сделать, разложив в ряд Тейлора медленно меняющиеся функции в интегральных выражениях (1.10 — 1.13) и взяв первые неисчезающие поправки. Можно видеть, что это приведет к чисто дифференциальной системе с производными по координате вплоть до третьего порядка. Однако полученная таким образом система также приводит к некорректности. Действительно, при анализе устойчивости это равносильно замене формфактора F (kR) на его длинноволновую асимптотику 1 — (kR)2/10, которая приводит к неограниченному росту возмущений при kR 1. Отсюда ясно также, что дифференциальное приближение любого порядка для достаточно больших к будет неустойчиво и соответствующая задача Коши некорректна. Здесь мы встречаемся с ситуацией, принципиально отличной от однофазной газодинамики.

Отметим возможность обобщения метода на плавные трехмерные течения. Действительно, определения a, f^ и G легко записать в инвариантном виде как интегралы по объему частицы, не выделяя заранее какого-либо направления в пространстве. Переменные плавного течения можно понимать как среднеобъемныеони же будут средними по площадке, перпендикулярной распространению возмущения. Можно ожидать, что и в трехмерном случае сглаживание коротких волн обеспечит корректность задачи.

Результаты данного раздела изложены в нашей работе «Об уравнениях механики двухфазных сред», ПМТФ, № 6, 1983. Независимо близкий подход был опубликован в работах В. М. Фомина и С. П. Киселева [19, 20] применительно к течениям с резкими изменениями концентрации частиц.

Практическая численная реализация модели была бы достаточно сложной. По нашему мнению, она и необязательна. Поскольку установлена корректность задачи и близость нелокальной и локальной постановок для плавных течений, можно решать дифференциальные уравнения типа (1.1), а для их регуляризации применять известные методы, перечисленные в начале данного раздела. Именно так мы будем действовать при решении задач о конвективном горении (гл. 2).

Таким образом, двухфазная неустойчивость оказалась «ложной», так как поддается регуляризации путем уточнения межфазных соотношений. Она не ведет к появлению каких-либо мезоструктур (но установить это удалось, рассматривая мезопроцессы).

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

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

Постановка задачиПеред фронтом горения в двухфазной системе существует область течения без химической реакции (воздушная пробка). Трение между газом и частицами в этой области преодолевается напором свежих продуктов горения.

Приходим к следующей постановке задачи. В пористой среде по заданному закону движется «жидкий» поршень, проницаемый для частиц и непроницаемый для газа. Требуется найти движение газа перед поршнем.

Здесь q — теплообмен. В уравнениях учтены пульсационный перенос импульса (р5и2) и энергия пульсационного движения К. Обычно при расчетах этими флуктуациями пренебрегают. Однако, когда объемная доля твердой фазы не мала, газ движется по очень извилистым каналам, и флуктуации будут порядка самих величин. Поэтому в принципе учет пульсаций необходим во всех случаях, когда надо учитывать инерционные члены и кинетическую энергию среднего движения.

Очевидно, что изотермический случай соответствует достаточно малым скоростям, когда успевает пройти теплообмен. Это приближение хорошо известно в теории фильтрации. Наоборот, пренебречь теплообменом можно при больших скоростях течения.

Изотермическая задачаПри постоянной температуре р = с2р (с — изотермическая скорость звука). Эту связь следует использовать вместо уравнения энергии. Для сравнительно небольших скоростей имеет смысл рассмотреть двучленный закон сопротивления:/ = CDu2/d + hu.

Проверка выполнения предположения (1.26) на полученном решении дает д Гтитт U{U + a) дтноткуда В = 1/(1 + а). Тот же результат получается из условия сохранения массы сгребаемого газа: Lоо (R-l)d? = T ос точностью до слагаемых 0(1).

На рисунке 1.5 приведены профили U и R при, а = 0,01 (малый вклад линейного сопротивления) для моментов времени г = 10- 102- 103 (кривые 1−3 соответственно). С течением времени профиль скорости, сохраняя в основном свою форму, «логарифмически» сдвигается вперед, так что размер возмущенной области близок к In т. В этой области плотность падает9экспоненциально<�координатой и растет пропорционально времени.

Учет инерционных слагаемых в уравнении импульса приводит к появлению скачков, если скорость волны превышает с. Сохранение изотер-мичности при таких условиях сомнительно, и этот случай подробно не разбираем.

Нулевой теплоотвод Рис- 1−5.

Для системы (1.22) ищем решение в виде волны, в которой р = const, и = v. Предполагаем, что пульсационные слагаемые — функции средних скорости и и плотности р, например, согласно (1.23). Первое уравнение выполняется автоматически, и система упрощается: Рх = -CDpv2ld, pt + 7 vpx = 0.

Видно, что скорость фронта волны D = jv, а зависимость р (х) линейна. Плотность и давление в волне выражаются формуламиР = —^Ро, р = ^rPQCDv2lVt' Х + Ро ¦ (1.27)7—1 7—1 dРешение существует в области между поршнем х = vt и фронтом волны х = fvt. Максимальное давление на поршне линейно растет со временем. Отметим сходство (1.27) с классическим газодинамическим течением — ударной волной. Здесь также имеем скачок, создающий область постоянных плотности и скорости, но вместо давления постоянен его градиент.

Продолжение формул (1.27) до х = 0 дает решение задачи о вдуве газа в полупространство с постоянным расходом. Стационарное решение можнополучить, если двигать со скоростью yv поршень, частично проницаемый для газа.

Решение (1.27) применимо для достаточно длинных волн (vt d). Вблизи фронта существует переходная область шириной порядка диаметра частицы, в которой происходит перестройка течения от состояния покоя к установившемуся решению (1.27). При дозвуковой скорости фронта волны переход плавный, а при сверхзвуковой возникнет скачок. Подробное исследование этой ситуации нецелесообразно, пока не установлено зависимости для флуктуаций в зоне резких изменений течения. В этой же области вырабатывается постоянное слагаемое ро Pqv2, которым можно пренебречь для длинных волн. Результаты этого пункта с очевидными изменениями сохраняются для любого закона сопротивления вида f (p, u).

Изотермическое решение представляет определенный интерес для классических задач фильтрации. Адиабатическое решение применяется в следующей главе к анализу конвективных волн.

1.2.2 Компактирование непрочной пористой среды газовым поршнемРассматривается воздействие газа, находящегося при высоком давлении, на прилегающую пористую среду. Газ, проникающий в поровое пространство, вовлекает в движение частицы среды, что приводит к ее сжатию. Ниже рассмотрен предельный случай среды низкой прочности, для которой в определенный момент достигается практически полное уплотнение. Динамика фильтрации и сжатия пористого скелета исследована аналитически (на качественном уровне) и численно.

Пример явления, для которого существенно взаимодействие фильтрационного течения и компактирования, — инициирование пористого взрывчатого вещества газовой детонацией или электрическим разрядом. Такие же процессы характерны для начальных стадий взрыва в пористом грунте [22, 23, 24]. Подобная ситуация возможна и при взрывном компактировании порошков.

Постановка задачиВ начальный момент времени t = 0 полупространство х < 0 занято газом с начальным давлением ро и плотностью ро • При х > 0 располагается среда с начальной открытой пористостью <ро • Давление газа в порах мало по сравнению с ро. Исходный диаметр частиц среды d, плотность частиц ps значительно превышает ро. Прочностью частиц пренебрегается, т. е. порошок уплотняется без сопротивления.

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

Затем смещение частиц под действием потока газа становится существенным. В отсутствие прочности оно приводит к образованию компактной «пробки» вблизи поверхности пористой среды. По крайней мере, в начале разгона частиц можно находить их движение, считая состояние газовой фазы известным из решения фильтрационной задачи. Разумеется, при заметных деформациях среды такое приближение не работает, но для оценок интегральных параметров, таких как время компактирования, оно будет пригодно.

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

Стадия фильтрацииНа этом этапе пренебрегаем движением твердой фазы, и достаточно рассмотреть уравнения для газа. Задача эта, хотя и существенно нестационарна, во многом подобна рассмотренной в предыдущем пункте 1.2.1. Как уже говорилось выше, система типа (1.22) сводится к более простой системе (1.24). Так как сопротивление трения в пористой среде оказывается преобладающим фактором, в уравнении импульса газа допустимо отбросить переносные члены и считать, что градиент давления уравновешивается силой трения: <�рдр/дх = — /.

Действительно, слагаемые dp^pu/dt и dpipu2/dx имеют порядок pipu2/L, где L — характерный масштаб течения. Сила трения / Сври2/(1. Осредненное описание течения возможно, когда L d. Следовательно, оба инерционных слагаемых малы по сравнению с /.

Кроме того, ввиду приближенности подхода, ориентированного на оценки, разумно дальнейшее упрощение задачи: вместо уравнения энергии газа будем использовать простейший адиабатический закон р р1. При этом 7 следует понимать не как показатель адиабаты газа, а как эмпирический параметр, варьирование которого позволяет перекрыть область возможных состояний газа.

При 7 = 1 асимптотика получается более сложной. Прямым вычислением можно проверить, что в этом случае вблизи нуляПри этом производная дР/д? имеет логарифмическую особенность.

Для расчета движения частиц понадобится переход к эйлеровым коорtдинатам. Из выражения для координаты х = r+f u (r, t) dt легко получитьоформулу пересчета для автомодельной эйлеровой координаты г] = kx/t2у/3 :dPdt,= J Plhdi.

Из приведенных выше асимптотик вытекает, что максимальное значение Tjj — г](£ = 0), т. е. автомодельная координата фронта ту = 2fA. Следовательно, фронт волны фильтрации движется по закону х = 2[А^Ъ jk. При 7 = 1 имеет место расходимость вида т] — 1/f). Практически при численных расчетах этот особый случай не вызывает затруднений, так как расходимость крайне слабая. Графики Р{г)) показаны на рис. 1.6,5.

Функция F{tj) имеет особенность вида ту-½ вблизи нуля. Она происходит из такой же особенности производной безразмерной скорости dV (jj)/dr), в которой легко убедиться, исходя из (1.32). На рис. 1.7 показаны графики V{t]) и G{t]) = y/rjF (tj) для крайних значений 7 = 1 и 3. Видно, что в этом интервале 7 при 77 О F «0,15/^.

Ч> оили trdtpoaops (1.34)4 psOi fd. c «сл/CdPoВидно, что реально зависимость времени tc от начального давления (т. е. плотности ро) обратно пропорциональная.

Интерес представляет также количество газа, проникающего в поры до момента коллапса. На единицу площади границы раздела приходится масса тс = poR (tc) = Ро? о&3/к. В принятых приближенияхЧисленное моделирование.

Аналитический подход по необходимости носит качественный характер. Более детальное представление о процессе дают прямые численные расчеты.

Система уравнений (2.5) (подробно описанная в следующей главе), без учета массообмена фаз, с начальными условиями, описанными выше, решалась в эйлеровых координатах методом Лакса-Вендроффа, модифицированным для учета недивергентных и алгебраических членов с сохранением второго порядка аппроксимации. Шаг по времени выбирался из условия Куранта, а также условия стабильности по правым частям. Слабая двухфазная неустойчивость подавлялась сглаживанием параметров течения, что равносильно небольшой искусственной вязкости [25, 26].

С достижением максимального сжатия зона компактирования становится непроницаемой. Начинается рост ее толщины за счет подпора газа (рис. 1.8,б, в). Почти весь градиент давления приходится на область разгона частиц перед волной. Оценка из законов сохранения скорости волны сжатия D на этом этапе дает D = -y/p/psao^o > что хорошо совпадает с полученной в расчете.

Стадия чистой фильтрации оказывается непродолжительнойпрактически фильтрация и компактирование не разделены во времени. Несмотря на это, оценка хорошо согласуется с численным расчетом: из (1.34) время tc = 4,8 мкс. В расчете за время коллапса можно принять начало формирования непроницаемой пробки, когда, а вырастает до максимального значенияэто происходит приблизительно к моменту t = 5 мкс (рис. 1.8,6). С учетом предыдущих замечаний можно сказать, что упрощенная двух-стадийная модель с разумным приближением описывает передний фронт волн v (x, t) и a (x, t) и движение газа до момента коллапса.2оУчет давления скелета не вносит заметной разницы, если при максимальном сжатии ps < ро. С другой стороны, учет теплообмена фаз необходим. При выключении теплообмена q температура газа растет из-за трения до нереальных значений. Это, в свою очередь, изменяет картину течения: фильтрация идет быстрее, а компактирование подавляется. Интересно, что простая модель (1.30) гораздо ближе к реальности, чем более полная система (2.5) с единственным недостатком — отсутствием теплообмена (хотя, казалось бы, это упрощение естественно для «адиабатического» уравнения состояния).

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

1.3 Неустойчивость вытесненияВ механике многофазных течений детерминированный континуальный подход, рассмотренный выше, далеко не всегда дает удовлетворительные результаты. В данном разделе рассмотрены задачи, в которых существенна случайная компонента течения в пористой среде, и свойственные таким течениям неустойчивости.

1.3.1 Вытеснение плотной жидкостиВажный пример такой ситуации — вытеснение жидкости из пористой среды. Течение предполагается быстрым, так что число Рейнольдса, вычисленное по характерному размеру пор, велико. Если вытесняющий флюид менее плотный (например, газ), то граница раздела неустойчива. Ситуация аналогична хорошо известной неустойчивости вязкого пальца (viscous finger instability), с тем отличием, что управляющим параметром является не вязкость, а плотность. Здесь мезопроцессы в порах непосредственно формируют мезоструктуры, во многом подобные классическим вязкимпальцам.

Вязкие пальцыИмея в виду сравнение с классической проблемой вытеснения, вначале опишем вязкий режим. В пористой среде имеется граница раздела между вязкой жидкостью (например, нефть) и сравнительно невязкой (вода). Особый интерес представляет предельный случай, в котором можно пренебречь вязкостью толкающей жидкости. В вязкой области при малой скорости течения и выполняется линейный закон фильтрации: К. и =-Vp.(1.36)Здесь р — давление, К — проницаемость среды, у. — вязкость жидкости. Из условия несжимаемости V • и = 0 следует уравнение Лапласа: Ар = 0.(1.37)В пределах невязкой области давление однородно. Следовательно, и на движущейся границе раздела давление постоянно. Это граничное условие позволяет найти решение уравнения Лапласа в вязкой области. Тогда скорость движения границы определяется из (1.36), что позволяет рассчитать процесс вытеснения.

Если вытесняется более вязкая жидкость, то граница раздела подвержена неустойчивости Саффмана-Тейлора [27]. Когда на границе появляется выступ (рис. 1.9), то в его окрестности градиент давления увеличивается, растет скорость течения, а вместе с ней и скорость границы. В результате выступ разрастается. «Пальцы» невязкой жидкости прорываются вперед. РисВо многих опытах, начиная с описанных в [27], для изучения вязких пальцев использовалась ячейка Хеле Шоу. Соответствие этой модели, представляющей собой течение в узком зазоре между двумя параллельными плоскостями, движению в пористой среде нельзя назвать полным. В ячейках Хеле Шоу образуются широкие, гладкие пальцы, так как коротковолновая неустойчивость подавлена поверхностным натяжением. В прямоугольной ячейке лидирующий палец препятствует росту своих соседей. В результате выживший палец заполняет примерно половину ячейки. Патерсон [28] показал, что в радиальной ячейке Хеле Шоу пальцы разветвляются, достигнув критической ширины. Эта ширина также определяется поверхностным натяжением, и форма пальцев по-прежнему гладкая.

В работах Федера и др. [29] использовалась более реалистичная модельслой стеклянных шариков между стеклянными пластинами. В этой ячейке получены совершенно другие картины вытеснения (рис. 1.10). Здесь воздух вытесняет эпоксидную смолу из модельной пористой среды. Газ подается через отверстие в центре структуры, так что растекание преимущественно радиальное. Поперечник пальца определяется размером пор. Пальцы разрастаются, вблизи концов выпускают ветви, и таким образом развиваются причудливые разреженные структуры. Вытеснение в прямоугольной ячейке исследовалось в [30], и были получены сходные формы. В тех случаях, когда удается исключить поверхностное натяжение, подобные картины наблюдаются и в ячейках Хеле-Шоу [31]. Эти результаты подтвердили аналогию, отмеченную в [32], между вязким вытеснением и диффузионно-лимитированной агрегацией (DLA). Концепция DLA, введенная Виттеном и Сандером [33, 34], рассматривает образование кластера из малых частиц, диффундирующих к поверхности кластера. Стационарная диффузия описывается тем же уравнением Лапласа, а граница кластера растет со скоростью, пропорциональной градиенту концентрации. В процессе DLA образуются весьма разветвленные структуры, которые с успехом воспроизводятся в стохастических моделях роста.

Аналогичные задачи возникают при исследовании электрического пробоя, с заменой давления на потенциал. Роль невязкой жидкости играет эквипотенциальная система каналов разряда. Применение стохастических моделей DLA к пробою предложено в [35].

Сейчас линейные задачи типа (1.36 — 1.37) достаточно исследованы. Во всех этих ситуациях неустойчивость границы ведет к развитию сложных (фрактальных) структур. В частности, кластер, образованный струями невязкой жидкости, описывается основным соотношением: М RD. ЗдесьМ — масса кластера (или занятый им объем), R — его характерный размер. Показатель D — фрактальная (дробная) размерность, которая меньше размерности вмещающего пространства. Из-за этого большие кластеры выглядят разреженными. Например, размерность структуры на рис. 1.10 равна 1,62.

Сравнительно мало внимания в литературе уделялось нелинейным задачам. В работах Нимейе-ра, Пьетронеро и др. моделировался электрический пробой в предположении, что скорость разрастания границы разряда пропорциональна электрическому полю в некоторой степени 77. При нелинейном («Л Ф 1) условии роста на границе в пространстве по-прежнему решалось уравнение Лапласа. В раРис. 1.10.боте [31] рассматривалось вытеснение неньютоновской жидкости, для которой и |Vp|m, где 77i > 1.

Невязкие пальцыОбратимся к противоположному предельному случаю быстрой фильтрации. Когда число Рейнольдса, вычисленное по характерному размеру пор, велико, вместо линейного соотношения (1.36) выполняется квадратичный закон сопротивления [36, 37] (режим Форхгеймера):pu|u| = -fVp. (1.38)Здесь р — плотность жидкости, / - постоянный коэффициент, по порядку величины близкий к размеру пор. В случае несжимаемой жидкости из уравнения неразрывности получаем для давленияV- (Vp|Vp|" ½) =0. (1.39)Если жидкость вытесняется газом с малой плотностью, то в газовой области градиент давления мал, и в первом приближении на границе раздела давление постоянно. Нормальная скорость границыu|Vp|½. (1−40)Качественно нелинейная задача (1.39 — 1.40) аналогична вязкой постановке, и следует ожидать неустойчивости вытеснения: роста выступов и их разветвления. По аналогии с неустойчивостью вязкого течения можно говорить о неустойчивости «плотного пальца», или о невязкой неустойчивости (inviscid finger instability).

Ситуации, в которых течение достаточно быстрое для применения простейшего подхода (1.39 — 1.40), довольно обычны. Для нас существенно, что до некоторой степени эта постановка описывает фильтрацию продуктов горения в пористом энергетическом материале. Высокая температура продуктов соответствует малой их плотности, а роль вытесняемой жидкости играет газ (воздух) в поровом пространстве. Неустойчивость границы раздела приведет к распространению пламени в виде прорывающихся вперед струй. Очевидна важность таких явлений для развития процесса.

Линейная стадияРассмотрим начальную (линейную) стадию неустойчивости границы. Пусть форма границы задается уравнением (рис. 1.11)y = h cos (fcz), (1−41)причем амплитуда возмущения мала: kh 1.

Вытесняемая жидкость находится выше линии (1.41). В невозмущенном течении задан вертикальный градиент давления, по величине равный |Vp|0. Представляя давление в виде р = ро — у | Vp|0 + р', (где ро приложенное давление, р' - малая добавка), в первом порядке по р' из (1.39) имеем- 1^ = 0 дх2 2 ду2Это уравнение сводится к уравнению Лапласа изменением масштаба по оси У в раз.

Таким образом, для начальной стадии неустойчивости переход к уравнению (1.39) не вносит качественных изменений. Как будет видно из дальнейшего, направление количественного влияния нелинейности — относительное замедление роста выступов и увеличение градиентов возмущений вблизи растущих участков — сохраняется и для развитой неустойчивости.

Развитые возмущения в нелинейной задачеБолее интересна и важна для формирования картины течения развитая нелинейная стадия роста возмущений. В качестве модельной постановки рассмотрим палец высоты h и поперечного размера 2 г, который существенно меньше h (рис. 1.12). Оценим увеличение скорости роста конца пальца по сравнению с фоновой скоростью для ровной поверхности.

Вначале рассмотрим ситуацию в двух измерениях. В линейной постановке характерный размер? спадания поля Vp вблизи вершины пальца оказывается порядка yfkr (а не г, как может показаться на первый взгляд).

График F (h/r) приведен на рис. 1.13. Здесь же показаны результаты численных расчетов. Крестики — расчет уравнения Лапласа по стандартной разностной схеме методом установления, точки — моделирование диффузии.

Последний метод опишем подробнее. На квадратной решетке в плоскости XY задается распределение частиц, каждая из которых за шаг по времени совершает скачок единичной длины в одном из четырех направлений. На верхней границе задается постоянная концентрация. На нижней, включающей выступающий палец — условие поглощения.

По истечении некоторого времени релаксации 5 концентрация частиц будет подчиняться уравнеFнию Лапласа. Вычисляя поток частиц на верхушку пальца, получим множитель F. В пределах ожидаемой погрешности, основные источники которой — дискретизация и конечные размеры области, наблюдается совпадение расчетов и теории. Практически диффузионный расчет оказался менее быстрым.-Л Л* т^-х— <Хf — О8 h/r 20 Рис. 1.13.

Фактор роста для нелинейной задачиУравнение (1.39) решалось численно в области, показанной на рис. 1.12. Детали вычислений изложены в работе [38]. Результаты расчета нелинейной задачи разностным методом показаны на рис. 1.14. Здесь коэффициент /3 равен 1 для двумерной задачи и 0,81 — для трехмерной.

Для плоского случая экспериментальные точки хорошо аппроксимируются зависимостьюF = (1 + h/r)0'6.

Для линейной задачи (уравнение Лапласа) получен показатель степени 0,491 при теоретическом значении 0,5. Поскольку фактор увеличения градиента имеет порядок (/i/r)0,6, характерный размер? спадания поля пропорционален h0,4r0'6, то есть меньше, чем в линейном случае, а скорость роста верхушки и wo (^/r)0'3.

Для трехмерного случая кажется разумным принять, что вблизи конца пальца единственный характерный размер — это поперечник г. Результатыln (l + J3h/r) Рис. 1.14.численного моделирования согласуются с таким предположением. Хорошая аппроксимация получается приF = (1 + 0.81/i/r)1'01.

Тогда скорость роста и щ (Ь/г)0,ъ.

Таким образом, для скорости роста ип = —Ay/Vp показатель степени фактора роста в плоском случае равен 0,3, а в трехмерном случае 0,5.

Эти зависимости сравним с линейной моделью. В электростатике — возрастание поля, как известно, в основном пропорционально fhjr в двумерном случае и /г/г, со слабой логарифмической поправкой — в трехмерном. В нелинейной модели пробоя диэлектрика [35] при степени т] = 0,5 получаем практически такое же ускорение роста: (h/r)0,25 и (h/r)0,5 соответственно.

Если допустить, что фрактальные свойства течения в основном определяются фактором роста ветвей, то известные данные (т] = 0,5) — модели можно было бы использовать как оценку для нелинейной задачи. Однако такое приближение не работает, по крайней мере в плоской геометрии.

Последние два результата согласуются с [35, 39, 40]. Наиболее близкая процедура применялась Эвертцем [39], и наши данные хорошо соответствуют его результатам (Do, 5 = 1,89, D = 1,64 для решетки размером 64-£>о, 5 = 1,918, D = 1,663 для бесконечной решетки). Это показывает, что систематические ошибки расчета незначительны. Кроме того, возможные неточности гасят друг друга при рассмотрении различий между вариантами а) — с). Поэтому можно утверждать, что «нелинейная» размерность 1,78 четко отличается от обоих контрольных значений Z? o, 5 и D.

В случае (Ь) по сравнению с © растущие кластеры имеют более толстые ветви и, соответственно, их фрактальная размерность больше. Это естественно, так как для малых т] контраст вероятностей на растущей границе сглаживается. В пределе 77 = 0 получается модель Идена (рост границы в случайно выбранном месте), для которой D = d.

Нелинейный случай (а) как по значению размерности, так и по внешнему виду кластеров — промежуточный между (Ь) и ©. Различие с линейным вариантом — естественное следствие относительно более медленного роста верхушек пальцев.

В обоих случаях длинные пальцы растут практически с одинаковой скоростью. Однако в (Ь) поле сильнее проникает в растущий кластер (меньше экранируется) и возникающая структура несколько плотнее. Распределение вероятностей роста в «полулинейной» модели — это просто сглаженное Лапласово распределение: модель не меняет характерного пространственного масштаба Для нелинейного поля распределение более контрастно из-за большего градиента на выступах и меньшего — в «фиордах». Таким образом, фрактальная размерность зависит не только от степени ускорения роста пальцев, но и от всего распределения поля на границе.

1.3.2 Эксперимент. Влияние анизотропииОбразуются ли в реальном быстром течении структуры возрастающей сложности, как это происходит с вязкими фракталами [29], или доминирует противоположный тип неустойчивости, как в [27], когда геометрия течения прогрессивно упрощается? Ниже представлены экспериментальные результаты для нелинейного вытеснения, а также теоретические оценки. Особое внимание уделяется вытеснению в анизотропной среде. Анизотропия почти неизбежна для двумерной модели пористой среды, представляющей собой набор каналов и препятствий: в двумерном слое легче всего возможна визуализация. Важное значение имеет анизотропия для процесса горения. Распространение пламени в решетках каналов рассматривалось, например, в [41]. Наконец, любая реальная пористая среда в малом масштабе анизотропна.

Ячейка заполнялась жидкостью. В ресивер, соединенный с камерой, накачивался воздух. Ресивер был снабжен клапаном в виде разрываемой мембраны. После открытия клапана сжатый воздух проникал в камеру и начинал вытеснять жидкость. На нижней границе камеры отверстия соединяли ячейку с объемом, заполненным жидкостью и закрытым пробкой. Когда начиналось вытеснение, пробка вылетала и жидкость могла выходить наружу. Сечение отверстий было достаточно, чтобы пренебречь их гидравлическим сопротивлением. Избыточное давление воздуха (0,5 — 1 атм) более чем на порядок превышало давление столба жидкости в ячейке. Деформация камеры внутренним давлением предотвращалась жестким креплением достаточно толстых (30 мм) стенок. Во избежание возможных утечек сборка тщательно уплотнялась.

Капиллярный скачок давления на поверхности раздела составлял десятые доли процента от избыточного давления. Однако капиллярные эффекты могли быть существенны в некоторых особых случаях. Эти ситуации подробно рассмотрены ниже.

Скорость течения получалась порядка метров в секунду. Число Рей-нольдса, вычисленное по поперечнику канала или по характерному размеру поры, порядка нескольких тысяч. Это значит, что закон сопротивления близок к квадратичному [42]. Ширина ячейки I = 20 см, длина L — 50 см. Объем ресивера (1 л) позволял поддерживать давление на границе, близкое к постоянному (например, для слоя малых шариков спад давления менее 10% при вытеснении половины столба жидкости). Это условие нарушалось для слоя шариков диаметра 12,2 мм, так как объем ячейки был гораздо больше. В этом случае при устойчивом одномерном вытеснении была бы практически постоянна скорость течения (градиент давления), т. е. падение давления приблизительно компенсировалось бы уменьшением высоты столба жидкости.

Форма границы регистрировалась фото-, киноили видеосъемкой, в отраженном свете либо на просвет. В качестве рабочей жидкости использовались подкрашенный спирт (в опытах со стальными шариками) или подкрашенная вода.

Результаты экспериментаНа рис. 1.17 показана кинограмма одного из опытов с квадратной решеткой каналов. Порядок кадров в каждом ряду слева направо, нижний ряд является продолжением верхнего. Скорость съемки 150 кадров/с. Избыточное давление воздуха 0,7 атм. Ориентация решетки — продольная (схематически она показана на последнем кадре). Воздух виден как светлая область, Рис. 1.17.а заполненные водой поры — темные. Вначале вытеснение довольно равномерное. Примерно через 1/15 секунды (10 кадров) на фронте вытеснения хорошо видны неоднородности, которые развиваются в продольные струи. Количество и форма струй меняются от опыта к опыту.

Анализ снимков показал, что при продольной ориентации решетки жидкость не только остается в промежутках между струями, но не вытесняется также из поперечной подсистемы каналов (течение имеет выраженную вертикально-полосатую структуру). В этом заметную роль играет поверхностное натяжение. Для горизонтального направления относительный масштаб капиллярных эффектов дается отношением перепадов давлений:(Аа/5) /(ddp/dx|) где сг — коэффициент поверхностного натяжения. Для вертикального градиента др/ду при характерных условиях опыта эта величина была бы порядка 1/10. Поскольку горизонтальный градиент значительно меньше, поверхностные мениски способны удерживать жидкость в поперечных каналах между соседними вертикальными, уже заполненными воздухом. Следуетотметить, что это не влияет на течение в области чистой жидкости, где нет границ раздела.

Интересно, что «застрявшая» жидкость может влиять на результат последующего эксперимента. На рис. 1.18 показан случай, когда над уровнем сплошного заполнения имелась область частично смоченных пор (видная как полоса промежуточной яркости под светлой сухой областью на первом кадре). Эта полоса без заметных изменений видна б течение всего процесРис. 1.18.са, то есть жидкость в поперечных порах в основном остается на месте. Однако в таких ситуациях узкие пальцы развиваются гораздо быстрее (ср. с рис. 1.17). По-видимому, остаточная жидкость создает дополнительные возмущения, например в виде локальных препятствий для воздуха.

На рис. 1.19 показан типичный результат «выстрела» с диагонально ориентированной решеткой, Здесь затруднен примитивный прорыв воздуха вдоль канала, так как каждый канал упирается в стенки. Остальные условия опыта те же, что на рис. 1.17 и 1.18. Диагональное вытеснение боРис. 1.19.лее медленное, однако искажения фронта возникают гораздо раньше. Из них также развиваются струи. Форма струй не удлиненная, а скорее клиновидная, угол при вершине клиньев приблизительно прямой. Эта тенденция усиливается с ростом амплитуды возмущений.

Эти результаты резко отличаются от тех, которые наблюдались при медленном вязком режиме вытеснения, когда образуются ветвистые структуры [29. 31, 30]. Причина различия не в нелинейности сопротивления (которое влияет только на степень разветвления, [38]), а в анизотропии среды. Поэтому ряд экспериментов был проведен с другой моделью пористой среды — слоем шариков. Эта модель представляется сравнительно более изотропной (во всяком случае, характер ее анизотропии явно другой).

На рис. 1.21 показаны «статический» кадр (слева) и одна из стадий течения (справа) в слое шариков диаметра 2,5 мм. Воздух (темная область) перемещается в кадре сверху вниз. Избыточное давление воздуха 0,7 атм. Расположение ячейки вертикально. Стрелка показывает направление течения. Светлые полосы, видные в верхней части кадров — препятствия, выравнивающие давление воздуха над жидкостью.

Рис. 1.21, Рис. 1.22. Рис. 1.23.

Это подтверждается экспериментами с более крупными шариками (рис. 1.22). Здесь изображены две различные стадии течения. Съемка на просвет, темные пятна соответствуют прохождению света по «воздушным» порам. Видно, что прозрачная область не сплошная. Видеосъемка показала, чтооткрывшиеся поры обычно закрываются, затем опять открываются, и т. д.

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

Здесь форма области вытеснения, по крайней мере качественно, близка к той, которую можно ожидать для изотропной среды. На границе раздела жидкости и воздуха видны типичные ветвящиеся струи. При съемке той же ячейки снизу картина совершенно иная — наблюдается небольшое количество узких струй. Различие объясняется тем, что при наклонном или горизонтальном положении ячейки воздух преимущественно распространяется по верхней части слоя.

В известных опытах [30, 31], соответствующих линейноТаким образом, качественная аналогия с вязким режимом вытеснения частично подтвердилась. Наблюдается неустойчивость фронта и развитие струй. Вместе с тем, течение в плоских модельных ячейках оказалось более разнообразным, чем можно было ожидать. Течение носит «отпечаток» симметрии анизотропной структуры модели. Только в шаровом слое при почти горизонтальном расположении, когда эффект плавучести блокирует переключение струй между стенками, удается получить картины, напоминающие фрактальные кластеры вытеснения (хотя для вязкого режимаму режиму трения, получаются более разреженные структуры. Для сравнения на рис. 1.24 изображена картина медленного вытеснения очень вязкого масла воздухом из работы Ленормана и др. [30]. Течение происходит в модельной пористой среде, инжекция газа — вдоль линии (нижней границы).

Рис. 1.24.течения характерны гораздо более разреженные структуры с тонкими ветвями [29, 30]). Не исключено, что течение в вертикальном шаровом слое сходно с почти горизонтальным случаем, а наблюдаемые различия — чисто внешние и связаны со склонностью струй частично прятаться на невидимой стороне ячейки. В горизонтальном течении эта тенденция подавляется плавучестью.

Модель анизотропного теченияОграничимся простейшим случаем квадратной решетки. Вначале рассмотрим продольную ориентацию, когда каналы направлены по осям X и Y. Здесь и далее, поскольку сила тяжести в расчетах не учитывается, удобнее рассматривать направление вытеснения вдоль оси Y снизу вверх.

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

При типичной скорости течения > 1 м/с число Рейнольдса, рассчитанное по ширине канала, превышает 2000. Такое течение является достаточно быстрым для установления турбулентного режима. Оценка по формуле Пуазейля показывает, что при ламинарном продольном течении развивались бы скорости более 20 м/с. Поэтому при анализе учитывалось только квадратичное («инерционное» по фильтрационной терминологии) трение.

Существенно, что средние вертикальные скорости во входном и выходном сечении различны. У втекающей из горизонтального канала жидкости вертикальная скорость практически нулевая, а у вытекающей — порядка V. Чтобы поддерживать такой вынос импульса из контрольного объема, и необходим перепад давления Ар. Войдя в горизонтальный канал, жидкость потеряет вертикальный импульс (передаст его стенке).

Согласно системе (1.43), градиент давления не антипараллелен скорости. В общей ситуации, когда компоненты одного порядка, первые (перекрестные) члены вносят заметный вклад в сопротивление. При Ъ > а они отклоняют градиент давления к биссектрисе квадранта, противоположного тому, в котором лежит вектор скорости, а при Ъ < а градиент «прижимается» к оси координат, ближайшей к направлению скорости. Следует отметить, что анизотропия сохраняется и при, а = Ь. Хотя в этом случае градиент точно противоположен скорости, но величина сопротивления зависит от направления течения (максимально сопротивление при течении вдоль диагонали).

На первый взгляд, анизотропию среды можно было бы учесть, введя тензор сопротивления, компоненты которого зависят от модуля скорости (по аналогии с линейным анизотропным случаем, когда применяется тензорное обобщение закона Дарси [36, 37]):Vpi = -RikVk.

Однако для квадратной решетки указанный тензор сводится к скаляру. Действительно, при естественном выборе осей — вдоль каналов решетки — матрица сопротивления должна быть диагональна (иначе при течении вдоль осей градиент давления имел бы составляющую поперек скорости). Из равноправия направлений х и у следует, что диагональные элементы равны и тензор Rik — RSik, где R — скаляр. Поскольку единичный тензор не меняется при вращениях системы координат, то при любой ориентации осей такой закон сопротивления имел бы вид Vpi = —Rvi, т. е. анизотропия полностью выпадает.

При продольном вытеснении можно связать градиент давления со скоростью струй V: AP/z = 0,62apV2/S. Эффективный коэффициент сопротивления каналов 0,36а = 0,13. Это в несколько раз больше, чем можно ожидать из известных гидравлических корреляций для гладких труб (0,02 при квадратном поперечном сечении [42]). По-видимому, заметный вклад вносит «проницаемость» стенок — поперечные каналы, на срезе которых возможен обмен массой и импульсом даже в отсутствие макроскопического поперечного течения.

Сравнительный масштаб сопротивления из-за поперечной инжекции массы (6 — члены в (1.43)) и из-за собственно трения дается отношением b/а «0,75. Это отношение можно считать мерой анизотропии решетки.

Относя это условие к оси х, находим решение (1.45):р = |Vpo| hcos (kx) ехр (— у/ + Ь/а • ку).

Возмущение в вертикальном направлении спадает более быстро по сравнению с вязким течением (ехр (-ку), [27]), как и в изотропном случае при квадратичном законе трения (ехр (—л/2ку), см. выше).

Подобным же образом для продольного течения из (1.43), предполагая и v, получимдрV (а (<Э</ +½ «=' bvdxдхСнова используя уравнение неразрывности, получаем в первом порядке по возмущению р'(L46)Очевидно, что анализ этого уравнения значительно сложнее из-за присутствия неаналитичного коэффициента (sgn[др/дх), равного +1 или —1 в зависимости от знака др/дх). В (1.46) с тем же успехом можно написать sgn [др'/дх), поскольку в невозмущенном течении др/дх = 0.

Регуляризация уравнения (1.46) требует дополнительных, довольно ограничительных предположений. Кроме того, нет уверенности, что (1.46) правильно описывает неустойчивость продольного течения, поскольку из-за явной неполноты вытеснения процесс более сложный. Поэтому в продольном случае инкремент получить не удалось.

Значение 9Х близко к ожидаемому, особенно если учесть грубость модели и значительный разброс экспериментальных данных. К сожалению, данные не позволяли провести более информативную обработку. Во-первых, спектральный интервал, допускающий такое рассмотрение, не особенно широк. Применимость континуальной модели ограничена условием kd = 2mrd/l < 1, что соблюдалось при п < 5. Во-вторых, фаза линейного роста непродолжительна, что затрудняет более точное определение скорости роста.

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

При продольном вытеснении на оси растущего в направлении у пальца сопротивление сводится к чистому трению: др/ду — av2. Разрастанию же боковой границы пальца в направлении х противодействует, кроме трения, поперечный эффект: др/дх & —bVu — аи2. Скорость V индуцируется движением «пьедестала» пальца, то есть основным фронтом вытеснения, и имеет тот же порядок, что скорость роста верхушки. Кроме того, поперечный градиент давления в общем меньше продольного, и скорость и течения в боковом направлении мала. Компоненты градиента одного порядка только вблизи верхушки пальца, заметно обогнавшего своих соседей. Из-за сравнительной медленности развития продольной неустойчивости такая ситуация для каждого из поперечных каналов непродолжительна. В результате заметного бокового разрастания пальца не происходит.

При диагональной ориентации рост неустойчивости быстрый, и выступы быстро достигают клиновидной стадии. После этого их форма становится гораздо более стабильной, так как границы имеют движение, близкое к продольному. (Следует заметить, что это распространение не в точности продольное, поскольку жидкость вытесняется полностью, см. рис. 1.20). Из-за общего смещения границы на ребрах клина не успевают развиться неустойчивости, так как «старые» участки поглощаются границей. Эксперимент показывает, что клинья не имеют заметной тенденции к ветвлению. Скорее наоборот, они сливаются в более крупные образования.

В конечной стадии формируется один большой выступ в центральной части ячейки. Это можно объяснить влиянием границ. Для боковых участков сопротивление больше, так как каналы упираются в стенки и течение вынуждено менять направление. В центре нижней части ячейки каналы «чувствуют» присутствие выходных отверстий, от чего сопротивление уменьшается.

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

Расположение слоя шариков было ближе к диагональному случаю. Однако из-за большей геометрической сложности порового пространства движение в этой модели более разнообразно по сравнению с регулярными решетками, в частности, из-за бифуркаций струй между стенками. Количество остаточной (не вытесненной) жидкости заметно больше, чем в диагональной решетке. В свете данных, полученных с квадратной решеткой, изаметной анизотропности шарового слоя даже некоторым сюрпризом выглядит картина течения на рис. 1.23, имеющая «почти» фрактальный вид. Можно полагать, что в пристеночном слое анизотропия шаровой модели не так велика, чтобы нарушить естественный для изотропного случая процесс ветвления при вытеснении. Чтобы ветви были видны, течение должно быть прижато к плоскости наблюдения.

ВыводыТаким образом, быстрое вытеснение плотной жидкости газом неустойчиво. В анизотропной среде скорость роста возмущений малой амплитуды зависит от направления теченияэксперимент и расчет демонстрируют качественное совпадение. Отметим, что два крайних случая — диагональная и продольная ориентации анизотропной среды — дают верхнюю и нижнюю границы для инкремента возмущения в изотропной среде. Развитие неу-стойчивостей приводит к формированию картины струй, форма которых определяется симметрией среды. Фрактальные структуры, типичные для линейного режима фильтрации, не наблюдались в квадратных решетках каналов. В самом общем виде можно связать это с инерционной природой сопротивления при высокоскоростной фильтрации: на мезоуровне преобладает сопротивление формы, а не трения. Инерция мезо-течения (струек в порах) не благоприятствует разворотам, необходимым для разветвления. Тем не менее вытеснение резко неоднородно, и выступы заметно опережают средний фронт. Следует ожидать, что в изотропном течении «невязкие пальцы» могут быть выражены слабее, чем в продольном вытеснении, но во всяком случае сильнее, чем в диагональном.

Предложенный закон (1.43) хорошо иллюстрирует природу квадратичной «инерционной» составляющей и для сопротивления в изотропной пористой среде. Течение в порах можно представить как набор струек со средней скоростью V. Случайная поперечная компонента скорости в общем случае также порядка V. Поперечный массообмен между струйками транспортирует часть жидкости из «приблизительно продольных» каналов в поперечные, где и тормозится направленное течение. Такие потерипродольного импульса могут быть одним из основных механизмов, порождающих квадратичное сопротивление.

Соответствие между стохастическими расчетами в изотропной модели и экспериментом лучше для слоя шариков. В почти горизонтальном слое получены ветвящиеся структуры, хотя и более плотные, чем расчетные кластеры в случае изотропного режима вытеснения. Отметим, что в трехмерной регулярной засыпке практически нет прямых каналов, и ее анизотропия для быстрого «инерционного» режима фильтрации заметно меньше. Поэтому следует ожидать, что реальное трехмерное вытеснение будет ближе к фрактальному режиму. В случайной пористой среде из-за местных неоднородностей присутствие каналов возможно, и течение приобретет некоторые черты канального вытеснения.

В реагирующей среде, как в процессе конвективного горения, на эти процессы может влиять добавочное газовыделение. Если струи горячего газа будут инициировать реакцию, это может как усиливать неустойчивость (из-за роста давления), так и сглаживать ее (из-за роста местного сопротивления). Для нас важно, что начальные стадии неустойчивости во всяком случае имеют место и в реагирующей среде.

1.4 Дискретные подходы к двухфазным течениямПри моделировании физического процесса важно выбрать правильный подход. Для двухфазных задач, в которых существенна структура среды и случайная компонента течения, а также неустойчивости, заслуживают внимания модели, отражающие такие свойства. В данном разделе рассмотрим метод клеточных автоматов (cellular automata), называемый также методом решеточного газа (lattice gas) — сравнительно новый подход к моделированию газодинамических процессов. Движение сплошной среды моделируется прямым образом расчетом эволюции специального «микромира» частиц, которые перемещаются по фиксированной плоской или пространственной решетке и сталкиваются между собой в узлах решетки. По существу, речь идет о предельно упрощенном варианте молекулярнойдинамики.

Решеточная газодинамикаЧисленное моделирование газодинамики развивается с 40-х годов и основано главным образом на применении разностных схем. Такие традиционные методы, называемые в дальнейшем детерминированными, с успехом применялись при решении многих проблем. Вместе с тем детерминированным методам присущи недостатки, затрудняющие их применение в ряде приложений.

Один из примеров такого рода — неустойчивые задачи. Разностные схемы сами по себе могут проявлять паразитную «счетную» неустойчивость. Если же их применять в задачах с неустойчивостью физического происхождения, то нередко удается описать лишь начальные стадии ее развития. Продвижение в нелинейную область связано с потерей точности.

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

Полное гидродинамическое описание процессов в каждой индивидуальной поре и невозможно, и не имеет особого смысла. С другой стороны, полное «размазывание», когда от всей геометрии остается только один параметр — пористость — недостаточно продуктивно. Метод клеточных автоматов (далее КА) представляет собой пример компромисса, позволяющего качественно учитывать флуктуации и применяемого во многих задачах газодинамики. Этот метод описывает движение на микроскопическом уровне, но с минимальной подробностью.

Модель клеточных автоматовОпишем простейшую стандартную модель, впервые изложенную в работеФриша, Хасслахера и Помо (FHP) [43]. Рассматривается газ, состоящий из одинаковых «атомов» — частиц, способных двигаться по жестко заданной плоской решетке. За единицу времени атом перемещается на единичное расстояние — из одного узла решетки в соседний. В узлах атомы могут взаимодействовать, рассеиваясь друг на друге. Так происходит передача импульса. Полный цикл расчета состоит из перемещения атомов по связям между узлами и рассеяния в узлах.

На первый взгляд, естественно рассматривать квадратную решетку. Но она имеет существенный недостаток — отсутствие в среднем вращательной изотропии. По-видимому, по этой причине более ранние работы [44] с квадратной решеткой не привлекли заметного интереса.

Как показано в [43], более разумный вари- /УУуУУЛ /ант — покрытие плоскости равносторонними VVVVA/ /треугольниками, так что из одного узла выходддддд —=> дит шесть связей (рис. 1.26). В каждом узле / может находиться до шести атомов с различлными направлениями скорости. Более одно-^го атома с данной скоростью не допускается (принцип Паули). Возможно наличие одного РисI-26неподвижного атома. На рис. 1.26 показаны также правила рассеяния для конфигураций, которые изменяются при взаимодействии в узлах. Остальные конфигурации не взаимодействуют. Разумеется, аналогичные правила применяются для конфигураций, повернутых на угол, кратный 7г/3. При лобовом столкновении результат выбирается из возможных случайным образом. Так в расчет вносится стохастичность. На первый взгляд, кажется странным, что не взаимодействуют встречные частицы на связях. Однако их упругое столкновение приводило бы к тому же результату — обращению скоростей — как и пролет без взаимодействия.

Правила, определяемые рис. 1.26, не единственно возможные. Вообще, допустимы любые конечные состояния при условии сохранения числа частиц и импульса при рассеянии. От набора взаимодействующих конфигураций зависит скорость обмена информацией между подсистемами частиц, имеющими разные направления скорости. Рассеяние «зацепляет» между собой эти подсистемы и таким образом приводит к установлению локального равновесия.

Авторы [43] предназначали модель для имитации практически несжимаемых течений. В этом случае первое уравнение не нужно, а во втором при постоянной плотности возможно масштабирование скорости: v = д (р)и. Это улучшает форму уравнения движения, хотя нефизическая «скоростная» составляющая давления остается.

Используя разложение (1.49) для nJ имеемтт и2 с чHifc = -jOik + —{щик — ywfcj > откуда и следует второе уравнение системы (1.48). Заметим, что j не зависит от правил рассеяния (т. е. другие правила и возможные процессы не изменяют полученное значение [45]). В следующих приближениях по скорости вращательная инвариантность нарушается.

Легко видеть, что это состояние равновесно. На квадратной решетке рассеяние происходит при лобовых столкновениях, причем состояния 1,3 переходят в 2,4 и наоборот. Но такие процессы из симметрии не меняют заданного распределения по направлениям.

Компоненты средней скорости их = иу = v, поэтому динамическая часть потока импульса обращается в нуль. В этом можно убедиться, рассмотрев поток импульса, скажем, в направлении X. Частицы типа 1 уносят положительный импульс п1 через единичную площадку, а частицы типа 3 приносят отрицательный импульс (—п3). В результате полное изменение импульса пропорционально п1 + п3 = р/2 за единицу времени. Действуеттолько давление р = р/2, а от скорости поток импульса вообще не зависит. Это и показывает неудовлетворительность квадратной решетки.

Подобно тому, как в «настоящем» классическом газе после небольшого числа соударений устанавливается локально равновесное распределение по энергиям, в решеточном газе вырабатывается распределение типа Ферми — Дирака (поскольку возможны заполнения состояний только 0 и 1) [43, 45]. В вероятности столкновений входят произведения множителей вида /J = nJ/(l — nJ). В равновесии такие произведения должны быть равны для начальных и конечных состояний. Это выполняется, когда логарифм рявляется линейной функцией сохраняющихся при столкновении величин, т. е. числа частиц и импульса. ОтсюдаИз этого распределения также можно получить уравнения (1.48) разложением по степеням скорости.

В трех измерениях не существует подходящей решетки, обеспечивающей одновременно заполнение пространства и вращательную инвариантность уравнений даже в квадратичном по скорости приближении. Однако изотропность достигается для трехмерной проекции четырехмерной гране-центрированной решетки (face-centered hypercubic, FCHC [46]).

Сравнительная характеристика модели КАОчевидна чрезвычайная простота расчетного алгоритма. Исключительно экономно расходуется память компьютера: в простых вариантах, аналогичных рассмотренному выше, состояние узла занимает менее одного байта. Велика скорость вычислений, так как все операции проводятся с небольшими целыми числами. Легко достигается распараллеливание задачи. Могут применяться процессоры специальной архитектуры [47].

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

Недостатки подхода являются «продолжением» достоинств. Сравнительная грубость модели приводит к -^сшу,.чта она аппроксимирует не точные уравнения газодинамики, а несколько искаженную систему (1.48). Не только фактор д{р), но и нефизическое давление, пропорциональное квадрату скорости, проявляются в расчетах [48]. В простейшей модели, по существу, отсутствует понятие температуры.

Дискретность скачков частиц приводит к появлению нефизических интегралов движения (staggered momenta), связанных с возможностью разбиения среды на две подрешетки [49]. Например, легко убедиться, что на квадратной решетке разность у — импульсов в точках с четными и нечетными у — координатами на каждом шаге меняет знак. Следовательно, эта величина, умноженная на (—1)', сохраняется. При специально подобранных начальных условиях возможно влияние таких инвариантов на гидродинамические моды. В несколько усложненных вариантах метода (скажем, при допущении хотя бы двух величин скорости частиц) эти недостатки значительно ослаблены. Вязкость накладывает ограничения на число Рей-нольдса Re достигаемое в реальных расчетах (см. ниже).

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

Разумеется, результаты подхода будут в основном качественными. Но это же относится к гораздо более рафинированной молекулярной динамике Леннарда — Джонса и во многих случаях к детерминированным разностным методам. Как указывают авторы [43], часто для получения результатов весьма скромной точности ведутся масштабные промежуточные вычисления с большим числом знаков. Резюмируя, можно сказать, что клеточныеавтоматы имеют свою область применимости, в которой традиционные методы мало полезны.

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

Первые применения модели касались несжимаемых течений в каналах и задач двумерного обтекания. Отметим расчеты трехмерного обтекания плоского диска, расположенного поперек течения [46]. При Re = 190 на сетке 128×128×256 наблюдалось возникновение вихрей, срывающихся с кромки, и развитие существенно трехмерного течения. Картина близко согласуется с классическими результатами экспериментов.

В работе [50] авторы рассмотрели две несмешивающихся решеточных жидкости, которые можно назвать «красной» и «синей». Чтобы поверхность раздела между ними не размывалась со временем, введено поверхностное натяжение. Оно моделируется модификацией правил рассеяния так, что частицы данного цвета посылаются в направлении, где этот цвет преобладает, при сохранении импульса и числа частиц. Кроме того, предложен способ восстановления галилеевой инвариантности введением большого числа М покоящихся частиц, что приближает фактор д (р) к единице. К сожалению, оптимальное число М зависит от плотности, что ограничивает этот рецепт несжимаемыми течениями. Учитывалась также плавучесть: с некоторой вероятностью частицы одного цвета приобретали импульс вниз, а другого — вверх. В результате в [50] решена задача о неустойчивости Релея — Тейлора. Начиная с неравновесного состояния (тяжелая жидкость выше легкой), прослежено развитие неустойчивости с образованием характерных грибообразных струй, распадом их на капли и, наконец, достижением устойчивой стратификации.

В [51] моделируется фазовый переход жидкость — газ. Жидкая фаза в основном состоит из связанных пар частиц, эквивалентных одной покоящейся частице двойной массы. В это состояние частицы газа могут перейтипри лобовом соударении. Связанные пары в соседних узлах решетки взаимодействуют между собой с потенциальной энергией, равной —1. Распад пары происходит с вероятностью, зависящей от энергии АЕ ее взаимодействия с окружением, и от температуры. Температура задается независимо и входит через функцию ехр (—АЕ/Т) в вероятность. Изотермы в описанной среде демонстрируют критическое состояние и двухфазную область. Зародыши жидкой фазы в области сосуществования сливаются в капли, а затем образуют сплошную жидкость. В [52] рассмотрена трехфазная система.

Различные течения со свободными границами рассмотрены в [53]. Здесь также введен цвет частиц и изучен ряд модификаций модели. Прежде всего моделировалась неустойчивость Кельвина — Гельмгольца в сдвиговом течении. Прослеживается развитие гидродинамических возмущений, усиливающее взаимную диффузию.

Чтобы иметь возможность представлять границу раздела несмешиваю-щихся жидкостей, необходимо подавить диффузию. Для этого в [53] применен «химический» подход. Тройные соударения с участием разноцветных частиц приводили к реакции изменения цвета. В результате возникает тонкий переходной слой между жидкостями, так как частица, углубившаяся в чужую область, меняет цвет. Моделирование на решетке 1024×256 точек при Re «330 демонстрирует развитие типичных волновых и нелинейных вихревых структур. Исследовалась также релей — тейлоровская неустойчивость. В общем, результаты работ [50, 51, 52, 53] иллюстрируют перспективность КА в неустойчивых проблемах.

Заметная часть работ посвящена применению КА к течению в пористой среде. В [54] изложены результаты масштабных двухи трехмерных вычислений. На моделях каналов воспроизводилось течение Пуазейля. Отсюда определялась вязкость решеточного газа. На твердых границах задавалось отсутствие проскальзывания, т. е. частица, попавшая в граничный узел, рассеивалась назад. В регулярных пористых средах хорошо подтверждается закон трения Дарси. В виртуальной среде с фрактальными свойствами изучена зависимость проницаемости от пористости и фрактальнойразмерности. Моделировался также реальный эксперимент — задавалось та же геометрия, как в реальном пористом образце, полученном травлением стекла. Вычисленная и экспериментально измеренная проницаемости совпали с точностью 10%.

В [55] рассмотрена неустойчивость Саффмана — Тейлора, возникающая при вытеснении из порового пространства вязкой жидкости. Вводилось случайное распределение узлов — рассеивателей, размещенных с небольшой плотностью в расчетной области. Рассеиватели моделировали трение, заменяя скорости частиц на обратные. В невязкой жидкости можно задавать отсутствие рассеивателей, но в [55] вместо этого вводилась компенсирующая трение массовая сила, что уменьшает градиенты плотности. Граница раздела задавалась явным образом как «цепь» узлов, способная деформироваться в результате ударов частиц каждой жидкости. Моделирование качественно воспроизводит известные экспериментальные результаты, включая разветвление пальцев при малом поверхностном натяжении.

В [56] моделировались диффузионные и химические процессы, происходящие при фильтрационном течении в пористой среде. Кроме скоростей частиц, в каждом узле задавалась концентрация примеси. Эту концентрацию частицы переносят в другие узлы, где после каждого перемещения результат получается усреднением. Реакции на поверхности пор изменяли концентрацию, переносимую отскочившими частицами. После некоторого числа таких отскоков граничный узел становился внутренним, т. е. увеличивался объем поры. Модель демонстрирует течение в довольно сложной области, сопровождающееся растворением стенок.

В [53] моделировался фронт горения в газе. Для представления исходной смеси и продуктов горения введены два типа частиц. Несгоревшие частицы имеют двойную массу и единичную скорость. Продукты горения состоят из частиц единичной массы с двойной скоростью. Следовательно, величины импульсов одинаковы, что позволяет пользоваться теми же правилами соударений. Температуры таких газов отличаются в 2, а массовые плотности — в 4 раза (при равенстве давлений).

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

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

Волны в решеточном газеУдарные волны с резким изменением плотности не рассматриваются обычно как подходящая область применимости КА. Это мнение имеет определенные основания, вытекающие из формы континуального аналога (1.48). Приводимые ниже результаты, хотя и достаточно просты, представляют некоторый методический интерес.

На рис. 1.27 видно различие между классической и решеточной газодинамикой. Допустимы значения скорости волны в довольно ограниченном интервале, для условий рис. 1.27: 1 < D/c < 1,034. Это связано с обращением д в нуль при р «3. Вблизи этогозначения имеем приблизительно акустику, а при произвольной плотности — небольшие отклонения от звуковой скорости. Заметим, что имеется два пересечения адиабаты и прямой, т. е. возможны два скачка с одной и той же скоростью.

Кривая 2 на рис. 1.27 соответствует начальной плотности ро = 5, большей значения, при котором д (р) обращается в нуль. Сверхзвуковыми здесь будут волны разрежения.

Описанные явно нефизические режимы наблюдаются в численном эксперименте. На рис. 1.28 показаны расчеты распада разрыва при малой и большой плотности в секции высокого давления. При соотношении плотностей 3/1 (рис. 1.28,а) течение имеет в общем обычный газодинамический вид с ударной волной, идущей вправо, и плавной волной разрежения, распространяющейся влево. Метка на оси х — начальное положение разрыва. Расчет проводился на сетке 100×181 без усреднения по ансамблю. Совсем другая картина (см. рис. 1.28 б) получается при соотношении плотностей 6/1 (поскольку использована 7-битовая модель с покоящейся частицей, плотность 6 еще не является предельной). Здесь влево идет лидирующая ударная волна разрежения. От нее отстает (из-за сноса течением) вторая подобная волна. Видно, что волновые движения в плотных решеточных газах отклоняются от газодинамики даже качественно. Изложенное касалось наиболее ранней, простейшей модификации метода. Несколько лучшие результаты получены в [51], где фактор д остается приблизительно посто* = 40 t = 30х хРис. 1.28.янным в заметном диапазоне плотностей. Достигнута качественно разумная форма ударной адиабаты. Характер волновых движений может быть заметно приближен к реальности некоторым усложнением модели. Уже введение двух скоростей значительно уменьшает описанные неприятности [57, 58].

Перспективы КА в физике горения и взрываПростейшие варианты решеточных методов типа FHP можно применять к реагирующим течениям в пористой среде. Выделенность системы отсчета, связанной с пористым скелетом, до некоторой степени ослабляет роль нарушения галилеевской инвариантности. Скорость течения из-за трения невелика, так что плохое описание динамической части потока импульса также не критично. Легко внести, например в духе [56], рождение частиц, описывающее горение на поверхности. Для учета изменения температуры необходима модель хотя бы с несколькими скоростями частиц и уровнями энергии [57, 58, 59]. Однако в некоторых постановках, например в задаче о детонации в пористой среде из гранул унитарного топлива, горящих с поверхности, изменения температуры малосущественны. Следует отметить, что неодномерность фронта реакции и распространение его в виде «языков» получаются сами собой. Учет случайности характеристик среды также происходит автоматически.

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

Сейчас неформальной оценкой состояния вопроса является то, что работы, касающиеся КА, в основном сосредоточены в физических журналах. Пока что метод имеет довольно ограниченное применение среди механиков и не рассматривается как стандартный. Отметим, что традиционная вычислительная гидродинамика развивается уже более 50 лет. Решеточные же методы пока находятся в начале своего развития. Можно ожидать, что модификации подхода существенно расширят круг решаемых задач.

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

Применения КА к задачам данной работы изложены в следующей главе.

Глава 2ДВУХФАЗНАЯ ДЕТОНАЦИЯВозможно, одна из самых интересных задач физики и механики двухфазных сред — это проблема распространения волн реакции в активных системах. Известно множество качественно различных режимов, многие из которых изучены недостаточно. Сравнительно давно исследуется горение пористых энергетических материалов, но до сих пор некоторые процессы, в особенности переход горения в детонацию, не до конца ясны. Детонация аэровзвесей понятна гораздо в большей степени1, чем детонационно-подобные режимы в, казалось бы, сходных системах типа «газ-пленка» в пористой среде. В данной главе рассмотрены вопросы распространения и инициирования быстрых квазидетонационных волн в пористых активных средах. Именно в таких системах наиболее важны мезопроцессы на уровне пор.

Сразу следует оговорить терминологию. Обычно принято называть детонацией самораспространяющийся процесс, который ведется ударной волной, в согласии с моделью Зельдовича — Неймана — Деринга. На определенном расстоянии от ударного фронта происходит переход от дозвукового к сверхзвуковому течению (в системе фронта). Сверхзвуковое втекание вещества во фронт и сверхзвуковой отток от плоскости Чепмена — Жуге обеспечивают стационарность волны.

Однако в последнее время наблюдались процессы, в которых распрост1 Анализ этого круга вопросов проведен, например, в работах С. А. Ждана [60], А. А. Борисова, Б. Е. Гельфанда и др. [61, 62], Р. И. Нигматулнна, П. Б. Вайнштенна и др. [63], В. М. Фомина, А. В. Федорова и др. [64, 65]раняющим агентом заведомо не является ударная волна. При этом сохраняются такие свойства «нормальной» детонации, как стационарность (или квазистационарность) и нечувствительность к внешним воздействиям. Хорошим примером является «вакуумная детонация» [66, 67]. Вопрос о том, привязывать ли понятие детонации только к ударноволновым режимам (считать основным способ передачи возбуждения) либо к более широкому кругу явлений (считать основными свойства быстрого самораспространения и устойчивости). Решение этого вопроса — дело соглашения, которое еще в литературе не установилось. Можно называть неклассические режимы квазидетонацией, чтобы отличить их от обычных ударноволновых. Но, поскольку авторы экспериментальных работ, цитируемых ниже, склоняются к термину «детонация» для волн, скорость которых превышает скорость распространения возмущений в исходной среде, и ради краткости, мы будем часто употреблять этот термин в расширительном смысле.

2.1 Конвективные волны — экспериментальные данныеКак известно, при поджигании пористого взрывчатого вещества или пороха вначале типично возникновение медленного послойного горения. В достаточно прочной оболочке послойный режим переходит в более быстрое конвективное горение. Ускорение нестационарного волнового комплекса и начало реакции в твердой фазе впереди фронта пламени принято считать основными факторами развития детонации [68, 69, 70], хотя в некоторых ситуациях наблюдается переход другого типа, с генерацией в горящем ВВ вторичной волны (см., например, [71]). Стадии нормальной детонации обычно предшествует низкоскоростная детонация. Время перехода может 1 составлять несколько миллисекунд. Естественно, подобные медленные процессы происходят в условиях, обеспечивающих удержание (например, в толстых стальных трубах). При малой прочности оболочки ее разрушение обрывает этот классический сценарий перехода горения в детонацию (ПГД).

2.1.1 Двухфазные режимы — литературные данныеОписанный классический сценарий ПГД — не единственно возможный. Сейчас известны качественно другие процессы. Опишем накопившиеся на данный момент экспериментальные наблюдения.

Низкоскоростная детонация пористого ВВВ работах JI.A. Лукьянчикова, В. В. Андреева и др. показано, что если иницирующий импульс достаточно резок, удается «обойти» медленную стадию послойного горения, и детонация развивается гораздо быстрее (например, в масштабах 10 мкс/10 мм). Применялись такие способы возбуждения, как электрический разряд или взрыв проводника внутри порошкового ВВ [72, 73], впрыск горячих продуктов горения из отдельной камеры при прорыве мембраны [74]. Возбуждалась волна с характерной скоростью около 1 км/с и уровнем давления около 2 кбар. При малом диаметре заряда (менее 4 мм) и при очень слабой оболочке (оргстекло толщиной 1−3 мм) эта волна довольно стабильна2. Схематическая картина течения и профиль давления показаны на рис. 2.1. Особенности указанного волнового режима — это сложный фронт волны, представляющий собой случайный пульсирующий рельеф выступов и впадин, и плавный рост давления, со временем нарастания 1−2 мкс. По уровню давления можно оценить долю сгоревшего вещества — несколько процентов.

Волна сжатия той же крутизны и амплитуды не инициирует процесс [73]. Отметим, что перечисленные методы возбуждения — существенно двухфазные: во всех случаях поток горячих газов резко внедряется в порошок. Скорость инициирующего потока должна быть достаточно велика: > 500 м/с. Мембранное инициирование может произойти при давлении в камере всего порядка 100 атм.

Двухфазность резко меняет условия распространения волны. Хотя из-за малой прочности ВВ некоторая деформация пористого ВВ имеет место, само по себе сжатие недостаточно для возбуждения реакции. Более того, даже при максимальном сжатии (около 25%) сохраняется заметная про2При диаметре 4 мм и выше происходит ускорение процесса вплоть до развития нормальной детонации. ницаемость скелета. Согласно [75], при этом давление в твердой фазе не превышает сотен атмосфер, то есть преобладающий вклад в измеренное давление волны вносит газовая фаза. Совокупность наблюдательных данных говорит о том, что механизм распространения волн — конвективный, или струйный. Воспламенение впереди «среднего» фронта обеспечивают струи горячих газов, прорывающиеся из зоны горения.

Низкоскоростные режимы возбуждались также воздействием на порошок «ударом» газовой детонации [76]. Газовая детонация интересна как инструмент, позволяющий «тонкую» регулировку инициирующего воздействия. Недавнее подробное исследование [77, 78] дало качественные результаты, подобные уже описанным, и выявило новые свойства конвективных волн. В частности, удалось выделить процессы быстрого конвективного горения и низкоскоростной детонации, которые в более грубых ранних опыi тах смешивались. (Это было одной из причин определенной размытостиIтерминологии. Процесс со скоростью порядка 1 км/с назывался как скоростным конвективным горением, так и низкоскоростной двухфазной детонацией [79]).

В [80] исследовались смеси гексогена с мипорой с начальными плотностями ро = 0,2 -г-1,0 г/см3. При минимальной плотности получена скорость детонации 2,4 км/с, при давлении около 5 кбар. Авторы [80] также считают, что в таких режимах реакция протекает в виде горения зерен ВВ с поверхности, а реакция распространяется фильтрацией горячих продуктов между зерен. Вещество переходит из исходного состояния в состояние, соответствующее окончанию реакции, не через скачок параметров, а в результате постепенного их нарастания, связанного с последовательным выгоранием ВВ.

Детонация в жесткой пористой средеАналогичные выводы сделаны в работах В. В. Митрофанова, А. В. Пинаева, Г. А. Лямина и др. где получены своеобразные детонационно-подобные течения в жестких пористых средах, со скоростями волн порядка 1 км/с. Здесь также наблюдался плавный рост давления и струйные образования. Активной компонентой среды может быть газовая смесь, заполняющая поры [81] или слой взрывчатого вещества на поверхности пор [66]. В [82, 83, 84, 85] исследована более сложная система «газ-пленка», в которой газообразный окислитель в порах реагирует с пленкой топлива на поверхности поровой структуры. Прочность пористого скелета (например, засыпки стальных шариков) совершенно исключает инициирование смеси гомогенным сжатием среды. Из-за потерь на трение скорость волны недостаточна и для возбуждения газовой фазы за фронтом лидирующей ударной волны (как это происходит при детонации аэровзвесей) — отметим, что самоподдерживающаяся волна возможна и в вакуумированной системе («вакуумная детонация ВВ» [66, 86]). В этих последних работах концентрация ВВ составляла от единиц до десятков мг/см3, что несколько меньше концентрации выгоревшего вещества в [72, 73, 79].

2.1.2 Детонация сильноразбавленного ВВВ дополнение к уже описанным экспериментальным данным приведем наши результаты, полученные на несколько другой системе. Исследовалась детонация смесей тэна с кварцевым песком при содержании ВВ 10 — 30%. Здесь низкоскоростной режим стабилизируется разбавлением ВВ. Значительно большая прочность наполнителя позволяет надеяться на получение существенно двухфазных режимов детонации. Кроме того, уменьшение содержания ВВ приводит к снижению давления и скорости детонации и, следовательно, бризантного воздействия. Такие смеси могут применяться в случаях, когда требуется сравнительно «мягкое» воздействие на материал.

В работах [87, 88, 89, 90, 91], посвященных исследованию детонации при сильном разбавлении, обнаружены в основном стандартные режимы с ли- 1 дирующей ударной волной (УВ) и более или менее выраженным химпиком. Отметим, впрочем, что в [89] отмечается размытие ударного фронта при: увеличении размера частиц добавкив зарядах низкоплотной смеси вспененного полистирола, аммиачной селитры и алюминия с гранулами около 1 см также обнаружен плавный фронт без химпика [92]. Ниже показано, что в смеси тэн-песок вблизи предельного разбавления возможны струйные режимы.

Методика экспериментаРазмер зерна песка составлял S = 90-f700 мкм. Песок промывался в ацетоне, высушивался и просеивался. Использовался мелкодисперсный тэн, получаемый из раствора в ацетоне добавлением воды, от чего растворимость тэна резко падает [93]. Мелкодисперсный тэн представлял собой игольчатые кристаллы длиной 50 и толщиной 10 мкм. Система нужной концентрации готовилась тщательным перемешиванием ВВ и наполнителя. Связь частиц ВВ с зернами песка была достаточно устойчивой и допускала насыпное заряжание без заметного расслаивания. При содержании ВВ (3 = 20 — 25% фоторазвертки детонации не выявили колебаний скорости фронта, а съемка выхода детонации на торец заряда демонстрирует однородный (с точностью до размера пор) фронт свечения, радиус кривизны которого около 2 диаметров заряда.

Детонация возбуждалась электрическим разрядом. Скорость детонации D измерялась оптическим методом, а также электроконтактным вместе с массовой скоростью. Результаты измерений усреднялись по 2 — 3 опытам. Изучалось влияние размера зерен песка, диаметра заряда d и процентного содержания ВВ (3. В опытах d = 4 -г- 32 мм. Применялись оболочки из оргстекла, стекла и винипласта. Заметного влияния материала оболочки в опытах не отмечено. Далее размеры оболочки обозначаются в виде дроби: внешний диаметр/внутренний диаметр (в миллиметрах).

Схема опыта приведена на рис. 2.2. В качестве датчика массовой скорости частиц наполнителя (песка) v использовалась рамка из алюминиевой фольги шириной 1 мм и толщиной 0,07 мм, к концам которой приваривались медные выводы. В двухфазной системе возможно относительное движение твердой и газовой фаз, и плечо рамки движется, вообще говоря, со скоростью, промежуточной между v и скоростью газа и. Однако из-за малого содержания ВВ рнс. 2.2. эта средняя скорость ближе к v. Давление в условиях наших опытов практически определяется скоростью плеча датчика. Оценка влияния различия скоростей приведена ниже.

Для измерения и применялся разомкнутый датчик: два отрезка медной проволоки диаметром 0,5 мм, параллельных оси заряда. Базовое расстояние между ними составляло 0,5 см. Сопротивление продуктов детонации между этими контактами R & 250 Ом. Поскольку это сопротивление не постоянно, сигнал снимался в режиме холостого хода через большое балластное сопротивление, что уменьшало амплитуду более чем на порядок. В частности, по этой причине и потребовалось магнитное поле большой напряженности.

В момент прихода детонационной волны на разомкнутые электроды датчика появлялся сигнал даже при отсутствии магнитного поля (по-видимому, из-за влияния инициирующего разряда). От этого сигнала удалось избавиться, передавая детонацию через пластину из оргстекла 3 мм), изолирующую цепь подрыва. Установлено, что при не слишком малых разбавлениях изолирующая прокладка заметно не влияет на процесс, так как величины v и D существенно не изменялись. О том же говорят данные фоторегистрации, описанные ниже. При (3 10% детонация не передавалась, поэтому содержание ВВ в инициирующей части повышалось до 20 — 25%.

Ниже или выше по течению разомкнутого датчика ставилась алюминиевая рамка, что позволило в одном эксперименте определять массовые скорости газа и скелета, а также D (по интервалу между сигналами). По сигналу с алюминиевой рамки оценивалась ошибка, вносимая производной dH/dt в момент измерений. Предварительно снимался сигнал с датчика и замкнутого проволочкой, для определения его эффективной площади. Задержка подрыва подбиралась так, чтобы при dH/dt — 0 начинал работать датчик скорости газа. В целом ошибка из-за dH/dt мала.

Проведены контрольные опыты с чистыми тэном и гексогеном насыпной плотности в зарядах 21,5/18. Получено удовлетворительное совпадение с известными данными. Например, для тэна получено давление р =68,8 кбар (64,5 кбар в [95]). Пиковое значение сигнала с П-образного датчика было на 20 — 25% выше, чем с разомкнутого датчика, так как последний усредняет скорость вещества в проводящей зоне. Средняя скорость, естественно, ниже максимальной в детонационной волне с химпиком.

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

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

Р РПоскольку рс<�р р, w v даже при заметной разнице скоростей фаз (Z w).

Преобразуя (2.2) к переменным w и Z, имеемPGVPS®гу2 /о ор-ро = pHDw—Z. (2.3)РЗдесь рн = Posao + Ро’Ро начальная плотность смеси. Выражение (2.3) заменяет обычную связь давления и скорости. Второе слагаемое в (2.3) представляет собой поправку, относительный вклад которой порядка [pG.

Таблица 2.1.

При большом содержании ВВ двухфазность системы не сказывается в качественном смысле. Так, в опытах по измерению массовой скорости газа и для оболочек 21,5/18 при /3 = 25% и = 400 ± 100 м/с. Измерение скорости частиц песка П-образным датчиком дает v = 500 ± 90 м/с. Следовательно, можно утверждать, что скорость газа во всяком случае не рИС- 2.3. превышает скорости частиц: (и < v), и при таком содержании ВВ процесс практически односкоростной. Реакция идет обычным образом за ударной волной. Фронт импульса v крутой, наблюдается некоторое подобие хим-пика (рис. 2.4,а).

0,1 0,2 0,3ач.

Рнс. 2.5. Передача детонации через преграду из оргстекла. Заряд 21,5/18, /3 = 25%. а) — сплошная пластина, б) — перфорированная (18 отверстий).

Другой характер имеет волна при низком содержании ВВ (менее 15%). Наблюдается плавное нарастание скорости частиц за время 1 мкс (рис. 2.4,6). Скорость газа и несколько больше v, что говорит в пользу распространения детонации не УВ, а конвективным образом, т. е. потоком газов, постепенно увлекающим зерна песка. Низкая чувствительность смеси не позволила провести сопоставление по механизму передачи детонации. Отметим, что двухфазные эффекты проявляются при содержании ВВ, близком к критическомувоспроизводимость таких опытов была заметно хуже.

Граничное содержание тэна можно оценить из следующих соображений. При плотной упаковке частиц песка (пористость <р = 0,35) и малом содержании ВВ более мелкие зерна ВВ расположены в порах и заполняют их при массовом содержанииРе'-Ро/(рефО + POS^O) 0,2.

Здесь ре = 0,9 г/см3 — насыпная плотность ВВ, pos^o 1,3 г/см3 — насыпная плотность песка, <ро = 0,4 — начальная пористость песчаной засыпки.

При такой концентрации ВВ между зернами песка непременно имеются прослойки тэна (а реально они появляются при несколько меньшем содержании). Сжатие системы УВ приводит к деформации ВВ в прослойках, что может повлечь развитие реакции. Если же содержание ВВ мало, то зерна наполнителя контактируют между собой и защищают ВВ от воздействия волныдетонация передается конвективным способом. Практически граница между двумя режимами размывается из-за флуктуаций засыпки частиц и их возможного дробления в волне. Данные табл. 2.1 (немонотонное изменение плотности смеси с добавлением ВВ) в общем подтверждают эту оценку. Максимум плотности смеси, соответствующий заполнению пор, достигается между 12,5 и 15%, что указывает на смену типа упаковки смеси в этом интервале. Лучшее совпадение с приведенной выше оценкой получится при более разреженном состоянии тэна: ре < 0,9 г/см3, чего и следует ожидать при укладке ВВ внутри мелких пор песка.

Таким образом, получены данные, указывающие на возможность различных режимов детонации в смеси тэн — кварцевый песок. Двухскорост-ной режим наблюдается при малом содержании ВВ, недостаточном для заполнения пор. При большем содержании ВВ процесс односкоростной и ведется ударной волной."Взрывчатый песок" занимает промежуточное положение между чистым порошкообразным ВВ [72, 73] и ВВ в жесткой прочной структуре [66, 86]. Таким образом, во всем спектре таких пористых систем можно возбудить низкоскоростные волны с конвективным механизмом распространения. Для них существенным условием является незначительность сжатия вещества, с сохранением проницаемости порового пространства.

2.2 Континуальная модель двухфазного теченияВ этом разделе кратко рассмотрим применяемые в литературе принципы моделирования конвективных волн в пористой среде и изложим вариант континуальной модели, использовавшийся в наших расчетах.

Начиная с работ Куо, Саммерфилда и др. [96], выполнено многорасчетно-теоретических исследований конвективных волн как этапа перехода горения в детонацию. Наиболее изучен случай развития процесса в прочных оболочках, близкий к одномерному. Характерные масштабы явления: время 100 -г 1000 мкс, длина 50 500 мм (в зависимости от типа пороха или ВВ и размера частиц). Горение начинается в районе воспламенителя. Горячие газы проникают в пористую среду, нагревая зерна ВВ. Через некоторое время зерна начинают гореть, также производя горячие газы. Пористая среда сжимается в результате действия потока газа. Между экспериментальными и расчетными работами достигнуто некоторое согласие. Следует упомянуть зарубежные работы Криера [97, 98], Гофа и Шварца [99], Баера и Нунциато [100], Прайса и др. [101]. Среди отечественных исследований большой интерес представляют работы Института химической физики (А.А. Сулимов, А. А. Борисов, Б. С. Ермолаев и др.) [102, 103, 104] и Института механики МГУ (Р.И. Нигматулин, И. Ш. Ахатов, П.Б. Вайн-штейн) [25, 26].

Для больших плотностей совпадение хуже. Отметим, что в эксперименте деформация превосходит <ро, что невозможно для несжимаемых зерен. В то же время давление недостаточно для заметного сжатия сплошного материала. По-видимому, экспериментальные данные в области высоких плотностей менее надежны. Зависимость (2.6) предполагает практическую несжимаемость сплошного вещества (давление стремится к бесконечности при <р ipo), что представляется разумной асимптотикой. Для наиболее интересных режимов, описанных ниже, сжатие не достигает области полного закрытия пор, и поэтому неопределенность в давлении скелета несущественна.

Уравнение (2.6) выполняется при ip < ipm{n < <ро • Здесь.

Зная напряжения, легко выразить упругую часть Ер энергии твердой фазы, то есть работу, высвобождаемую при разгрузке. Довольно громоздкое выражение здесь не приводится. Отметим, что остальная часть работы сжатия порошка (для полной разгрузки соответствующая площади петли на рис. 2.6) идет на нагрев твердой фазы. В явном виде тепловой эффект сжатия введен в работах [25, 26], где показана его важная роль в процессе ПГД. В этом отношении данная работа несколько развивает модель [25, 26], поскольку здесь тепловая доля работы сжатия не предписывается заранее, а следует из состояния вещества и истории его деформирования. Межфазная сила, теплои массообменДля силы трения / в большинстве работ используется формула Эргана [105]. По более поздним данным [106] величина силы трения примерно вдвоеменьше. Поэтому в уравнение Эргана вводился поправочный коэффициентгде d — размер частиц пористого скелета, Re^ - число Рейнольдса, определяемое по относительной скорости фаз и диаметру частицы. При скорости обтекания в сотни м/с число Рейнольдса порядка сотен тысяч, так что сопротивление в основном квадратично по скорости.

При достижении некоторого условия (например, заданной температуры твердой фазы) начиналось горение вещества. Считалось, что зерна газифицируются с поверхности: J = psAS ¦ ип (2.9)причем интенсивность массообмена J зависит от параметров среды через ип — скорость регрессии поверхности. Весь поступивший в поры газ сгорает сразу (в правой части уравнения энергии газа (2.5) входит произведение JQ). Реакция в твердой фазе не учитывалась. Обычно в расчетах горения принимают, что начавшаяся реакция выключает теплообмен (тепло, переданное частице от газа, при горении возвращается в газовую фазу). В данной работе также принят этот переход «теплообмен-горение», но, чтобы избежать резких переключений, он размазывался на небольшую область управляющего параметра.

Стандартное выражение для скорости регрессииип = Bp. (2.10)Например, для тэна константа В дает ип = 10 см/с при давлении 1 кбар [108]. Использовались и другие законы горения, которые обсуждаются ниже.

В систему (2.5) входит только средняя температура твердой фазы Ts. Такое приближенное рассмотрение использовано в [25, 26], тогда как в [102, 103] учитывается распределение температуры внутри зерен, и горение начинается при заданной температуре поверхности. Эти подходы не так далеки, что показывает следующая оценка. При постоянном тепловом потоке к поверхности частицы q ее средняя температура возрастает на ATs = 6qt/(dpcp) за время t. Прирост же температуры поверхности можно оценить как ДТ/ = q^/t/(Лрср). ОтсюдаATSS/ATS = d/(6 Vat), где, а = Х/(рср).

Для тэна температуропроводность, а = 8,4 • 104 см2/с [109]. Соотношение прироста средней и поверхностной температур для времен прогрева 1 мкс и 10 мкс при размере зерна 0,3 мм показано в таблице 2.2. В [25, 26] приТаблица 2.2.

ДТа, К 3,5 7,0 14,0ДТ/, К (t = 1 мкс) 600 1200 2400ДТ/, К {t = 10 мкс) 190 385 770нимается ATS = 8,2 К либо 55 К, а в [102, 103] ДТ/ = 300 К. Нас будут интересовать главным образом времена в интервале 1−10 мкс. Тогда значение ATS = 7 К будет разумной оценкой «средней» температуры вспышки, поскольку и поверхностная температура попадает в разумный диапазон.

В пользу усреднения температуры частиц говорят (кроме упрощения модели) два обстоятельства. Во-первых, при дроблении возникают частицы самых различных размеров, и учет распределения температуры в каждой фракции был бы весьма трудной задачей. Во-вторых, поверхностнаятемпература вспышки вряд ли может считаться известной с приемлемой точностью, и ДТ/ на практике приходится назначать из соображений согласования с экспериментом. Как первое приближение, усредненный подход представляется наиболее приемлемым. В части расчетов ДТ5 варьировалось вокруг среднего значения 7 К в обе стороны: до 3,5 К или 14,0 К.

Задача окончательно определяется заданием начальных и граничных условий. Обычно на границах расчетной области ставилось условие непротекания. Наиболее часто используемый в литературе «мягкий» вариант начальных условий — зажигание определенной части заряда при низком давлении. Это приводит к медленному росту давления при разгорании вещества и постепенному развитию волны компактирования порошка. Более соответствует экспериментальным работам [72, 73], которые мы и имели целью описать, условие «резкого» инициирования, когда вначале имеется область горячего газа с достаточно высоким давлением (сотни атм). Ближе всего эта постановка к условиям [74], где такие давления развивались при горении ВВ в предварительной камере, откуда инициирующий поток прорывался в основной заряд после разрушения мембраны. При резком инициировании в порошок внедряется быстрый высокоскоростной поток горячего газа, что значительно ускоряет все существенные процессы (теплообмен, компактирование порошка).

Что касается способа решения уравнений, мы в основном следовали работам [25, 26]. Использовался эйлеров метод Лакса-Вендроффа, модифицированный для учета недивергентных и алгебраических членов, с сохранением второго порядка аппроксимации. Шаг по времени выбирался из условия Куранта, а также условия стабильности по правым частям. Слабая двухфазная неустойчивость подавлялась сглаживанием параметров течения, что равносильно небольшой искусственной вязкости. В гл. 1 вопрос о неустойчивости разобран достаточно подробно.

Таким образом, двухфазная двухскоростная модель (2.5 — 2.9) использует общие балансовые соотношения. Отличительные стороны моделииспользование близкого к эксперименту закона деформирования тэна (вбольшинстве расчетных работ используются гораздо более жесткие уравнения состояния, описывающие трудносжимаемые пороха), возможность учета роста поверхности при дроблении, нетрадиционных механизмов газификации и, наконец, резкое инициирование. Влияние этих факторов демонстрируется ниже при обсуждении конкретных расчетов.

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

2.3.1 Стационарная волнаНапомним основные характеристики стационарного процесса по данным [79]. При скорости около 1 км/с давление достигает 2 кбар, массовая скорость 200 м/с, а время нарастания фронта — около 2 мкс. Максимальная деформация скелета е = AV/Vo по этим данным равна 20%, то есть пористость падает до 0,275.

При таком сжатии давление скелета менее 0,1 кбар. Следовательно, преимущественный вклад в давление вносит газовая фаза. Максимальному давлению в волне соответствует плотность газа р «0,2 г/см3 при температуре 3−4 тыс. К. При пористости 0,275 газ составляет примерно 4 -5% массы вещества, что и дает оценку степени выгорания. При диаметре частиц 0,3 мм и времени горения 2 мкс требуется средняя скорость горения (скорость регрессии поверхности) около 1 м/с. Это на порядок больше нормальной скорости горения при среднем давлении 1 кбар [108].

Требуемую скорость реакции могло бы обеспечить дробление вещества, с уменьшением среднего размера частиц на порядок. Это (для квазистационарной волны) представляется маловероятным из-за довольно слабого сжатия. Другая возможность — ускорение горения из-за обдува частиц высокоскоростным потоком газа (эрозионное горение). Напомним, что существенная скоростная неравновесность — непременная черта описанных экспериментов.

Эрозионное горениеЭрозионное горение (по старой терминологии — раздувание) рассматривается в ряде работ [110, 111, 112, 113, 114]. Ускорение процесса объясняется турбулизацией пограничного слоя и соответствующей интенсификацией теплового потока к горящей поверхности. К сожалению, экспериментальные данные в основном имеются только для скоростей потока < 300 м/с. Расчеты, основанные на теории горения Зельдовича и ее аналогах, дают коэффициент эрозии е, представляющий собой отношение скорости горения ип к ее величине в отсутствие обдува и^. При малой скорости обдува V эрозия незначительна, но при больших V эффект е растет приблизительно пропорционально параметру Вилюнова Vi = yf^pVj psu^ где? -локальный коэффициент сопротивления горящей поверхности. Достигаемые максимальные значения в указанной области параметров е «2 — 3. Не очевидна возможность применения установленных закономерностей к более жестким условиям нашего процесса. Однако, если экстраполировать теоретические оценки к значениям скорости потока до 1 км/с, следует ожидать роста е еще в несколько раз, возможно, до 10.

Неустойчивость слояОпишем механизм горения, который возможен в наших условиях, при экстремальных параметрах потока. На поверхности частицы ВВ, обдуваемой потоком газа, образуется слой испаренного вещества с почти нулевой тангенциальной скоростью. На этом контактном разрыве возникает неустойчивость Кельвина-Гельмгольца с развитием волн на поверхности раздела. В пренебрежении вязкостью наиболее важны возмущения с длиной волны порядка 2irh, где h — толщина испаряющегося слоя. Когда амплитуда волны станет порядка h, произойдет срыв гребней волн. Начнется интенсивное перемешивание испаренного вещества с горячими продуктами реакции. В результате следует ожидать быстрого воспламенения. Горячие газы начнут прогревать следующий слой вещества. Циклическое повторение такогопроцесса обеспечивает большой градиент температуры в тонком поверхностном слое и, следовательно, высокую скорость горения.

Инкремент неустойчивости возмущения с волновым вектором к можно представить в виде [115]kvy/pth (kh)/f/ 7 1 + pth.(kh)/f/ 'jгде р, v — плотность и скорость внешнего потока, р — плотность испаряющегося слоя.

Наиболее важны неустойчивости, у которых kh 1, для которых время ростаг = 1/7 у/p/p'h/v.

Толщина испаряющегося слоя лимитируется теплопроводностью, а тангенциальная скорость в нем — вязкостью. Кинематическая вязкость газа v «а, поэтому число Рейнольдса, рассчитанное по толщине слоя h, порядка 1. Известно, что ламинарный пограничный слой теряет устойчивость при числе Рейнольдса около 420 [42]. Это критическое число получено для ламинарного профиля Блазиуса. Однако вдув способен резко понизить устойчивость слоя, поскольку может порождать профили скорости с перегибом [116], неустойчивые при любых числах Рейнольдса. В том же направлении будет действовать шероховатость поверхности, которая по предположению должна возникать в описанном процессе, и заведомо высокая степень турбулентности в обтекающем потоке. С другой стороны, отвод тепла в частицу повышает устойчивость. Учесть все эти сложные процессы не представляется пока возможным. Поэтому изложенный механизм горения следует расценивать как гипотезу.

Остановимся на сравнении классических зависимостей и формулы (2.11). В работах Зельдовича и Вилюнова [110, 112] предполагается, что вдув не меняет существенно пограничный слойв частности, используются параметры течения для непроницаемой пластины. Это не согласуется с нашими условиями: например, при скорости обтекания 1 км/с и диаметре частицы d = 0,3 мм число Рейнольдса Re порядка 5 • 105, и толщина классического пограничного слоя была бы около 0,37c?/Re0'2 «0,027с? [42]. С учетом же вдува при скорости регрессии ип = 1 м/с, при которой параметр вдува F = psun/pV «0,017, имеем толщину порядка 5Fd & 0,085d. Видно, что влияние вдува чрезвычайно сильное. Можно сказать, что по сравнению с [110, 112] зависимость (2.11) описывает обратный предельный случай, когда обтекание не просто усиливает несколько теплопередачу, а меняет механизм горения. Можно ожидать, что закон (2.11) будет в определенном смысле предельным, поскольку он предполагает минимальность «теплового сопротивления» газовой прослойки между горючим веществом и продуктами реакции.

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

2.3.2 Инициирование. Роль дробленияДля моделирования резкого инициирования в небольшой части заряда задавалась область высокого давления (сотни атм.), как результат горения при постоянном объеме. Затем этот разрыв распадается (как при раскрытии мембраны) и начинается истечение газа в пористую среду. Поток разгоняет частицы и сжимает вещество. Поскольку тэн — довольно непрочный материал, он легко сжимается почти до плотности монолита и образуется непроницаемая «пробка». Этот процесс из-за больших начальных градиентов происходит довольно быстро (10 мкс) и может прекратить ПГД, задерживая течение газа. Через некоторое время легкая оболочка разрушается и процесс затухает.

В этой постановке необходим учет дробления частиц. Для существенного сжатия необходимо дробление, поскольку само ВВ практически несжимаемо при давлениях порядка сотен атм. Поры должны заполняться фрагментами исходных зерен. Это ведет к росту удельной поверхности и скорости реакции. Оценка показывает существенное увеличение поверхности, а моделирование демонстрирует резкий рост давления вблизи момента максимального сжатия. Это можно интерпретировать как формирование детонационной волны.

Рост удельной поверхностиВыбор коэффициента увеличения поверхности S — наиболее сложный момент в данной задаче. Приведем минимальную оценку этого эффекта.

При компактировании S растет: дробление производит частицы меньших размеров. Поскольку замена частицы, скажем, на две части не изменяет заметно S и (р основное внимание будем обращать на малые обломки размера d/k (происходящие из малых контактных областей), которые заполняют поры между исходными зернами. Значение масштабного фактора к явно больше 1, но, по-видимому, не превышает 10.

Оценку (2.12) можно назвать консервативной. Более резкому росту поверхности соответствуют данные [68]. С разумной точностью их можноописать зависимостью1 + 4Ч> Ч> о1 1(2.13)Например, при ip = 0,1 S 50.

Естественно, что зависимостями (2.12,2.13) следует пользоваться при сжатии порошка. При разгрузке удельная поверхность уменьшается: т. е. удельная поверхность при расширении пропорциональна объемной доле твердой фазы.

Начальные условияИсходно небольшой участок заряда (5 мм) заполнен продуктами горения под высоким давлением (несколько сот атм, 4000 К, зерна нагреты до ATs). В остальной части заряда — нормальные условия (293 К, 1 атм, пористость <ро = 0,42). Полная длина заряда обычно была 30 мм. Начальные скорости обеих фаз нулевые. Затем начинается течение газа, вызывающее движение твердой фазы. Это моделирует внезапное открытие диафрагмы, разделяющей две секции. Такое инициирование можно назвать «резким», поскольку начальный разрыв давления газа вызывает весьма быструю инжекцию горячего газа в пористую насыпку. Типичный расчет занимал менее 50 мкс.

Неопределенность расчетных параметров диктует необходимость рассмотрения большого числа вариантов. Для идентификации варианта примем сокращенные обозначения. Увеличение поверхности при дроблении будем обозначать Н — high, сильное, формула (2.13), С — консервативное, см. (2.12), либо N (нет дробления). Далее следует давление в секции инициирования (в сотнях атм) и ATs. Наконец, вариант резкого инициирования обозначается S. Например, первый вариант, рассмотренный ниже, обозначается C/3/14/S (консервативный/300 атм/14 К/резкий).

Результаты расчетовВ начале основной процесс — это фильтрация газа в порошок. Поток газа постепенно разгоняет частицы ВВ. При движении частиц создается тонкаяS = Smax ((1 — ф)/(1 — (pmin)) зона сильно сжатого тэна на некоторой глубине вблизи границы раздела. Эта стадия рассмотрена аналитически в предыдущей главе. В пренебрежении прочностью ВВ время коллапса обратно пропорционально инициирующему перепаду давления (около 10 мкс при 300 атм). В конце стадии коллапса давление в твердой фазе достигает уровня инициирующего давления.

Эволюция течения для варианта C/3/14/S показана на рис. 2.7. Профили объемной доли твердой фазы, а = 1 — давления р и скорости газа и показаны на момент 0,5 мкс (пунктир) и 10 мкс (сплошные линии). Метки на графиках, а и р обозначают место воспламенения, которое в этом расчете отстает от пика компактирования. Для момента 10 мкс показаны также профили скорости твердой фазы v, «твердого» давления ps и фактора роста поверхности S. Для последнего использована логарифмическая шкала. Перед волной 5 = 1- одно деление по вертикали соответствует увеличению 5 в 4 раза. Видно, что в течение 10 мкс почти достигнут коллапс пористой укладки в тонкой зоне заряда, где наблюдается сильная фрагментация. Эта картина мало зависит от входных данных и деталей модели. Необходимо только слабое сопротивление сжатию: прочность порошка должна быть мала.

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

Случай 200 атм (в центре) — промежуточный. Здесь пробка сформировалась, однако газ, профильтровавшийся в порошок перед этим и толкаемый вперед пробкой, в течение времени расчета успевает зажечь ВВ во фронте сжатия. Рост давления тоже быстрый, но начинается он поздно и с довольно высокого уровня, достигнутого в процессе медленного ускорения. В слабой оболочке на стадии ускорения следует ожидать срыва ПГД.

Без дробления (варианты.N, т. е. при выключении роста удельной поверхности) быстрый переход в детонацию в расчетах не происходил.

Влияние температуры воспламенения показывает рис. 2.12 для умеренного роста поверхности и инициирующего давления 300 атм (варианты C/3/(3.5,7,14)/S). Среднее значение 7 К (средний кадр) представляется оптимальным. При повышении ATs до 14 К (справа) образуется пробкаэтот вариант показан также на рис. 2.7 и 2.8. Интересно, что при малой ATs = 3,5 К (слева) развитие процесса сравнительно медленное. Здесь раннее воспламенение приводит к выделению газа в количестве, достаточном для замедления сжатиясоответственно эффект дробления ослаблен.

20 t, мкс 40 020 мкс 40 0 Рис. 2.12.

20 t, мкс 40Аналогичные вычисления для случая сильной фрагментации дают быстрый ПГД во всех случаях (рис. 2.13, варианты Н/3/(3.5,7,14)/S). Более того, при ATs = 14 К переход наиболее резкий: при малом содержании газа сжатие сильное, и поверхностный эффект значителен.

В некоторых расчетах начальный разрыв по давлению газа заменялся плавным изменением на размере 10 мм. Начальное течение газа и соответственно компактирование были заметно медленнее. При умеренном росте поверхности типично образование пробкиначало быстрого роста давления сдвигалось на 5 — 10 мкс. Для сильного роста поверхности реализовывался режим без пробки. Переход запаздывал на 5 — 7 мкс, но становился более резким. В общем, для сглаженных начальных условий быстрый ПГД менее вероятен.

Не удалось наблюдать в расчете быстрый ПГД при малых инициирующих давлениях. Например, при десятках атм резкие изменения начинаются после «100 мкс.

Таким образом, быстрый переход происходит при совместном действии двух инициирующих агентов: потока тепла от горячих продуктов к частицам ВВ, который имеет место и в классическом сценарии, и резкой механической активации среды в результате коллапса пористой структуры, которая существенно ускоряет процесс.

Степень увеличения поверхности — это основной фактор, определяющий возможность быстрого перехода. Консервативная оценка (2.12) срабатывает в наиболее жестких условиях. Формула (2.13) приводит к быстрому переходу почти рутинно. Реальная ситуация находится между этими крайними случаями, что показывают недавние результаты [117]. В этой работе получено, что удельная поверхность растет линейно с увеличением а. При сжатии до о. = 0,9 фактор S возрастает приблизительно до 6,5. Это примерно в 1,5 раза более сильная зависимость, чем (2.12). Отметим, что разные методы измерения дали существенно разные (на порядок) величины удельной поверхности, но пересчет к относительной величине S заметно уменьшил разброс.

2.3.3 Роль эрозионных эффектов и неустойчивости на поверхности зернаОписанная модель позволяет исследовать в расчете поведение системы при «включении» описанного выше механизма быстрой регрессии. Различие с традиционными закономерностями оказалось довольно резким. Опишемрезультаты расчетов для случая инициирования пересжатой газовой детонацией. Этот процесс изучался в [76], а в последнее время также в [78]. Принимаем, что в начальный момент область высокого давления свободна от частиц, а заполнена горячим газом с параметрами, соответствующими отражению пересжатой детонационной волны от жесткой стенки. Такая постановка при задании стандартного закона горения (2.10) и тех же начальных давлениях дает результаты, очень близкие к случаю «мембранного» инициирования, который рассмотрен в предыдущем параграфе. Экспериментально газовое инициирование удобно тем, что процесс гораздо более воспроизводим.

При типичной постановке (диаметр зерна 0,3 мм, скорость горения соответствует (2.11)) практически независимо от начального давления в порошок распространяется волна со скоростью 800−900 м/с. Давление в этой волне близко к инициирующему, то есть развития и ускорения процесса не происходит. На рис. 2.14 показаны результаты расчета на момент времени 21,5 мкс при инициирующем давлении 2000 атм. Это больше, чем достигается в экспериментах (например, при начальном давлении газовой среды 1 атм давление инициирования около 270 атм). Однако этот случай интересен тем, что конечное состояние представляет собой квазистационарную волну с параметрами, близкими к наблюдавшимся в работах JI.A. Лукьянчикова и В. В. Андреева. Расчетная скорость волны 890 м/с, уровень давления около 2 кбар. На врезке рис. 2.14 показаны экспериментальные осциллограммы: массовой скорости (а) и давления, измеренного пьезодатчиком (б), приведенные к одному масштабу по давлению. Между ними помещена расчетная осциллограмма массовой скорости (в), также пересчитанная на давление. Амплитуда и время нарастания сигнала недалеки от эксперимента. В экс20 X, ммРис. 2.14.перименте сигнал довольно резко спадает после максимума из-за бокового разлета, тогда как одномерный расчет этот эффект не описывает.

Сжатие скелета в волне незначительно (см. профиль а). Поэтому невелик и эффект дробления: степень увеличения поверхности несущественна. Зона реакции соответствует участку роста давления: именно в этой области заметно различие скоростей фаз, которое определяет скорость регрессии согласно (2.11). Причина столь резкого отличия от стандартной постановки — в том, что для механизма (2.11) сразу при внедрении горячего газа в порошок начинается интенсивное газовыделение. Поэтому среда сопротивляется сжатию, и коллапс не происходит.

При уменьшении размера частиц до 100 мкм скорость волны заметно не изменяется. Однако давление в ней постепенно растет. Подобное же влияние оказывает ускорение регрессии вдвое. Интересно, что рост происходит за счет зарождения возмущений позади фронта волны. Малое возмущение скорости проскальзывания фаз включает реакцию в этой зоне, на профиле давления возникает положительный «зуб». Эта неустойчивость, развиваясь, догоняет фронт, повышая уровень давления. Поскольку реальные эксперименты проводились с зарядами малого диаметра в слабой оболочке, такой механизм развития является нефизическим. Можно сделать вывод, что в предположении справедливости (2.11) «абляционный» процесс может осуществляться в квазистационарной волне, но он недостаточен для ее развития.

Напротив, эрозионные эффекты (которые рассчитывались по Вилюно-ву) при классическом механизме горения не приводили к качественно новым эффектам. Процесс мог, например, ускориться на 1 мкс при типичном времени 15 мкс. Это можно объяснить тем, что в стандартной постановке значительная скорость проскальзывания имеет место в той части волны, где еще не началось горение.

На рис. 2.15 показаны расчетные осциллограммы массовой скорости при варианте расчета, отличающемся от C/3/7/S, подробно обсужденного в предыдущем параграфе, только учетом эрозии (и отсутствием ВВ и инициирующей секции). Цифры под кривыми означают заглубление датчика в ВВ (верхний ряд) и время, соответствующее началу каждой записи (нижний ряд). Результаты качественно согласуются с полученными в эксперименте (см. работы [77, 78], откуда взято расположение датчиков). Сравнение результатов проведено в таблице 2.3. Уровни массовой скоросТаблица 2.3.h, мм 0 1,5 3 8,7 10,8Vmax раСЧ., м/с 96 102 113 500 500vmax эксп., м/с [78] 140 170 140 520 235tf расч., мкс — 5 7 0,7 1tf эксп., мкс [78] 1,4 3,7 3,4 1,4 2,1ти отличаются от эксперимента не более чем в полтора раза, а времена нарастания — не более чем вдвое. Задержка между сигналами с датчиков, расположенных на глубине 1,5 и 8,7 мм, равна 15 мкс в эксперименте и 17,75 мкс в расчетедля датчиков при 3 и 10,8 мм имеем соответственно 17,5 мкс и 16,58 мкс.

Такое согласование представляется вполне приемлемым, если учесть экспериментальный разброс (например, скорость при 10,8 мм вдвое ниже полученной в другом эксперименте при 8,7 мм). Таким образом, двухфазное инициирование качественно описывается моделью коллапса, учитывающей дробление ВВ.

В то же время квазистационарная двухфазная волна лучше описывается «абляционной» моделью. Модель коллапса со стандартным законом регрессии приводит либо к быстрому инициированию, либо к формированию пробкиквазистационарный режим не получается. Пока не ясно, как увязать эти предельные режимы: необходимы дополнительные исследования. Очевидно, что важную роль играет боковое расширение (двумерность).

0 1,36 3,92 19,1 20,5 мкс Рис. 2.15.

2.4 Конвективные волны в жесткой пористой средеВ этом разделе рассмотрим квазидетонационные волны, наблюдавшиеся в экспериментах В. В. Митрофанова и др. [81], А. В. Пинаева и Г. А. Лямина. Изученные этими авторами реагирующие системы можно рассматривать как предельный случай, в котором исключена деформация пористого скелета. Это существенно упрощает анализ, и можно надеяться на продвижение в понимании и более сложных процессов.

Вначале мы рассмотрим основную континуальную модель течения, затем применим ее к детонации ВВ в жесткой пористой среде [66, 86], а затем iIк более сложному случаю детонации газ-пленка [82, 83, 84, 85]. Кроме того, ^ будут изложены результаты дискретного подхода — модели решеточного s газа.

2.4.1 Континуальная модель. ЗаторможенностьВ традиционном представлении детонация — это одномерная квазистационарная волна. Хотя такое описание неверно даже для многих гомогенных систем (например, для реальной многофронтовой газовой детонации), оно может применяться как некоторая средняя картина. Осредненные газодинамические уравнения с успехом используются и для гетерогенных сис-. тем. Этот подход имеет недостатки, особенно очевидные для волн с кон- ' вективным механизмом распространения: исчезают инициирующие струи. Безусловно, применимость осреднения при описании течений, в которых настолько существенна случайная компонента, ограничена. Тем не менее, при определенных условиях использование осреднения оправдано. Кроме того, континуальная модель является хорошей отправной точкой, позволяющей перейти к дискретной модели, рассматриваемой в следующем разделе.

Будем предполагать, что волна имеет нестационарный предвестник, где происходит поджигание топлива и распространение горения на все сечение заряда. Для этой части волны осредненное описание недостаточно. К предвестнику примыкает основная зона реакции, для которой уже можноприменять осредненный подход. Разумеется, и здесь присутствуют флуктуации течения («пористая турбулентность»), но в зоне реакции они являются второстепенным фактором.

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

Здесь р — суммарная плотность газа, и — скорость течения, р — давление, Q — теплота реакции горения при постоянном объеме (на единицу массы топлива). Так как пористость постоянна, она не входит в уравненияпри этом массоприход j сила трения /, внутренняя энергия смеси Е и теплоотвод q определены на единицу объема пор. Переменная г обозначает часть плотности газа, порождаемую поступлением со стенок в поры топлива (по предположению сразу и полностью сгорающего). В частном случае детонации ВВ можно считать, что г ^ р, а при газопленочной детонации р> г > 0.

Будем считать, что показатель адиабаты 7 газа в порах постоянен. Если поры заполнены смесью продуктов сгорания и непрореагировавшего окислителя, считаем оба газа идеальными с одинаковыми 7. Такое предположение резко упрощает анализ и одновременно близко к истине (по крайней мере, неточность его представляется несущественной на фоне гораздо более грубых приближений, особенно касающихся межфазных взаимодействий). Тогда Е и для смеси выражается через давление р одинаковое для обоих компонент, независимо от распределения объема между газами и их перемешивания: Е = р/(7 — I).

Уравнения (2.14) записаны в системе отсчета, связанной с пористым скелетом. Поступающее топливо в этой системе имеет нулевую скорость ине вносит вклада в импульс газа. Можно показать, что в зоне реакции теп-лоотвод q не успевает заметно повлиять на течение: потери оцениваются в несколько процентов. Поэтому, рассматривая зону реакции и ее окрестность, далее будем пренебрегать потерями тепла.

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

Укороченная модельНачнем с «нулевого» приближения, т. е. полностью опустим конвективные слагаемые (и = 0). Сравнивая изменения массы топлива и энергииrtaj и EtajQ, (2.15)получаемE^rQ + Eo, (2.16)где Ео = ро/(у— 1) — начальная тепловая энергия. Смысл равенства (2.16) в том, что реакция идет практически при постоянном объеме [81].

Следовательно, равенство (2.16) с хорошей точностью выполняется также и в первом порядке по и. Для давления получаем соотношениеp"(7-l)rQ + p0, (2.18)которое заменяет уравнение энергии.

Температура же продуктов горения в зоне реакции меняется в гораздо более узких пределах. Действительно, в однофазной детонационной волне выделение тепла реакции приводит к нагреву среды, чем и определяется рост температуры. В рассматриваемой же системе основной разогрев происходит при горении в тонких пограничных слоях на поверхности частиц, которые исключаются из рассмотрения. В поровое пространство поступают уже «готовые» продукты горения. Этот процесс обеспечивает в зоне реакции своего рода термостат. Влияние на температуру изменяющихся условий (давления, плотности) в зоне реакции невелико. Из-за затрудненности движения в пористой среде внутренняя энергия не может перейти и в кинетическую.

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

За зоной реакции температура будет спадать из-за теплоотвода в скелет. Однако этот спад сравнительно медленный. В зоне горения роль теплоотвода имеет порядок отношения объема пограничных слоев к объему пор.

Сразу отметим, что сходство графиков с классической картиной в большой степени внешнее. Прямая (2.22) не является прямой Релея — Ми-хельсона: по мере поступления массы состояние газа отнюдь не меняется вдоль нее. Скорее линия (2.22) — это аналог детонационной адиабатына нее точка, изображающая состояние газа, попадает, только набрав определенные значения J и F. Поэтому указанную прямую будем называть детонационной адиабатой. В процессе эволюцииизображающие точки движутся по гиперболе состояния р = c2/V. Но и эта линия не будет аналогом ни детонационной адиабаты, ни линии Релея, так как на ней находится вообще любое состояние за волной, в том числе и вне стационарной зоны.

Тем не менее физически задача подобна классической. Приток массы J, в результате которого волна является самоподдерживающейся, аналогичен теплу реакции Q, тогда как трение F определяет потери. Разумеется, употребление термина «адиабата» для описанной модели в известной степени противоречиво.

К VРис. 2.16.

Из (2.26) для Fe > О D < 2с. Поскольку для рассматриваемых волн скорость и мала, реально около 0,1с, (2.28) дает наилучшую оценку скорости волны D т. е. D близка к с. Следовательно, при доминирующей роли трения скорость волны практически фиксирована и пользы от выражения (2.27) для D немного. Это отличает ситуацию от классической, когда условие касания определяет единственную скорость D, через которую определяются остальные величины. Вместо D хорошим параметром будет Fe. Из (2.26) для данного притока массы Je определяется значение Fe и тем самым неявно уровень скорости и, достигаемый в волне.

В практически интересном случае Je c/Vq или D/Vq (масса продуктов велика по сравнению с начальной массой газа), из (2.25,2.26) р cJe ро V и с/ Je <С V0 Fe& р. Из (2.26) D & 2cFe/ Je (что по существу означает D «с). Эти результаты не зависят от вида функции j (при естественном условии j > 0).

Изотермическая модель предсказывает величину скорости детонации D=/(y-l)Q. Такой результат для режима, предельного по потерям, получен В. В. Митрофановым [118], причем использовалось уравнение энергии. В данной работе показано, что режимы детонации в стесненных условиях автоматически близки к предельным, а также отмечена приближенная изотермичность зоны реакции, что значительно упрощает анализ. Близкий результат в несколько других предположениях получен также С.С. Рыба-ниным [119].

Характерная величина скорости D для типичных ВВ близка к 1 км/с. Экспериментальные данные имеют заметный разброс вокруг этого значения, чего и следует ожидать от приближенной модели.

Легко проверить, что в уравнении импульса (2.14) отброшенные при оценках члены (pu)t и (ри2)х малы по сравнению с основными — градиентом давления и силой трения. Действительно, (ри2)х pu2/L <С /, а (pu)t рис/L рхи/с <С рх. Заметим, что «заторможенное» приближение подтверждено постфактум. Это не случайно: приближение работает на квазистационарных решениях. Можно представить себе ситуацию, когда оно не выполняется, например если состояние фаз однородно в пространстве, но имеется значительная скорость проскальзывания.

Эволюция нестационарной волныПоскольку правые части уравнений зависят от решения, результаты предыдущего раздела говорят лишь о возможности описанного режима. Продемонстрировать его реальность можно в полной нестационарной постановке. Для анализа структуры волны и ее развития из заданных начальных условий численно решались уравнения (2.20). Для простоты, чтобы избежать явного рассмотрения струйного предвестника, в данных расчетах горение включалось ударной волной в порах.

Применена эйлерова разностная схема Лакса — Вендроффа. Неоднородные члены аппроксимировались с сохранением второго порядка точности. Фронт УВ выделялся и рассчитывался с учетом условий на разрыве и на догоняющей характеристике. Точки, примыкающие к фронту волны, которые не поддаются процедуре Лакса-Вендроффа, рассчитывались по специальному алгоритму [120]. В расчетах принималось с = fbjl (как в решеточной модели, излагаемой далее) — расчетная область 0 < х < 200. Пространственный шаг h = 1 понимался также как характерный внутренний размер пористой среды d.

Скорость реакции задавалась в виде j = j (tj) (tj (x) — время, прошедшее с момента, когда УВ прошла через данную точку х) в двух вариантах. В первом (квазипостоянная скорость реакции) j = Cj/tr (0 < tf < Q, Str), j = (Cj/tr) • cos2(1,25tTtf/tr) (0,8tr < 1,2tr).

Здесь tr — характерное время реакции, Cj — коэффициент, регулирующий интенсивность массоприхода. Такая зависимость моделирует постоянную скорость выделения массы за УВ в течение времени tr сглаженную в конце зоны реакции (при 0,8tr.

Рис. 2.18.

Вариант расчета выбран таким образом, чтобы продемонстрировать и ударную волну, и рост давления за ней. Более реальна постановка с малой начальной плотностью (или интенсивным горением), при которой скачок в ударной волне незначителен. Роль ударной волны тогда сводится к включению реакции, что на практике обеспечивается неодномерностью течения (струйным механизмом).

В расчетах горение имело «безусловный» характер. Нетрудно ввести срыв реакции, скажем, по достижении некоторого значения напряжений или теплового потока в пограничном слое. Тогда в зависимости от трения можно получить условие развития самоподдерживающейся волны. Учет теплоотвода приведет к более быстрому спаду давления.

2.4.3 Детонация системы газ-пленкаПерейдем к случаю «дважды гетерогенной» системы, в которой газообразный окислитель, заполняющий объем пор, реагирует с пленкой топлива на поверхности поровой структуры [82, 83, 84, 85]. Этот случай — более сложный по сравнению с детонацией унитарного топлива, рассмотренной выше. В частности, здесь не проходит концепция «вакуумной» детонации [67] (подробнее обсуждаемая ниже, в п. 2.5).

В континуальном подходе основной системой уравнений будет (2.19), а ключевым параметром волны — скорость «звука» с = s^/r/p, где s = у (7 — 1) Qдавление в волне согласно укороченной модели следит за составом: р = pQ—s2r. Для сравнения с укороченной системой будет рассмотрена и стандартная модель (2.14).

Теплота сгорания реальных топлив примерно на порядок больше, чем у ВВ (главным образом из-за того, что Q рассчитывается на единицу массы топлива, тогда как ВВ содержит «встроенный» окислитель и балласт типа атомов азота). При Q = 11,5 ккал/г и 7 = 1,3 константа s 3,8 км/с.

Однако выражение для с содержит теплоту Q, отнесенную к единице массы газа в целом. Поэтому для газопленочной системы с заметно меньше, чем s так как увеличение плотности газа в волне невелико. Например, для стехиометрической смеси топлива брутто-состава СН2 с кислородом масса горючего (г) составляет 14/62 «0,23 от суммарной массы продуктов (р), так что с < s/2. В результате конечное (максимальное) значение се 1,8 км/с — несколько больше, чем при горении ВВ.

Структура волныДля экспериментов [82, 83, 84, 85] характерно плавное нарастание давления в волне. Рост давления определяется постепенным выгоранием топлива.

Для структуры волны существенно, что в газопленочной системе величина с отнюдь не постоянна в зоне реакции — она растет по мере выгорания топлива (формально от нуля до максимального значения, которое s/2). В этом состоит главное отличие от однокомпонентной системы. (Разумеется, растет и скорость звука для стандартной модели С = y/ур/р, которая в зоне горения примерно в y/у раз больше).

Рассмотрим стационарную волну, имеющую скорость D 1 км/с, предполагая, что скорость течения и <С D. Поскольку скорость звука за фронтом волны растет, начиная с малого значения, течение классической структуры — дозвуковое в зоне реакции с переходом к сверхзвуковому в точке Жуге — оказывается невозможным. В начале зоны реакции обязательно будет D — и > с, а переход через звук, если и происходит, то в «обратном» направлении. (В предыдущем разделе для изотермического случая при постоянной с решение удается построить за счет малых изменений и).

Обсудим допустимые структуры. Во-первых, скорости и и D могут быть близкими. Разумеется, при этом приближение «заторможенности» не работает, и следует рассматривать полную систему (2.14). Если на фронте волны и D, то к фронту может примыкать дозвуковая область. При последующем спаде и в принципе не исключена стандартная структура течения. Этот случай аналогичен классическому решению, например, для газовой детонации, где по мере выделения тепла скорость звука также растет. Однако, так как скорость и велика, трение будет значительным. Необходимость «проталкивать» пробку газа должна ограничить скорость волны. Поэтому следует ожидать полностью дозвукового решения: D — и будет меньше С даже при полной остановке течения, т. е. при и = 0. Такой режим можно назвать «медленной волной с быстрым течением».

Во-вторых, сохраняя предположение о заторможенности, можно вообще отказаться от дозвуковой зоны, примыкающей к фронту. Вполне мыслима инвертированная структура. Если распространение волны происходит за счет проникновения опережающих струй, то скорость ее определяется не средней динамикой, а мезопроцессами в порах. По существу, такая волна аналогична известному примеру недосжатой детонации, в которой принудительным инициированием задается высокая скорость фронта. В таком случае можно говорить о «быстрой волне с заторможенным течением».

В обоих случаях в результате выделения тепла конечная скорость звука велика (около 2 км/с, тогда как скорость волны в экспериментах не превышала заметно 1 км/с). Следовательно, в конечном состоянии течение дозвуковое. Отметим, что особая роль звуковой точки и режима Жуге имеет газодинамическое происхождение. Для пористой среды правильнее говорить о фильтрационном течении, для которого газодинамические закономерности могут быть недействительны. Это относится и к способу распространения процесса, и к влиянию возмущений снизу по потоку, действие которых в пористой среде быстро затухает.

В начальный момент в небольшой области вблизи начала координат инициировалась реакция. Рост давления приводил к формированию практически стационарной волны по истечении 2−3 характерных времен реакции.

Р.ГПл р, г/см3 0,080,04/?, ГПа р, г/см38 ЛГ, СМи, А км/с 1,20,4 08 ЛГ, СМРис. 2.19. а — dp = 0,3 мм, t = 176 мксб — dp = 2,5 мм, t = 105 мкс.

На рис. 2.19 показаны профили волны, рассчитанные по стандартной модели. Здесь и далее профиль j (для удобства приподнятый) отмечает расположение зоны реакции. На рис. 2.19,а представлен вариант, моделирующий волну в порах песка (дисперсность в [82] была 0,15 — 0,5 мм, в расчете задано среднее значение 0,3 мм). На рис. 2.19,6 демонстрируется расчет для укладки металлических шариков со сравнительно крупными порами [84].

Рост давления в волне постепенный, кроме небольшого начального скачка. Такая зависимость соответствует эксперименту. Плотность газа резко увеличивается в ударной волне, а далее спадает, несмотря на поступление массы топлива. Как уже отмечалось, в стандартной модели скорость течения и порядка скорости волны D, уровень которой отмечен на рисунках. Для крупнодисперсной среды на фронте волны и «D. Заметим, что при dp = 2,5 мм длина зоны реакции составляет лишь несколько диаметров частиц. При этом, разумеется, трение не так существенно, как для мелкозернистой среды.

Для крупных пор расчетная скорость волны 807 м/с недалека от экспериментальной (940 м/с для укладки шариков [85]). Однако для мелкодисперсной среды отклонение больше (расчет 490 м/с, эксперимент «850 м/с для песка [82]), т. е. расчет дает слишком сильную зависимость от дисперсности.

Профили, А = С + и на фоне уровня скорости волны D подтверждают заключения, качественно обоснованные выше. Скорость волны сравнительно невелика, а скорость звука С = fypfp в конце реакции примерно та же, что для детонации газовзвеси такого же состава. На фронте также С + и > D как и должно быть за ударной волной. Поэтому все течение за фронтом дозвуковое. Следовательно, волны, полученные в расчете, отнюдь не являются режимами Чепмена-Жуге.

Рассмотрим теперь результаты укороченной модели. Для нее ударная волна как инициирующий агент не годится. При струйном механизме распространения отбор средней скорости волны определяется процессами на уровне отдельных пор. Сейчас такие процессы практически не изучены, и пока трудно увязать скорость волны с параметрами среды. Поэтому в данной модели скорость волны задавалась принудительно (в соответствии с экспериментом) равной 1 км/с. На фронте волны задавалась нулевая степень выгорания и нулевая (средняя) скорость и. Примеры расчетар, ГПа р. г/слг3 0,08.

Основные результаты работы опубликованы в следующих статьях:

1. Ершов А. П. Быстрые движения газа в пористой среде // Журнал прикладной механики и технической физики. 1983. № 1. С. 65−69.

2. Ершов А. П. Об уравнениях механики двухфазных сред // Журнал прикладной механики и технической физики. 1983. № 6. С. 79−86.

3. Ершов А. П., Куперштох A.JI. О температуре продуктов детонации при взрыве в камере // Физика горения и взрыва. 1986. Т. 22, № 3. С. 118−122.

4. Ершов А. П. Растущие волны в активной двухфазной среде //В кн.: Проблемы нелинейной акустики. Труды И Международного симпозиума IUPAP-IUTAM по нелинейной акустике. Ч. 2. Ред. Кедринский В. К. Новосибирск, 1987. С. 26−28.

5. Ершов А. П., Куперштох A.JI., Коломийчук В. Н. Образование фрактальных структур при взрыве // Письма в Журнал технической физики. 1990. Т. 16, № 3. С. 42−46.

6. Ершов А. П., Куперштох A.JI. Образование фрактальных структур при взрыве // Физика горения и взрыва. 1991. Т. 27, № 2. С. 111−117.

7. Ершов А. П., Куперштох A.JI. Экзотермическая коагуляция малых кластеров во фронте детонационной волны // Письма в Журнал технической физики. 1993. Т. 19, №. С. 76−80.

8. Ershov А.P., Kupershtokh A.L., Dammer A.Ya. Fingering in the fast flow through porous medium // Journal de Physique II (France). 1993. V. 3, N 7. P. 955−959.

9. Ершов А. П. Газодинамика клеточных автоматов // Физика горения и взрыва. 1994. Т. 30, Ш. С. 107−117.

10. Ершов А. П. Изотермическая детонация и ее стохастическое моделирование // Физика горения и взрыва. 1994. Т. 30, № 3. С. 112−124.

11. Ершов А. П., Иванов П. И., Андреев В. В. Детонация сильноразбавлен-ного ВВ // Физика горения и взрыва. 1994. Т. 30, № 3. С. 124−130.

12. Ershov А.P., Kupershtokh A.L., Dammer A.Ya. Structured flows in porous media modeling by stochastic methods // In Numerical Methods for Multiphase Flows. Proc. Int. Symp. on Multiphase Flows. Crowe C.T. et al (Ed). Lake Tahoe, NV, 1994. N.Y.: ASME, FED-Vol. 185,. P. 59−64.

13. Ershov A.P. Isothermal detonation // Combustion and Flame. 1995. V. 101, N 3. P. 339−346.

14. Ershov A.P. Convective combustion of granular explosive after sharp initiation — fragmentation effect // In Multiphase Flow '95. Proc. 2nd Int. Conf. on Multiphase Flow. Serizawa A. et al (Ed). Kyoto, 1995. P. CO-27 — CO-31.

15. Kupershtokh A.L., Ershov A.P., Medvedev D.A. Coagulation of carbon clusters in detonation front //In Shock Compression of Condensed Matter — 1995. AIP Conf. Proc. 370. Part 1. Schmidt S.C., Tao W.C. (Eds). Seattle, WA, 1995. AIP, 1996. P. 393−396.

16. Ершов А. П. Компактирование непрочной пористой среды газовым поршнем // Журнал прикладной механики и технической физики. 1996. Т. 37, т. С. 156−164.

17. Ершов А. П. Конвективная детонационная волна в пористой структуре // Физика горения и взрыва. 1997. Т. 33, № 1. С. 98−106.

18. Куперштох A.JI., Ершов А. П., Медведев Д. А. Модель коагуляции углеродных кластеров при высоких плотностях и температурах // Физика горения и взрыва. 1998. Т. 34, № 4. С. 102−109.

19. Ершов А. П., Куперштох A.JI., Медведев Д. А. Моделирование волн горения в пористой среде методом клеточных автоматов // В кн.: Динамика сплошной среды. Ред. Кедринский В. К. Вып. 115. Новосибирск, Ин-т гидродинамики, 1999. С. 74−79. Акустика неоднородных сред. Доклады 5-го семинара СНГ.

20. Ershov A.P., Dammer A.Ya., Kupershtokh A.L. Inviscid finger instability in an anisotropic porous medium // In Modern Approaches to flows in porous media. Internat. Conference dedicated to P.Ya. Polubarinova-Kochina. Abstracts. Moscow, 1999. Institute for Problems in Mechanics RAN. P. 49−51. (Современная теория фильтрации. Международная конференция памяти П.Я. Полубариновой-Кочиной (1899−1999)).

21. Ершов А. П., Сатонкина Н. П., Дибиров О. А., Цыкин С. В., Янилкин Ю. В. Исследование взаимодействия компонент гетерогенных взрывчатых веществ методом электропроводности // Физика горения и взрыва. 2000. Т. 36, № 5.

22. Ershov А.P. Modeling of the deflagration to detonation transition in porous PETN //In Proc. 11th Symposium (International) on Detonation. Snowmass Village, CO, 1998. Ampersand Publ. Group, 2000. P. 686−692. http://www.sainc.com/onr/detsymp/FinalPapers/DE183.pdf).

ЗАКЛЮЧЕНИЕ

.

В самом кратком виде результат данной работы — это построение физической картины детонационных процессов в двух важных случаях:

• для неидеальных квазидетонационных волн в пористых средах и в двухфазных средах типа газ-частицы;

• для процесса выделения избыточного углерода при детонации плотных ВВ и роста ультрадисперсных углеродных частиц.

Эти явления критически зависят от мезопроцессов и не могут быть поняты без их учета.

Перечислим более детально основные результаты работы:

1. Предложена уточненная модель механики двухфазных сред типа «газ-частицы», в которой регуляризована коротковолновая неустойчивость.

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

3. Введено понятие неустойчивости «невязкого пальца» при вытеснении плотного флюида менее плотным из пористой среды, указанная неустойчивость исследована в изотропном и анизотропном случае.

4. Экспериментально исследована детонация сильно разбавленного ВВ. Найдено, что при малых концентрациях ВВ (до 10%) процесс распространяется без ударной волны;

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

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

Т. Показано, что при синтезе алмаза из углерода ВВ существенна роль газодинамических процессов в стадии сильного расширения ПД.

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

9. Экспериментально изучено взаимодействие компонентов смесевых ВВ при взрыве. Установлена зависимость электропроводности гетерогенных ВВ от дисперсности. Численные расчеты взаимодействия демонстрируют роль неустойчивостей течения и находятся в согласии с экспериментом, а также с известными изотопными измерениями.

10. Показано, что динамические эксперименты по рассеянию синхротрон-ного излучения хорошо согласуются с концепцией быстрого формирования алмазных частиц за фронтом детонации.