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

| ISSN 2542-0542 Вестник Концерна ВКО «Алмаз - Антей» | № 4, 2017

УДК 623.468

Ю. С. Зуев, Н. И. Карманов

Моделирование процесса перехода горения во взрыв во взрывчатом веществе на основе октогена

Представлены результаты численного моделирования процесса перехода горения во взрыв в образцах взрывчатого вещества LX-10, заключенных в прочную сборку. Расчеты проведены в программе конечно-элементного анализа ANSYS Autodyn с использованием модели горения Powder Burn. Получено удовлетворительное соответствие расчетных и экспериментальных данных. Ключевые слова: горение, взрыв, ANSYS Autodyn, скорость разлета.

Введение

Изделия, содержащие взрывчатые вещества (ВВ), в процессе жизненного цикла могут подвергаться различным аварийным нерегла-ментированным воздействиям (падения, пожар, воздействие пуль и др.). Наиболее частой аварийной ситуацией является пожар. Известно, что при воспламенении ВВ возможны различные исходы: разгерметизация изделия и выгорание ВВ, возникновение конвективного горения, низкоскоростного взрывчатого превращения или детонации [1].

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

Прогнозирование реакции изделия можно осуществлять с помощью современных программ конечно-элементного анализа, например ANSYS Autodyn. В программе ANSYS Autodyn представлена довольно гибкая модель горения Powder Burn, которая подходит для описания процессов горения в широком диапазоне давлений, характерных для переходных режимов взрывчатого превращения (0,002...0,030 Мбар). Структура модели Powder Burn Модель включает в себя уравнения состояния твердого вещества и газообразных продук-

© Зуев Ю. С., Карманов Н. И., 2017

тов разложения, а также скорости выгорания твердого вещества [2]. Считается, что образец твердого вещества состоит из частиц определенной формы соответствующего характерного размера г0.

Уравнение состояния газообразных продуктов разложения может быть задано в виде уравнения JWL. Для описания поведения твердого вещества можно использовать уравнение компактирования, для которого необходимо задать две табличные функции: первая - зависимость давления в твердой фазе от плотности; вторая - зависимость скорости звука от плотности. Скорость выгорания вычисляется по формуле

dF

