Определение параметров уравнения состояния JWL продуктов детонации по данным цилиндр-теста: метод с понижением размерности

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

Экспериментальная методика цилиндр-тест [1, 2] является наиболее популярным методом оценки рабочих характеристик ВВ. В стандартной постановке используется медный цилиндр, заполненный ВВ, которое инициируется с одного торца цилиндра в начальный момент времени. В ходе эксперимента получают информацию о смещении наружной границы трубки в зависимости от времени, которую используют для расчета кинетической энергии ПД, а также для определения параметров уравнения состояния (УРС) продуктов детонации, из которых наиболее популярным является УРС Jones–Wilkins–Lee (JWL) [2]:

( p (V, E) = A 1 к

e ^R 1 ^V

+ B 1 к

R v J

e - R 2 V + ^{— E} V

и соответствующее ему уравнение изэнтропы p ( V )

p(V ) = Ae - ^R1V + Be - ^{R 2} V + CV ^{- (1 + —}) ,

где p - давление продуктов детонации, V = р ₀/ р - степень расширения продуктов детонации, а A, B , C , R 1 , R 2 , m - параметры, подбираемые для каждого взрывчатого вещества.

Уравнение изменения энергии на изэнтропе для данного уравнения состояния имеет вид

E ( V ) = j pdV = RA- e -^R - ^V + Be - ^R 2 V + CV - ^m .                        (3)

Классическая методика [2] предполагает итерационную подгонку параметров уравнения состояния JWL в гидрокоде до совпадения численного и экспериментального расширения цилиндрической трубки, которая получила развитие, например, в работах [3, 4]. Существуют аналитические подходы, основанные на использовании соотношения энергетического баланса, которое связывает энергию взрыва и совершенную работу со степенью расширения продуктов ВВ [5–11]. Также известны различные методы прямого расчета неизвестных параметров уравнения состояния JWL по экспериментальным данным [12–15]. Однако большинство из них требуют оптимизации в шестимерном пространстве параметров A, B, C, R1, R2, m, что сопряжено с вычислитель- ными трудностями и неоднозначностью решения.

В настоящей работе предлагается метод, сводящий задачу к поиску всего трёх параметров R ₁, R ₂, т с последующим линейным расчётом A , B , C . Это достигается за счёт использования условий в точке Чепмена–Жуге (CJ) и данных об энергии Гарни на трёх фиксированных степенях расширения. Дополнительно вводятся физические ограничения, обеспечивающие доминирование соответствующих членов уравнения в различных диапазонах давлений.

В качестве необходимых данных для метода требуются:

• 1.    Траектория разлета стенки, полученная из эксперимента и представляющая собой набор дискретных значений.

• 2.    Параметры материала и геометрии цилиндра: плотность, упругие константы, внешний и внутренний радиусы. Отметим, что в большинстве экспериментов в качестве материала используется медь, для которой достаточно точно определены требуемые параметры.

• 3.    Параметры самого ВВ, которые должны быть определены из дополнительных экспериментов: плотность при нормальных условиях, удельная энергия взрыва, скорость детонационной волны и один из параметров в точке Чепмена–Жуге.

Данный метод основан на следующих допущениях:

• 1.    Материал стенки цилиндра несжимаемый, а его предел текучести постоянный.

• 2.    Детонация является идеальной.

• 3.    Взрывчатое вещество мгновенно превращается в продукты взрыва.

• 4.    Продукты взрыва невязкие и расширяются изоэнтропически.

• 5.    Течение продуктов взрыва симметричное относительно оси цилиндра.

Экспериментальные данные и их предварительная обработка

Результатом эксперимента цилиндр-тест является траектория движения стенки цилиндра, которая представляется в дискретном виде как набор точек

[ t i ,( А r 2 ) i ] , i = 1,2,..., n ,                                                      (4)

где n - общее количество точек, измеренных в эксперименте, Аr2 = r2 - r20 - смещение внешней поверхности цилиндра, r20 – начальный внешний радиус.

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

Первая функция для аппроксимации имеет вид [9]

r ₂

u ₀

- a ) [ 1 - exp( - b ( t - 1 0 ) ] ,

где a , b , t ₀, u ₀ – параметры, подлежащие оптимизации.

Вторая зависимость имеет вид [8]

r2(t) =£ ai [bt -(- - e ')], i=1

где a i , b i ( i = 1,2) - параметры, подлежащие оптимизации.

Третья зависимость имеет вид [15]

• Г 2 ( t ) =              ^        ,

^v g(0) + g ( t ) a ₀

где a ₀ - начальное ускорение цилиндра, v_x - асимптотическая радиальная скорость цилиндра, а функция g ( t ) определяется как g ( t ) = (1 + 1 ) ^а - 1.

Таким образом, в данном случае требуется определить параметры a ₀, v_x , а .

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

Выбранная аппроксимация A r₂ = F ( t ) позволяет вычислить скорости и ускорения в эйлеровой системе координат r 2 (t ), r ^ (t ), r'( t ), r 2 (t ), r^t ), r ^ (t ).

Переход к лагранжевым координатам выполняется с учётом угла наклона образующей цилиндра 0 (рис. 1). Из кинематических соотношений [8] имеем:

va = ~c(t), tan 0 = —, v = 2D sin 0, a = aa cos3 0, a dt            D             2 a где D – скорость детонации, Fc (t) – смещение центральной поверхности цилиндра, определяемое из условия несжимаемости материала, для определения которого используется соотношение 1

r2 - rc = rc - r1 = 2 (r20 - r10 ), где r10, r20 – начальные внутренний и внешний радиусы соответственно, rc – радиус центра стенки.

Рис. 1. Схема расширения цилиндрической оболочки и кинематические параметры

Энергетический баланс и определение внутренней энергии ПД

Система, состоящая из заряда ВВ и цилиндрической трубки, удовлетворяет закону сохранения энергии. За время At через сечение проходит масса p0 DA At, и энергия распределяется сле- дующим образом [8]:

I                         2

= puA A t ,

Ро DA A t U - Q + - + Tkm + Adef \        2

где U – удельная внутренняя энергия ПД, Q – удельная теплота детонации, u – массовая ско- рость ПД в осевом направлении, Tkin – кинетическая энергия вследствие радиального движения цилиндра и газов (энергия Гарни) на единицу массы ВВ, Adef – работа деформации цилиндра на единицу массы заряда ВВ.

Переходя к удельным величинам на единицу объёма ВВ, получим

E - E 0 +     + E kn + W def = pU ,                            (5)

2D где E = p0U и E0 = p0Q - внутренняя энергия и теплота детонации единицы объема заряда ВВ, Ekin = Р0Tn и Wdef = Р0Adef - энергия Гарни и работа деформации на единицу объема.

Массовая скорость ПД определяется через степень расширения по формуле u = D (1 - V (гю /r1 )2) , где степень объемного расширения ПД определяется по формуле

r 1 ²

r 2 cos 9

.

Удельную энергию Гарни рассчитаем по формуле [ 16 ]

E kin

P

r ₂

In r ²

к r io 7

⁴ к r i 7

cos 9

( r ₂ ' )²,

где p m - плотность материала цилиндра.

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

Работа деформации рассчитывается по соотношению

W

Wdef

= ^a f

r 20

- 1

к r io 7

In r ' r '

r 10 ⁺ r 20

где a_f - предел текучести материала цилиндра.

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

Сведение задачи определения параметров уравнения состояния JWL к трёхмерной оптимизации

Подбор оптимальных параметров проводится исходя из условия наилучшего совпадения изэнтропы уравнения состояния E_JWL (3) с зависимостью для внутренней энергии ПД E (5), т. е. в ходе решения следующей задачи оптимизации:

II E jwl ( A , B , C , R i , R 2 , ю ) - E ||² ^ min.

При поиске оптимальных параметров A, B, C, R1, R2, ю необходимо учитывать, что уравне- ние состояния JWL должно одновременно удовлетворять трём условиям в точке Чепмена–Жуге, следующим из законов сохранения и термодинамики:

• 1.    Давление в точке Чепмена–Жуге должно быть равно экспериментально определенному Ae - ^R V^CJ + Be - ^{R 2 V cj} + CV J ^ю ) = p_CJ .                               (6)

• 2.    Уравнение линии Рэлея в любом состоянии детонации имеет вид:

  P

  F V )

  
  

  = - ^p o D ² .

  
  

Наклон изоэнтропы (2) в точке Чепмена-Жуге должен быть равен наклону линии Рэлея dP k d V J S

dp

. dV

Отсюда следует второе условие

R 1 Ae - ^r ' ^{V cj} + R ₂ Be - ^{R 2} V CJ + ( ю + 1 ) CV _cj ^{- ( ш + 2 )} = p ₀ D ².                       (7)

3. Давление, объем и внутренняя энергия должны соответствовать уравнению состояния

JWL в точке Чепмена–Жуге

A 1

ю

к    R 1 V cj )

e ^R ' V^CJ

(

+ B 1

ю

к    ^R 2 ^V CJ )

e ^{- R 2} V CJ + E cj = P cj ,

где внутренняя энергия определяется из соотношения Гюгонио

E CJ = ^E 0 ⁺ 2 ^pCJ ^{( 1} - V CJ ) .

Кроме того, из (3) для любых степеней расширения Vj должно выполняться равенство меж- ду энергией, вычисленной через (5), и аналитическим выражением (3). В классическом подходе [2] и работах [17, 18] обосновано использование трех значений Vj = 2,4, 4,4, 7,0, соответствую- щих характерным стадиям расширения ПВ (от малых до больших объёмов). Таким образом, имеем три дополнительных уравнения:

A-e - RV + Be ^- R ²V^j + CV = EV j ), j = 1,2,3.                      (9)

