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

Вычислительные технологии

Том 5, № 4, 2000

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

Ю.А. Гапоненко Институт вычислительного моделирования СО РАН

Красноярск, Россия e-mail: adrian@cckr.krasnoyarsk.su

An appearance of a cumulative jet at a detonation of two flat parallel charges is numerically investigated. For calculations, the modification of Godunov — Rodionov scheme has been employed. The model of the medium represents the mixture of perfect gases with a computed effective isentropic exponent.

Введение

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

Данный процесс был исследован в работах [3, 6] методом численного моделирования на задачах детонации зарядов с плоскими и цилиндрическими каналами. Из результатов, приведенных в этих работах, очевидно наличие кумулятивного течения продуктов детонации, двигающихся против основного течения продуктов взрыва со скоростью, превышающей скорость детонации взрывчатого вещества. Однако в вычислительной модели, используемой в этих работах, показатель адиабаты, несмотря на его различие для продуктов взрыва и окружающей среды, принимался постоянным во всей расчетной области. Такое упрощение не только усложняет задание граничных условий на фронте детонации, но и не раскрывает полностью структуры образующегося течения.

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

© Ю. А. Гапоненко, 2000.

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

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

1. Численная модель

1.1. Постановка задачи (граничные условия)

Рассмотрим постановку задачи [3, 6] для случая детонации двух плоских параллельных зарядов. Выбор плоской, а не цилиндрической, постановки задачи обусловлен тем, что в этом случае процесс формирования кумулятивной струи происходит менее интенсивно, что значительно облегчает исследование.

Взрывчатое вещество (далее ВВ) ABCE (рис. 1) расположено над поверхностью симметрии FG. Участок границы (фронт детонации) симметричен относительно границы OH (точка O находится посередине ), поэтому в расчетах учитывается только участок OA, на котором задаются условия вдува продуктов взрыва (далее ПВ) с параметрами точки Чеп-мена — Жуге детонирующего ВВ. На участке были поставлены условия непротекания, которые остаются справедливыми с физической точки зрения, пока отраженная от FG ударная волна (УВ) не достигнет линии OH. По этой причине размеры расчетной области подбирались таким образом, чтобы отраженная УВ выходила из расчетной области на участке HG. Участок HG — "фиктивная граница", являющаяся "прозрачной" для возмущений, выходящих из области, в связи с чем, на ней были заданы модифицированные временно-зависимые "неотражающие" условия Томпсона [9]. В данной постановке задачи, как ив [3, 6], не учитывалось взаимодействие непрореагировавшего ВВ с кумулятивной струей, поэтому на участке AE были поставлены условия непротекания.

в

и 1 A | 1 1

Y 1 1 D-► 1 J->X s" s

Рис. 1. Расчетная область.

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

В начальный момент времени область АОИОГЕ заполнена двигающимся со скоростью детонации О газом, состояние которого определяется следующими параметрами: давлением р0, плотностью р0 и показателем адиабаты 70. На участке фронта детонации ОА за-

даются параметры ПВ, соответствующие точке Чепмена — Жуге: давление ра, плотность Ра, показатель адиабаты 7а и скорость истечения ПВ па-

Были проведены расчеты для детонации ВВ (гексоген) с параметрами в точке Чепмена-Жуге [1, 3]: давление ПВ ра = 213 кбар, скорость детонации Д = 7.3 км/с, плотность ПВ Ра = 1-96 • 103 кг/м3, скорость истечения ПВ па = 5.2 км/с, показатель адиабаты ПВ 7 = 2.63. Начальные параметры газа, окружающего заряд, таковы: плотность р0 = 1.16 кг/м3, скорость и0 = Д = 7.3 км/с, показатель адиабаты 7 = 1.4. Давление внешней среды р0 в начальный момент времени в расчетах изменялось от 1 до 100 бар. Изменялась также и толщина заряда.

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

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

др + д(ри) +д м = 0

д£ дх ду д (ри) д(ри2 + р) д (ри^) о д£ дх ду ' (1)

д (р^) д (ри^) д(р^2 + р) о

д£ дх ду

де д(е + р)и д (е + р)^ о

д£ дх ду

/ и2 + V2 \ р

е = р е + —й— Ь Ф'Р)

2 Г Р(7 - 1)'