- = G (1 - F УН (Ps),

(1)

где F - величина выгорания (меняется от нуля до единицы); t - время;

H (P ) - табличная зависимость скорости послойного горения частицы H от давления в газовой фазе Pg.

Здесь также константа G характеризует размер частицы. Например, для сферической частицы G = 3 / r0 (r0 - радиус), для кубической - G = 1 / r0 (r0 - длина ребра). Константы c характеризуют форму частицы. Например, c = 2/3 для сферической или кубической частицы, c = 1/2 - для цилиндрической. Экспериментальные данные Используя модель Powder Burn, попробуем воспроизвести результаты эксперимента [3] по переходу горения во взрыв в высокоплотных, содержащих октоген образцах ВВ LX-10, заключенных в прочную сборку.

Постановка опыта схематично показана на рис. 1. Пять цилиндрических образцов из

ф о о.

I-

Ü о

Ig

та

.

та m о ч

V ^

и о

о

<и

У S

s о о

о сч

OI

<

I

(0 те

s

0 ^

CQ те

1

о.

ф

£

и

V

со

сч ■clin

с?

сч ■clin сч

(П (П

Рис. 1. Постановка опыта [3]: 1 - воспламенитель; 2 - образцы ВВ ЬХ-10; 3 - стальная трубка ЛвтМв1 100; 4 - крышка

ВВ ЬХ-10 (октоген - 95 %, связка витон - 5 %) помещались в стальную трубку с внутренним диаметром 44,93 мм, длиной 203,2 мм и толщиной стенки 2,93 мм. В одном из образцов ВВ располагался воспламенитель.

После задействования воспламенителя произошел взрыв. Зафиксированные скорости разлета трубки и осколков V в плоскости, проходящей через воспламенитель и перпендикулярной оси сборки, показаны на рис. 2. Скорость нарастала в течение примерно 250 мкс (см. рис. 2), что служит подтверждением недетонационного характера взрыва. Существенный рост скорости трубки от 10 м/с до 400...700 м/с

V, м/с

W

10

ю1

10

10

,-1

4

л[/

3 5

и 2

-250 -200 -150 -100

-50

0 t, мкс

(до осколкообразования) происходил в течение 30...41 мкс. Значения скоростей осколков изменялись в диапазоне 800.2100 м/с. Исходные данные для расчетов Цель исследования - описание динамики разлета стальной трубки в опыте [3] с помощью модели Powder Burn. Для этого необходимы исходные данные, которые были заимствованы из работ [3, 4]. Например, в статье [4] приведены значения параметров для уравнения состояния JWL как для твердого ВВ LX-10, так и для газообразных продуктов детонации данного ВВ (табл. 1). Уравнение JWL выглядит следующим образом:

P(P, E) = A

1

wp R tp

- R P0

e p +

о

+B

wp R 2 P

- R

e 2 p + wpE,

о

где P - давление; p - плотность; E - внутренняя энергия; A, B, Rt, R2 - параметры; p0 - исходная плотность ВВ LX-10.

Для описания прочностных свойств ВВ LX-10 выбрана модель прочности по Мизесу с модулем сдвига, равным 0,05 Мбар, и пределом текучести Y0 = 0,002 Мбар.

Параметры для продуктов детонации из табл. 1 можно использовать сразу, так как в модели Powder Burn продукты разложения ВВ также описываются уравнением JWL.

Таблица 1

Параметры для ВВ LX-10 (г/см3)

Параметры Твердое ВВ LX-10 Продукты детонации

А, Мбар 7320 16,689

В, Мбар -0,052654 0,5969

Ri 14,10 5,900

Rl 1,410 2,100

w 0,8867 0,450

E0, Мбар 4,40110-3 0,104

Рис. 2. Скорости разлета стальной трубки и осколков, измеренные по доплеровской методике (PDV) и с помощью радаров (Radar) соответственно [3]:

1 - PDV3; 2 - PDV2; 3 - PDV1; 4 - Radar 1; 5 - Radar2; 6 - Radar3

В модели Powder Burn для описания поведения твердого вещества используется уравнение компактирования. Изначально оно было предназначено для моделирования поведения пористых веществ, таких как песок или порох

| ISSN 2542-0542 Вестник Концерна ВКО «Алмаз - Антей» | № 4, 2017

насыпной плотности. Известно, что вещества в насыпном состоянии при нагрузке существенно уплотняются, а при снятии нагрузки их плотность падает, но не до исходного состояния, а до некоторой величины, большей начальной плотности. Необходимо приспособить уравнение компактирования для рассматриваемой задачи, в которой используются изначально высокоплотные малопористые образцы ВВ. Для нахождения параметров уравнения ком-пактирования используем параметры уравнения JWL для твердого ВВ ЬХ-10 (см. табл. 1).

Для уравнения компактирования необходимо задать две табличные функции: функцию Р(р), определяющую зависимость давления в твердой фазе от плотности в процессе нагру-жения; функцию с (р(4), описывающую разгрузку ранее нагруженного до определенного давления вещества и определяющую зависимость скорости звука от плотности разгруженного вещества. Для определения значений табличных функций, фигурирующих в уравнении компактирования, следует построить график уравнения JWL для случая ударно-волнового нагружения с параметрами из табл. 1 (рис. 3).

Значения табличной функции Р(р) можно выбрать сразу. Кривая на рис. 3 получена из уравнения JWL для твердого ВВ, точки соответствуют табличной функции Р(р) . Возьмем

Рис. 3. Определение значений табличных функций Р(р) и с (р^>) :

-- JWL; 1 - fa; P); 2 - (р2; P2); 3 - (р3; P3);

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

Таблица 2

Значения табличной функции

№ P, Мбар р, г/см3 № P, Мбар р, г/см3

1 0 1,865 6 0,043 2,15

2 2,84840-3 1,900 7 0,094 2,29

3 6,88840-3 1,940 8 0,169 2,40

4 0,012 1,980 9 0,248 2,48

5 0,019 2,030 10 0,539 2,66

Для заполнения табличной функции с ) можно поступить так, как показано на рис. 3, где кривая соответствует уравнению JWL для твердого ВВ, точки - табличной функции Р(р) (см. табл. 2). Линия, проходящая, например, через точки (р4; Р4) и (р5; Р5), пересечет ось абсцисс в точке (р4°; Р = 0). Тангенс угла наклона этой прямой к оси абсцисс равен квадрату скорости звука с4. Вдоль этой линии будет происходить разгрузка ВВ, изначально нагруженного до состояния (р5; Р5). В результате получена пара значений ; с4). Аналогично можно получить остальные пары значений. В табл. 3 представлены значения табличной функции с (р(4), определяющей разгрузку ВВ.

Таблица 3 Значения табличной функции с(р(с))

№ с, см/мкс р(с), г/см3 № с, см/мкс р(с), г/см3

1 0,2870 1,8650 6 0,6036 2,0313

2 0,3178 1,8711 7 0,8257 2,1514

3 0,3575 1,8854 8 0,9937 2,2282

4 0,3742 1,8936 9 1,2715 2,3259

5 0,4472 1,9343 10 1,4718 2,4105

Для стали ЛетЫе1 100 использовано уравнение состояния Грюнайзена:

Р(р,Е) = Рн + ГоРо(Е - Ен), р _ РоСсУ(1 + Ю „ Рн (1+ Д) Д

[1 - (S -

^H

2ро 1 + Д

Ро

4 - (р4; Р4); 5 - (р5; P5); 6 - (р<?; с4) ; 7 - (р?; c5)

где Го = 1,84 ;

р0 = 7,923 г/см3;

ф о о.

I-

Ü о

I-

та

.

та m о

ч

ф

и о

о

<и у

s

s о о

о сч

01

<

I

(0 те

s

о ^

со те г о. ф

£

и ф

со

сч ■ci-io

с?

сч ■ci-io сч

(П (П

c0 = 0,4534 см/мкс; S = 1,5 [3].

Также введем модель прочности Штейн-берга - Гинена:

7 = Y0 (1 + ё'),

Gm - G0

1 + b-

n

1/3

h (T - To)

n

=_p_ Po

d-F = / (1 - F )b dt

f V

P-1 - a

VPc

+

0< F < Fg

(2)

+ G1(1 - F)cFdPy + G2 (1 - F)eFgPz

0< F < FG

Fg 2min <F <1

Здесь p - текущая плотность;

p0 - начальная плотность ВВ;

I, Gj, G2, a, b, c, d, e, g, x, y, z,

F , FG1 ,

igma^' Glmax '

F

G 2m

параметры.

Эта модель присутствует в расчетных программах ANSYS Autodyn и LS-DYNA. Для ВВ LX-10 целесообразно использовать ее в диапазонах давлений 0,015.0,5 Мбар. Для применения модели Powder Burn в данном диапазоне давлений необходимо подогнать параметры модели под параметры модели Ignition and Growth, значения которых приведены в табл. 4 [4].

Таблица 4

Кинетические параметры LX-10 для модели

Ignition and Growth [4]

a = 0 G1 = 3,1 G2 = 400 I = 7,43 • 1011

b = 0,667 y = 2 z = 2 F = 0 3 g max ?

x = 20 c = 0,667 e = 0,333 FG1 max = 0,5

d = 0,333 g = 1 FG2 min = 0,5

где Y - динамический предел текучести;

ёp - эквивалентная пластическая деформация;

GCR - динамический модуль сдвига;

T - температура.

Параметры для стали AerMet 100 [3] следующие. При температуре T0 = 300 K и давлении P = 0 начальный модуль сдвига G0 = = 0,748 Мбар. Начальный предел текучести Y0 = = 0,01016 Мбар при температуре T0 = 300 K, давлении P = 0. Максимальный предел текучести Ymax = 0,03 Мбар, в = 2, b = 0, n = 0,5, h = 0.

Известно, что для количественного описания процесса перехода относительно сильной ударной волны в волну детонации в высокоплотных ВВ используется модель Ignition and Growth [4], в которой уравнение выгорания (1) можно записать так:

Скорость горения H(Pg) в модели Powder Burn задается табличным способом, значит, можно для каждого значения величины давления Pg задать такую скорость сгорания H, которая позволит описать в среднем поведение кривой dF/dt, задаваемой уравнением (2). Для этого при конкретном значении давления найдем среднее значение функции dF / dt на интервале значений выгораний F = 0-1. Среднее значение функции dF/dt вычисляется по формуле

у v

dF_ dt

=J

I (1 - F )b

1 - a

vFo ,

+

0< F < Fig

G(1 -F)cFd Py G2(1 - F )eFg P

+

0< F < FG1

n < F <1

dF.

Примем, что в уравнении (1) величина G = 1, а c = 0,001. Тогда

G (1 - F )cH (Pg)

h (p)=f

Значение параметра c = 0,001 соответствует, например, кубическим частицам, горящим только на одной из шести поверхностей. Коэффициент G принят равным единице для удобства определения Н (р ). Его реальные значения для ВВ LX-10 лежат в диапазоне 100.200 1/см (размер частицы ВВ г0 ~ 1 / G = = 0,005...0,01 см). Таким образом, информация о величине G фактически будет содержаться в значениях табличной функции

Н () = (dF / dtТеперь можно выбрать необходимое число характерных значений давлений и для каждого этого значения определить

величину Н () = (dF / dt) (табл. 5).

Затем выберем значения функции Н(р) для давлений р < 0,01 Мбар. Скорости послойного горения, измеренные при Р„ = 0,03 Мбар

| ISSN 2542-0542 Вестник Концерна ВКО «Алмаз - Антей» | № 4, 2017

Таблица 5

Значения табличной функции Н(Рг) для высоких давлений (Р& > 0,015 Мбар)

№ 1 2 3 4 5 6

Pg, Мбар 0,015 0,030 0,050 0,080 0,120 0,450

H, мкс-1 0,019 0,077 0,216 0,613 2,792 1301

Таблица 6

Значения табличной функции Н(Рг) для низких давлений (Pg < 0,01 Мбар)

№ 1 2 3 4 5

Pg, Мбар 10-6 0,0001 0,0016 0,0040 0,007

H, мкс-1 10-8 10-4 0,0020 0,0055 0,010

и Р8 = 0,12 Мбар, равны V = 7 • 10-4 см/мкс и V = 0,03 см/мкс [5]. Если значения скоростей умножить на 100, то получим 0,07 и 3 соответственно, что близко к значениям функции Н (0,03 Мбар) = 0,077 мкс-1 и Н (0,12 Мбар) = 2,792 мкс-1 (см. табл. 5).

Формально получается, что в уравнении (1) должно быть О ~ 100 . Такое значение соответствует размеру частицы ВВ г0 ~ 1 / О = = 0,01 см и реальным размерам частиц окто-гена. Примем, что О = 100 и для давлений Р < 0,01 Мбар. Теперь умножим скорости 2-10-5 см/мкс, 5,5 • 105 см/мкс и Ы0-4 см/мкс [5] послойного горения при значениях Р, равных 0,0016 Мбар, 0,004 Мбар и 0,007 Мбар на 100. При этом примем константу О = 1. Получим значения функции Н (р ) для давлений Р8 < 0,01 Мбар (табл. 6).

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

R, см

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

Расчетная геометрия показана на рис. 4. Размер сетки - 0,1 см. Расчеты проведены в осе-симметричной ALE постановке в программе ANSYS Autodyn (лицензия № 774-2013-ША от 01.10.2013). Воспламенение начиналось в нескольких элементах, в которых в течение 30 мкс задавалось начальное давление 0,002 Мбар. В расчетах не учтено влияние воспламенителя, а также каналов в образцах ВВ, в которых прокладывались провода для воспламенителя.

В результате расчета с использованием функции H (Pg) (см. табл. 5, 6) и G = 1 получено, что скорость разлета трубки нарастала от 10 до 100 м/с в течение 16 мкс, а от 100 до 700 м/с - 15 мкс. При этом переход в детонацию произошел достаточно рано. В эксперименте скорость нарастала от 10 до 100 м/с и от 100 м/с до 700 м/с в течение 19 мкс. Если G = 0,7, то скорость разлета трубки нарастает от 10 до 100 м/с в течение 26 мкс, а от 100 до

-1

-з

-5

-7

L

-13 -11

-9

11 z, см

-7 -5 -3 -1 1 3 5 7 9 Рис. 4. Постановка расчетов: 1 - крышкаЛетЫе1 100; 2 - трубкаЛетЫе1 100; 3 - ВВ ЬХ-10; 4 - область воспламенения

V

о о.

I-

ü о

Ii та

.

та m о

ч

ф

^

и о

о

<и

У S

s о о

Рис. 5. Результат расчета в момент достижения скорости V = 700 м/с

о см

01

<

I

(0 те

s

0 ^

CQ те

1 о. ф

£

и ф

CQ

СМ ■clin

с?

см ■ci-io см

(П (П

700 м/с - 21 мкс. При этом детонация возникает ближе к крышке. Если коэффициент G = 0,5, то скорость разлета трубки нарастает от 10 до 100 м/с в течение 30 мкс, а от 100 до 700 м/с - 36 мкс. Перехода в детонацию при этом не происходит.

Расчеты при G = 0,7 наиболее хорошо согласуются с экспериментом. Результаты этого расчета показаны на рис. 5. Сравнение расчетных и экспериментальных данных приведено на рис. 6.

Расчетное и экспериментальное перемещение стенки трубки ~1 см. Максимальная расчетная скорость трубки ~2100 м/с, что соот-

V, м/с

W

10

ю1

10

10

-1

A K

n V

3 J

14 \ X., Ж /

>f1 f

-250

-200

-150

-100

-50

0 t, мкс

Рис. 6. Сравнение расчетных и экспериментальных

скоростей движения трубки: 1 - PDV3; 2 - PDV2; 3 - PDV1; 4 - расчетные данные; 5 - Radar 1; 6 - Radar2; 7 - Radar3

ветствует максимальным скоростям осколков, зафиксированным в опыте. В работе [3] отмечено, что разрушение трубки начинается при 15%-ном радиальном растяжении. В расчете при такой деформации скорость трубки достигла 500 м/с, что соответствует экспериментальным скоростям 400...700 м/с, при которых начинается осколкообразование. Заключение

С помощью программы конечно-элементного анализа ANSYS Autodyn и встроенной в нее модели горения Powder Burn удалось удовлетворительно описать результаты эксперимента по переходу горения во взрыв в образцах ВВ LX-10, заключенных в прочную сборку [3]. Скорость выгорания в модели Powder Burn задается табличным способом, поэтому появляется возможность описать энерговыделение для широкого диапазона давлений (Pg = 0,002... 0,5 Мбар), где может реализоваться переход горения в детонацию.

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

1. Переход горения конденсированных систем во взрыв / А. Ф. Беляев, В. К. Боболев, А. И. Ко-ротков и др.; под ред. П. Ф. Похила. М.: Наука, 1973. 293 с.

2. AtwoodА., Friis Е. K., Moxnes J. F. А Mathematical model for combustion of energetic powder materials // 34th International Annual Conference of ICT. Karlsruhe, Federal Republic of Germany, 24-27 June, 2003. Pp. 1-14.

I ISSN 2542-0542 Вестник Концерна ВКО «Алмаз - Антей» I № 4, 2017

3. McClelland K. A., Maienschein J. L., Howard W. M., Nichols A. L., de Haven M. R., Strand O. T. Measurement and ALE3D simulation of violence in a deflagration experiment with LX-10 and Aermet-100 alloy // 13 th International Detonation Symposium. Norfolk, VA, United States 23-28 July, 2006. Pp. 1-12.

4. Tarver C. M., Urtiew P. A., Chidester S. K., Green L. G. Shock compression and initiation of

LX-10 // Propellants, Explosive, Pyrotechnics. 1993. Vol. 18. Pp. 117-127. 5. Nichols A. L., Tarver С. M., McGuire Е. М. ALE3D statistical hot spot model results for LX-17 // Proceedings of Conference of the American Physical Society Topical Conference on Shock Compression of Condensed Matter. Portland, Oregon, July 20-25, 2003. Pp. 1-6.

Поступила 30.11.17

Зуев Юрий Семёнович - кандидат технических наук, старший научный сотрудник, заведующий кафедрой технической механики Снежинского физико-технического института - филиала государственного бюджетного образовательного учреждения высшего образования «Национальный исследовательский ядерный университет «МИФИ», г. Снежинск.

Область научных интересов: динамика и прочность машин, летательных аппаратов.

Карманов Николай Иванович - аспирант Снежинского физико-технического института - филиала государственного бюджетного образовательного учреждения высшего образования «Национальный исследовательский ядерный университет «МИФИ», г. Снежинск.

Область научных интересов: моделирование процессов перехода горения во взрыв во взрывчатом веществе.

Simulation of combustion-to-explosion transition processes in Octogen-based explosive

The study presents the results of numerical simulation of combustion-to-explosion transition processes in the samples of the explosive LX-10, enclosed in a solid assembly. The calculations were carried out in the ANSYS Autodyn finite element analysis program using the Powder Burn combustion model. The satisfactory agreement between the calculated and experimental data was obtained. Keywords: combustion, explosion, ANSYS Autodyn, separation velocity.

Zuev Yuriy Semenovich - Candidate of Engineering Sciences, Senior Research Fellow, Head of the Department of Technical Mechanics, Snezhinsk Physics and Technology Institute of the National Research Nuclear University "MEPhl", Snezhinsk.

Science research interests: aircraft and machines dynamics and strength.

Karmanov Nikolay Ivanovich - post-graduate student, Snezhinsk Physics and Technology Institute of the National Research Nuclear University "MEPhl", Snezhinsk.

Science research interests: simulation of combustion-to-explosion transition processes in an explosive.

<u о о.

I-

ü о

Ii те

.

<я m о

ч

ф

и о

о

<и

У

s

s о о