R 1        R ₂        a

Уравнения (6)–(9) образуют систему шести нелинейных уравнений относительно шести не известных A,B,C,R1,R2,a. Ключевая идея предлагаемого метода заключается в том, что если временно зафиксировать значения R1, R 2,ю, то уравнения (6)-(8) становятся линейными относительно A, B, C . Действительно, перепишем их в виде:

e ^{- R}V^CJ • A + e — ^{R 2 Vcj} • B + V J^ю ) • C = p_CJ ,

R1 e~RV • A + R2e-R2VCJ • B + (1 + a)VtJRa) • C = p0D2, e - R1Vj      -RR2Vj      V 7ю

• A + —— • B + - j- • C = E (V j ), j = 1,2,3.

R 1        R 2        a

Это – система из пяти линейных уравнений (одно из (9) для трёх Vj даёт три уравнения) с тремя неизвестными A, B, C . Система переопределена, на практике можно использовать любые три уравнения (например, (6), (7) и одно из (9)) для точного определения A, B, C , а остальные два уравнения служат для контроля точности. Такой подход гарантирует выполнение условий Чепмена–Жуге и наилучшее воспроизведение энергии Гарни.

Таким образом, процесс подбора параметров состоит в использовании метода наименьших квадратов в пространстве поиска R 1 , R ₂, ю с невязками в виде E - E(V ) JWL в точках V = 2,4, 4,0, 7,0 либо в диапазоне степени расширения продуктов детонации V е [2,7] с использованием большего числа точек. Минимизация невязки осуществляется алгоритмом Левенберга– Марквардта.

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