где р — плотность; и, V — компоненты скорости вдоль осей х и у; р — давление; е — полная энергия; е — удельная внутренняя энергия. Для нахождения эффективного показателя адиабаты вместе с (1) решалось уравнение:

д (Р7) + д(Ри7) + д (PV7) = 0 (2)

д£ дх ду

означающее отсутствие механизма диффузии (чистый "перенос") отдельных компонентов смеси.

1.3. Метод решения

Для численного интегрирования системы уравнений (1), (2) была применена разностная схема (РС) А. В. Родионова [7] второго порядка точности по времени и почти всюду по пространству типа предиктор-корректор.

Используемая в работе РС является распространением схемы [7] на случай двумерных нестационарных уравнений. Данное обобщение не представляет значительной сложности, и поэтому подробное описание РС здесь опущено. Основное отличие применяемой РС от схемы [7] состоит в способе определения приращений параметров на этапе предиктора при задании кусочно-линейной интерполяции на текущем временном слое. Выбор способа

определения приращений параметров обусловлен результатами решения ряда тестовых задач. В данной работе для их определения использовался монотонизирующий алгоритм, предложенный в [8] и примененный впоследствии, например, в [2]:

= )(ff-1/2,j), A^j = *(Sij )(f™-1/2),

f _ i n i „ „'»T fn _ fn fn fn _ fn cn _ fn /fn

где 1 — (p,u,v,P,I) , fi-1/2,j — fi,j — fi-1,j, Ii,j-1/2 — fi,j — fi,j-1, ri,j — fi+1 /2 j/ fi-1 /2,j, Sij — fiilj+1/2/fiilj-1/2. Здесь целыми индексами (i,j) помечены величины в центрах ячеек (средние значения). Индекс n соответствует текущему временному слою tn. Индексы x, y обозначают координатные направления приращений параметров.

Функция Ф(г), отвечающая за наклон линейного профиля в ячейке, одновременно исполняет роль ограничителя производной и монотонизатора решения. В зависимости от вида функции Ф(г), как показано в [8], может быть построено семейство РС, имеющих второй порядок точности по пространству в областях гладкости решения, а также обладающих высокой степенью локализации газодинамических разрывов и свойством "монотонности" решения на них. В данной работе во всех расчетах в качестве функции Ф(г), применялся ограничитель Ван Лира: Ф(г) — (r + |r|)/(1 + |r|) [8].

Дальнейшее вычисление величин на новом временном слое, включая все промежуточные этапы, описано в [7].

Шаг интегрирования по времени т выбирается в процессе счета на каждом временном слое из условий устойчивости для двумерной схемы [5]:

т — K-Tx ' Ty

Tx + Ty

где 0 < К < 1,тх,ту — временные интервалы, за которые волны, образующиеся в задаче о распаде разрыва, достигают грани ячейки по осям х,у соответственно.

2. Результаты расчетов

2.1. Влияние давления внешней среды р0 на процесс кумуляции

Первый расчет проводился для внешнего давления, которое соответствует детонации при нормальных условиях (р0 = 1 бар). Толщина заряда АВ = 4 мм ( = 2 мм). Толщина зазора между зарядом и плоскостью симметрии (ГС на рис. 1) ГЕ = 10 мм. Толщины зазора и заряда не менялись при изменениях давления. Расчет проводился на сетке Ах = Ду = 0.1 мм. Результаты расчета представлены на рис. 2, 3.

Рассмотрим изолинии показателя адиабаты 7 в различные моменты времени во всей расчетной области АОИСГЕ (рис. 2, а). Использование модели среды с эффективным (вычисляемым) показателем адиабаты делает возможным выделение ПВ в окружающем газе. Можно заметить, что процесс формирования кумулятивной струи протекает достаточно быстро, за 4-5 мкс. ПВ практически полностью перекрывают зазор набегающему потоку и распространяются под зарядом в направлении участка ГЕ в виде кумулятивной струи вдоль плоскости симметрии ГС. Часть ПВ, находящаяся между струей и зарядом, смешивается с набегающим потоком.

Рассмотрим изолинии плотности во время формирования кумулятивной струи (рис. 2, б, 3-5 мкс) и установившегося кумулятивного течения (рис. 2, б, 20 мкс). Показана только

Рис. 2. Течение при ро = 1 бар: изолинии показателя адиабаты 7 (а) и изолинии плотности р (б).

3 мкс ------4 мкс ---5 мкс -20 мкс