• -    все параметры уравнения состояния ( A , B , C , R 1 , R ₂, ю ) должны быть положительными для обеспечения положительности давления при любых степенях расширения. Кроме того, изоэнтропа должна быть монотонно убывающей функцией, что также выполняется при положительных

параметрах автоматически;

– первый член уравнения состояния должен преобладать над вторым в точке Чепмена–Жуге (высокие давления). В окрестности точки Чепмена–Жуге первый экспоненциальный член должен преобладать, чтобы правильно воспроизводить крутой рост давления при сжатии (рис. 2, кривая 1). Таким образом, условие имеет вид

Ae ^{- r} ' ^Vcj >  Be ^{- R 2 Vcj} ;

• -    второй член уравнения состояния должен преобладать над первым при V >  2 (средние степени расширения). На этом участке происходит основное ускорение оболочки, и второй экспоненциальный член определяет форму изэнтропы (рис. 2, кривая 2)

  —е ² R R 1

  
  

B

< —e ^{2 R 2}

;

R 2

– третий (степенной) член уравнения состояния должен преобладать над экспоненциальными при больших степенях расширениях ( V >  7 ), что обеспечивается выполнением условия (рис. 2, кривая 3)

Ae ^{- 7 R} + Be ^{- 7 R} 2 <  ^C 7 ^{- ю} ;