Рис. 3. Распределение параметров на плоскости симметрии при ро = 1 бар: плотность р (а) и

число Маха М (б).

часть расчетной области АОНСЙ'б'' (см. рис. 1). При истечении ПВ с фронта детонации образуется стандартная ударно-волновая конфигурация: ударная волна (рис. 2, б), контактная граница и следующая за ней волна разрежения. Процесс образования кумулятивной струи возникает в зоне взаимодействия контактной границы АМ с плоскостью симметрии (рис. 2, б, 3 мкс). Далее, в пространство, ограниченное деформирующимся отраженным (точка ) скачком, образованным взаимодействием "головной" ударной волны с плоскостью симметрии, и контактной границей АМ, начинают поступать ПВ (рис. 2, б, 4-5 мкс) с дальнейшим распространением вдоль плоскости симметрии, против направления движения набегающего потока. Деформация отраженного скачка происходит за счет торможения набегающего потока и распространения ПВ в кумулятивной струе. В установившемся кумулятивном течении (рис. 2, б, 20 мкс) можно выделить границу кумулятивной струи, контактную границу (тангенциальный разрыв) АМ и "висячий" скачок МЯ, разворачивающий большую часть ПВ в направлении, параллельном плоскости симметрии и противоположном кумулятивному течению.

На рис. 3 приведены графики распределения по плоскости симметрии значений плот-

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

Как следует из [3], при давлении р0 ^ 100 бар и выше кумуляции ПВ не существует. Поэтому дальнейшие исследования проводились следующим образом: давление р0 повышалось от значения при нормальных условиях (1 бар) до 100 бар, и изучались образующиеся течения.

На рис. 4 представлены результаты расчета при давлении р0 = 10 бар в виде изолиний плотности и показателя адиабаты на момент времени 20 мкс. Толщина заряда и толщина зазора ЕЕ оставались постоянными. По глубине распространения ПВ под зарядом можно заметить, что кумуляция здесь замедлена по сравнению с течением при р0 = 1 бар. Данное явление можно объяснить увеличением противодавления со стороны набегающего потока, связанного с повышением давления р0. В связи с этим изменился и угол наклона контактной границы (точка сместилась по потоку). При сравнении течений по изолиниям показателя адиабаты (см. рис. 2, а и 4, б) следует заметить, что расстояние между изолиниями при р0 = 10 бар увеличено. Это говорит о том, что ПВ при развороте смешиваются с набегающим потоком.

Для изучения влияния повышения давления р0 был выполнен расчет, результаты которого (при р0 = 40 бар) представлены на рис. 5. Толщина заряда и толщина зазора ЕЕ при этом не изменялись. Показаны изолинии плотности и показателя адиабаты на момент времени 15 мкс. Здесь имеет место стационарное течение продуктов взрыва. Положение ударной волны и контактного разрыва со временем не меняется. Продукты взрыва распространяются только в направлении набегающего потока без разворота (без образования кумулятивной струи). По изолиниям показателя адиабаты видно, что в

Рис. 4. Течение при ро = 10 бар: изолинии плотности р (а) и изолинии показателя адиабаты 7 (б).

Рис. 5. Течение при ро = 40 бар: изолинии плотности р (а) и изолинии показателя адиабаты 7 (б).

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

Рис. 6. Потоковые функции: Q(t) (а) и Кт(Ь) (б).

Дальнейшее увеличение давления р0 до значения 100 бар не меняет стационарного характера течения. В зависимости от величины давления изменяется только (в меньшую сторону) угол наклона ударной волны и контактного разрыва к плоскости симметрии ГС.

На рис. 6, а представлены графики массового потока газа через сечение ББ' (см. рис. 1) для течений с различным давлением р0 в зависимости от времени. Для совпадения координат значения массового потока отнесены к соответствующему давлению р0 и представляют собой величину:

н

Я = — (ря ия )(у,

Ро ) о

где ря ,ия — параметры в сечении ББ', И — толщина зазора ГЕ.

На примере графика для давления р0 = 1 бар можно проследить прохождение сечения ударной волной АК (£ = 5 мкс), кумулятивной струей с ПВ (£ = 6 мкс). Подобные графики приведены в [6]. Таким образом, очевидно, что с повышением давления р0 процесс кумуляции ПВ замедляется (графики смещены по времени) и подавляется набегающим потоком (значение массового потока уменьшается) с переходом в стационарное течение с отсутствием кумуляции (см. рис. 5).