R 1       R 2        a

Рис. 2. Вклад слагаемых УРС JWL при различных степенях расширения:

1 – первый член, 2 – второй член, 3 – третий член

Указанные ограничения реализуются через добавление штрафа к целевой функции при невыполнении условий.

Таким образом, предложенный метод обладает рядом преимуществ:

– понижение размерности с шести до трёх параметров резко сокращает вычислительные затраты и уменьшает риск попадания в локальные минимумы;

– линейное определение A , B , C гарантирует единственность решения для заданных R ₁ , R ₂ , ω и исключает необходимость в итерациях по этим коэффициентам;

– физические ограничения отсекают нефизичные комбинации параметров;

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

Результаты определения параметров уравнений состояния для эталонных ВВ

Разработанная методика была применена к четырём хорошо изученным ВВ: ТЭН (PETN), Состав В (Comp B), нитрометан (Nitromethane) и циклотол (Cyclotol). Исходные данные о плотности, скорости детонации, давлении в точке Чепмена–Жуге и теплоте детонации взяты из справочника [19].

В табл. 1 приведены параметры уравнения состояния JWL, полученные предлагаемым методом, в сравнении с литературными данными.

Таблица 1

Параметры уравнения состояния JWL

ВВ	Источник	A(ГПа)	B(ГПа)	C(ГПа)	R 1	R 2	ω
ТЭН (PETN)	[2]	796,53	19,241	0,665	4,8	1,2	0,25
	[12]	670,36	10,400	1,564	4,44	1,19	0,35
	Данная работа	737,656	15,397	0,895	4,634	1,169	0,276
Состав В (CompB)	[2]	524,229	7,6783	1,0818	4,2	1,1	0,34
	[8]	497,08	3,4246	1,1260	4,0489	0,7833	0,346
	Данная работа	502,987	4,131	0,976	4,07	0,818	0,32
Нитрометан (Nitromethane)	[2]	209,25	5,6895	0,7704	4,4	1,2	0,3
	[12]	203,532	3,605	0,984	4,286	1,0839	0,35
	Данная работа	224,003	7,473	0,889	4,568	1,387	0,326
Циклотол (Cyclotol)	[2]	603,414	9,9236	1,0753	4,7	1,1	0,35
	[12]	560,038	5,131	1,361	4,12	0,995	0,35
	Данная работа	595,028	10,182	1,523	4,304	1,23	0,412

На рис. 3 представлены графики внутренней энергии E ( V ) , вычисленной по энергетическому балансу (5) и аппроксимированной уравнением (3) с найденными параметрами. Для всех ВВ наблюдается хорошее совпадение во всём диапазоне расширения V е [1,7]. Наибольшее расхождение не превышает 5 %.

в)

г)

б)

Рис. 3. Зависимость внутренней энергии от степени расширения ПД для ТЭНа ( а ), состава В ( б ), нитрометана ( в ) и циклотола ( г ): 1 – расчет по формуле (2), 2 – расчет по формуле (4)

Численная верификация

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

Постановка задачи: В начальный момент времени медная трубка заполнена ВВ. Детонация инициируется с левого торца. Геометрические размеры трубки [2]: внутренний радиус R _in = 1,27 см, внешний радиус R out = 1,53 см, длина трубки L = 30,5 см.