2.2. Влияние толщины заряда на процесс кумуляции

Изучение влияния толщины заряда на процесс кумуляции проводилось для р0 = 1 бар при постоянной величине зазора ГЕ = 10 мм. Толщина заряда изменялась от 4 до 8 мм. Для вычисления доли ПВ, попадающей в кумулятивную струю, вычислялся КПД заряда К :

н

Кт = -ТГ~ (рЯия - роио)(у,

рлилИвв 3 о

где ря, ия — параметры в сечении ББ' (см. рис. 1), И — толщина зазора ГЕ, ИВв — толщина заряда АВ. Здесь из массового потока газа через сечение Б Б' вычтена составляющая натекающего массового потока окружающего газа р0 • и0, рл • ил • Ивв — полный поток продуктов взрыва с фронта детонации.

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

2.3. Сравнение с физическим экспериментом

Сравнение результатов численного расчета с физическим экспериментом [4] проводилось для случая синхронной детонации двух плоских параллельных зарядов. Толщина заряда (см. рис. 1) составляла 4.5 мм, величина зазора между зарядами — 6 мм (ЕЕ = 3 мм). Численный расчет проводился на сетке с Ах = А у = 0.1 мм.

Для сравнения определялась зависимость относительной скорости О/О от времени, где и — величина средней (по сечению зазора) скорости ударной волны К К' (см. рис. 2, б, 4), О — скорость детонации. На рис. 7 приведены результаты физического эксперимента [4] и численного расчета. Очевидно, что в начальный момент времени (2-3 мкс) в процессе формирования и начального развития течения кумулятивной струи ПВ скорость ПВ О выше скорости детонации О (О/О > 1). Соответствие на этом этапе результатов численного расчета физическому эксперименту свидетельствует о достаточной точности метода и правильности постановки задачи. Дальнейшее уменьшение скорости О физического эксперимента обуславливается, согласно [4], "формированием сгустка ударно-сжатого газа, что приводит к уменьшению угла разлета ПВ" и прекращению кумуляции. Однако, из результатов численного расчета подобного уменьшения скорости О, показывающего прекращение кумуляции, не видно. Наоборот, результаты численных расчетов при р0 = 1 бар (см. также пп. 1.1, 1.2) показывают усиление кумулятивного течения и, кроме того, не наблюдается каких-либо причин к его ослаблению и тем более прекращению.

Рис. 7. Зависимость относительной скорости и/О от времени.

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

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

Заключение

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

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

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

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

[1] Баум Ф.А., Станюкович К. П., Шехтер Б. И. Физика взрыва. Наука, М., 1975.

[2] БЕЗРУКОВ А. А. Применение квазимонотонных разностных схем второго порядка аппроксимации в моделировании разрывных газодинамических течений. Моделирование в механике, 5(22), №5, 1991, 22-29.

[3] Бывших А. И., Васильев C. M., Кирко B. И. Численное моделирование возникновения газовой кумуляции в зарядах взрывчатых веществ с плоскими и цилиндрическими каналами. ТВТ, 31, №2, 1993, 286-291.

[4] Васильев C.M., Кирко В.И., Пак Н. И. Газодинамическое течение при детонации заряда ВВ в плоском канале. Препринт №491 Ф, Ин-т физики СО АН CCCP, Красноярск, 1988.

[5] Годунов С. К., Забродин А. В., Иванов М.Я., Крайко А. Н., Прокопов Г. П.

Численное решение многомерных задач газовой динамики. Наука, M., 1976.

[6] Лобанов В.Ф., Фадеенко Ю.И. Кумуляция продуктов детонации полого цилиндрического заряда. ФГФ, №1, 1974, 119-124.

[7] Родионов А. В. Повышение порядка аппроксимации схемы С. К. Годунова. ЖВМиМФ, 27, №12, 1987, 1853-1860.

[8] SPEKREIJSE S. Multigrid solution of monotone second-order discretizations of hyperbolic conservation laws. Math. Comput., 49, No. 179, 1987, 135-155.

[9] THOMPSON K. W. Time dependent boundary conditions for hyperbolic systems. J. Comp. Phys, 68, No. 1, 1987, 1-24.

Поступила в редакцию 27 июля 1999 г.