Математическое моделирование проводилось в ПК ANSYS Autodyn с использованием ла-гранжево-эйлерового метода на базе метода конечных элементов (ALE). Размер расчетной области - 560 х 112 мм, размер ячейки: 0,35 х 0,35 мм, ячеек в области - 1600 х 320. Продукты детонации описывались уравнением (1) с параметрами из табл. 1, а медь – уравнением состояния в форме Ми-Грюнайзена c начальной плотностью р ₀ = 8942 кг/м³ и объёмным модулем упругости K = 135,4 ГПа. Для моделирования упругопластических свойств меди использовалась модель

Шестаковская Е.С., Определение параметров уравнения состояния JWL Стариков Я.Е. продуктов детонации по данным цилиндр-теста … пластичности Джонсона–Кука с модулем сдвига G = 45 ГПа и пределом текучести σ = 89,6 МПа.

На рис. 4 приведены рассчитанные зависимости Δr2 (t) для четырёх ВВ в сравнении с экспе- риментом [2] и с моделированием с использованием параметров из [2, 8, 12].

б)

в) г)

Рис. 4. Зависимость изменения внешнего радиуса трубки от времени для ТЭНа (а), Состава В (б), нитрометана (в) и циклотола (г): 1 – эксперимент [2], 2 – расчет с коэффициентами [2], 3 – расчет с коэффициентами данной работы, 4 – расчет с коэффициентами [8] для Состава B и – расчет с коэффициентами [12] для остальных ВВ

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

Таблица 2

Сравнение среднеквадратичных отклонений (единицы измерения – см)

	[2]	Данная работа	[8]	[12]
ТЭН	0,4141	0,6551	–	1,0227
Состав B	0,6859	0,7898	0,8336	–
Нитрометан	0,3748	0,3674	–	0,4633
Циклотол	1,4379	0,4161	–	0,7433

Из табл. 2 видно, что для ТЭНа и Состава В среднеквадратичное отклонение предложенного метода превышает значения, полученные в классической работе [2]. Это связано с принципиальным различием в постановке задачи: классический подход минимизируют кинематическую невязку – разницу между расчетной и экспериментальной кривой расширения цилиндрической трубки. Предложенный же метод минимизирует энергетическую невязку – отклонение значений внутренней энергии, рассчитанных по уравнению изэнтропы (3) от значений, полученных из уравнения энергетического баланса (5). Другими словами, наш метод находит параметры JWL, которые обеспечивают строгое выполнение закона сохранения энергии на всем протяжении расширения, тогда как кинематическая точность выступает в роли вторичного критерия. Для задач метания оболочек и расчета работы расширения энергетический баланс является критически важным, поэтому небольшое ухудшение среднеквадратичного отклонения имеет меньшее значение по сравнению с физической согласованностью параметров. Для нитрометана и циклотола предложенный метод показал более высокую точность по среднеквадратичному отклонению, что свидетельствует о хорошей приспособляемости алгоритма к разным типам ВВ.

Таким образом, сравнение с литературными параметрами показывает, что различия в отдельных коэффициентах могут быть значительными (например, для Состава B параметр B отличается почти в два раза), однако итоговые зависимости Δ r ₂( t ) и E ( V ) близки. Это свидетельствует о корреляции параметров JWL и о том, что разные наборы могут описывать одну и ту же изэнтропу с приемлемой точностью. Тем не менее для численного моделирования конкретных задач (например, метания оболочек) предпочтительнее использовать параметры, полученные с учётом энергетического баланса, как в данной работе.

Заключение

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

• •    автоматический выбор наилучшей аппроксимирующей функции для экспериментальной кривой расширения цилиндрической стенки;

• •    учёт кинематических поправок при переходе к лагранжевым координатам;

• •    вычисление внутренней энергии ПД из уравнения энергетического баланса;

• •    решение системы из пяти линейных уравнений для A , B , C при фиксированных R ₁ , R ₂ , ω ;

• •    наложение физических ограничений на параметры.

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