CC BY
339
электронное научно-техническое издание
НАУКА и ОБРАЗОВАНИЕ
Эя № ФС 77 - 305БЭ. Государствен над регистрация №0421100025.155Н 1994-0405_
Математическое моделирование функционирования взрывных устройств
77-30569/334177
# 02, февраль 2012 Колпаков В. И.
УДК 519.63:532:539.5
МГТУ им. Н.Э. Баумана Kolpakov@sm.bmstu.ru
Взрывные или импульсные устройства в настоящее время широко используются в боеприпасах и средствах поражения, в ракетной технике, в работах по интенсификации скважин при газо- и нефтедобыче. К ним, например, можно отнести фугасные, осколочные и кумулятивные заряды разного назначения. Однако, несмотря на то, что экспериментальные исследования играют ключевую роль в изучении импульсных устройств и технологий, использующих в своей основе взрывчатые вещества (ВВ), без глубокого теоретического анализа, как правило, не удается достичь требуемого результата. Кроме того, современные условия диктуют необходимость существенного сокращения количества испытаний при их отработке. В этой связи существенно возросли значение и практическая ценность исследований, проводимых на основе численных методов механики сплошной среды, которые, в свою очередь, предъявляют повышенные требования не только к качеству физико-математических моделей, но и уровню разработанных на их основе алгоритмов расчета. В рассмотренном ключе в настоящей работе излагается методика численного расчета, хорошо зарекомендовавшая себя при моделировании функционирования взрывных устройств различного назначения в течение последних пятнадцати лет.
1. ФИЗИКО-МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ
Изучение процессов метания оболочки осколочных макетов или обжатия кумулятивной облицовки (КО) продуктами детонации (ПД) с последующим формированием кумулятивной струи (КС) или поражающего элемента (ПЭ) как для осесимметричных, так и для удлиненных кумулятивных зарядов (КЗ) целесообразно проводить в двумерной постановке (в цилиндрической или декартовой системе координат соответственно). Есте-
ственно, в первом случае не учитывается асимметрия в изготовлении и инициировании заряда ВВ, во втором — протяженность и пространственная кривизна реальных КЗ. Кумулятивные заряды, формирующие ПЭ, принято также называть снарядоформирующими зарядами (СФЗ).
Типовые расчетные схемы рассматриваемых взрывных устройств показаны на рис. 1. Здесь 1 - ВВ, 2 - корпус КЗ или осколочного макета, 3 - точка инициирования (ТИ), 4 - КО, 5 - фронт детонационной волны (ДВ), 6 - ПД, 7 - линза (или экран), 8 - ре-перные точки (маркеры), используемые для идентификации течения в КО или осколочной оболочке.
а) б)
Рис. 1. Расчетные схемы взрывных устройств: (а) КЗ; (б) осколочный макет, где 1 - ВВ, 2 - корпус КЗ или осколочного макета, 3 - точка инициирования, 4 - КО, 5 - фронт детонационной волны, 6 - ПД, 7 - линза (экран), 8 - реперные точки (маркеры)
Так как для обеих схем постановка задачи примерно идентична, рассмотрим ее на примере КЗ. Для этого будем полагать, что в начальный момент времени ( = 0) в точке 3 (точка инициирования) осуществляется подрыв заряда ВВ с начальной плотностью рвв и теплотой взрывчатого превращения Q. От точки инициирования начинает распространяться фронт ДВ (кривая 5) со скоростью Овв с образованием ПД (зона 6). С течением времени ДВ начинает отражаться от поверхностей корпуса (позиция 2) и КО (позиция 4), на которые действуют давления порядка 20 ... 60 ГПа. Его величина зависит от свойств ВВ, угла подхода фронта ДВ к поверхности облицовки, материала и толщины облицовки. Под действием ПД кумулятивная облицовка начинает обжиматься с образованием КС или ПЭ. При этом для получения общих закономерностей или особенностей их формирования для конкретного КЗ, обусловленных формой облицовки, геометрией заряда, формой и месторасположением линзы и ТИ, физико-механическими свойствами используемого
состава ВВ или материалов облицовки и корпуса, целесообразно использовать модель сжимаемой идеальной упругопластической среды с баротропным уравнением состояния [1 - 4]. Вязкими свойствами материала облицовки в процессе формирования КС или ПЭ можно пренебречь. Использование баротропного уравнения состояния позволяет избежать интегрирования уравнения энергии в системе соотношений, описывающих поведение взаимодействующих сред. Причем если при формировании КС можно ограничиться газодинамической моделью, то для определения формы и кинематических характеристик ПЭ, особенно на поздних стадиях движения, является важным учет проявления упругопласти-ческих свойств материала. Для зоны ПД течение считается изэнтропическим, так как даже отражение фронта ДВ от абсолютно жесткой поверхности не приводит к заметному росту энтропии. Распространение детонации рассматривается вне общей системы уравнений и задает границу области, охваченную течением. С учетом сделанных допущений система уравнений, описывающая двумерное течение в переменных Эйлера, имеет вид:
д р д (р уг ) д (р ) а р V
д г
+
д г
+
р
d г
= р
д V„
д vr
—- + vг—г-д г д г
+ V,
ду^ д г
д г д а.
+
г д В
= 0
а
-+-агг + —(2 В + В )
а/., V агг аи!
г д г г
р
d V d г
р|
д
+ V..
д V д г
+ V.
д
д г
д а„ д В,
-+-- +а(в )
а/., V агг!
2 д г г
А ■ [р/р0 Уь -^ р/ро ^ 1 для КО, корпуса, линзы,
Р = р(р):
ат
^р/р о < 1
Врп + Сру
СВВ [(р/рвв )3 - 1
для КО, корпуса, линзы, для ПД, для ВВ;
а гг = Вагг - Р , а гг = Вая - Р , а гг = Вагг ^
ВВ„
вг
= 20
гдуг 1 dрл
—- +---
дг 3р dt
ВВ агг =
вг
О
ВВп7
_а
Вг +
= 20
г 1 dрл —- +---
дг 3р dt
дг дг
/ = 2 • (В2агг + В2т + В^ + ВаГГ • В агг )< | У2
(1) (2) (3)
(4)
(5)
(6)
(7)
Здесь р - плотность; р - давление; г - текущее время; уг, - компоненты вектора скорости по осям 0г и 0г выбранной системы координат (рис. 1); огг, о22 - нормальные и Огг - касательные напряжения; Вагг, Вогг, Вогг - компоненты девиатора напряжений;
г
*
-0(...) / Dt - производная в смысле Яуманна, учитывающая поправки к составляющим де-виатора напряжений, обусловленные поворотом фиксированного элемента среды как целого [2, 3, 5, 6]; О, У - текущие значения модуля сдвига и динамического предела текучести материала среды, принимаемых в областях разгрузки в виде функциональных зависимостей от их начальных значений О0, У0 и безразмерной плотности у = р /р0, причем при О0 = 0 и У0 = 0 исходная система уравнений естественным образом трансформируется в уравнения идеального газа, описывающие поведение ПД; а - коэффициент симметрии (для плоского случая а = 0, для осесимметричного случая а = 1).
В приведенной системе уравнений соотношения (1) - (3) представляют собой, соответственно, законы сохранения массы и импульса; (4) - баротропные уравнения состояние взаимодействующих сред (КО, корпуса, линзы, ВВ и ПД), конкретизируемые ниже. Соотношения (5) связывают компоненты тензора полных напряжений с шаровой и девиа-торной составляющими; (6) - закон Гука в дифференциальной форме; (7) - условие пластического течения Мизеса ог- < У (ог- - интенсивность напряжений). Если условие (7) нарушается, т.е. (/ > 2У /3) и материал деформируемой среды находится в состоянии пластического течения, то компоненты Вагг, Оаг2, Па22 умножаются на множитель У(2/3/)°'5 (процедура приведения напряжений к кругу текучести [6]).
Для описания разрушения металлических элементов конструкции зарядов под действием ПД использовалась комбинация критериев откольной прочности р = - а р (а р - от**
кольная прочность) и предельных пластических деформаций 8р < 8 р (в р - предельная пластическая деформация). Первый из них использовался на этапе нагружения элементов
конструкции ДВ, второй - при последующем пластическом деформировании. Причем при
*
р/р0 < 1 и р = - ар - величина давления ограничивалась в соответствии с условиями [7, 8]
р = р(рЖу) ; у = уо ф(у); о = оо ф(у),
(8)
где
1 для у>у1,
фМ=НУ- У 2 У(У1- У 2) для У 2 <У<У^ (9)
о для у<у 2;
71, 72 - постоянные материала, числовые значения которых, например, для стали по данным работ [9] составляют - у1 = 0.95, у2 = 0.90 для условий нагружения пластин и у1 = 0.975, у2 = 0.95 - для цилиндрических оболочек.
В качестве уравнений состояния материалов облицовки, линзы и корпуса заряда обычно используются ударные адиабаты в форме
р = р(р)=л [(р/р0 т -1,
где р0 - начальная плотность; Аь, пь - эмпирические константы [1], причем при пь Ф 1 это уравнение называется ударной адиабатой Тэта, а при пь = 1 оно вырождается в линейную баротропную зависимость, в которой Аь = К0, где К0 - модуль объемного сжатия.
Для продуктов детонации в качестве уравнения состояния рекомендуется использовать изоэнтропу в форме степенного двучлена
р = В рп + С ру .
Коэффициенты изоэнтропы В, С и п (у = 1.2.1.4) определяются по параметрам в точке Чепмена - Жуге [1 - 3]
(к-у)рс; 0_рШ-(У- 1) Ре; ес; / п - 11 С = ра - В Рс;
п = 1 +-ч 'ТС;-, В =
п - у
РС
рс; -(У-1) Р С; ес/ Р С;
где рС}, ре;, еС}- давление, плотность и энергия на фронте ДВ, к - показатель адиабаты (для конденсированных ВВ к~ 3):
р вв ^вв к +1 „ рс; ( 1 1
рс; = Г Г , Р с; 7 р вв, ~
к +1 к 2
+ б.
\Р вв Р,
Для расчета детонации используется искусственный прием, согласно которому наиболее полное выполнение условий на фронте ДВ достигается подбором изэнтропы ПД и уравнения состояния непрореагировавшего ВВ
р = свв [(р/рвв)3 - 1],
в котором коэффициент Свв определяется из условия
0.4 < р/р„ < 0.75.
Область разностной сетки, заполненная воздухом, в представленной постановке рассчитывается приближенно. Для этого в качестве изэнтропы условного «воздуха» используется изэнтропа ПД.
Граничные условия для рассматриваемых задач в рамках идеализированных расчетных схем (рис. 1), задаются на участках поверхностей взаимодействующих сред. Так, например, поверхности КО и корпуса, контактирующих с воздухом, и поверхность продуктов детонации, истекающих из корпуса, допустимо рассматривать как свободные от действия внешних поверхностных сил. То есть, пренебрегая силами атмосферного давления - Суп1 = 0 (рп/ = рпу = 0 для ПД), где п - вектор единичной нормали. На поверхностях КО, линзы и корпуса, контактирующих с ПД, накладываются ограничения на скорости движения индивидуальных точек в соответствии с условием непроницаемости
у(КО)п, = уПД)п, . улинза)п, = у(ВД)^ . у(корп)^ = „(ВД)^ ,
а также на напряженное состояние, реализующееся в этих точках в соответствии с третьим законом Ньютона
_ (КО) , = р (ПД). _ (линза) , = р (ПД). _ (корп) , = р (ПД) ^Ч Гт ' Ч Гт ' Ч Гт
При формулировке граничных условий на оси симметрии (ось 0г) необходимо учитывать, что при г = 0 частицы среды движутся только в осевом направлении (уг = 0), а осевые ускорения этих частиц должны быть ограничены. Из уравнений движения (2) - (3) следует, что это может быть реализовано только при отсутствии касательных напряжений на оси симметрии (оГ2 = 0). Давление, плотность и энергия во фронте ДВ задаются равными параметрам в точке Чепмена - Жуге.
Выделение контактных разрывов типа ПД - корпус или ПД - КО осуществлялось методом «концентраций» [2, 3], разработанного В.В. Кореньковым. Для этого в дополнение к основным характеристикам течения (р, р, уг, уг, ВоГГ, ВОГ2, В022) для неоднородной системы ПД (ВВ) - КО (корпус) определяют массовую концентрацию веществ следующим образом:
М, [1 для ВВ, ПД или воздуха,
щ = —1 = <
М [0 для корпуса. Здесь М1 - масса ВВ, ПД или воздуха в ячейке, а М - масса ячейки. Таким образом, для однородных ячеек, содержащих только первую компоненту щ = 1, для однородных ячеек со второй компонентой щ = 0 и, наконец, для смешанных ячеек 0 < щ < 1.
Локализация контактных границ осуществляется из анализа текущего распределения концентрации взаимодействующих веществ, которое определяется законом сохранения концентрации
д (р щ) д (р щ уг) д (р щ ) р щ уг
у + -г— + -У--+ а —-^ = 0. (10)
д ? д г д г г
Давление в смешанной ячейке вычисляется из условия аддитивности удельных
объемов и равенства давлений содержащихся в ней компонент
\Ур =(1 - №=1 + р щ=о; (11)
,Ро(р щ=0) = Р1(р щ =1 X где р, = 0, р» = 1 - текущая плотность смешанной ячейки и плотности содержащихся в ней
компонент; р0(рщ = 0) - ударная адиабата материала корпуса; р1(рщ = 1) - изэнтропа ПД или
уравнение состояния ВВ. По значениям поля концентраций определяется вид контактной
границы и на его основе рассчитывается поток из смешанной ячейки. На разрывах типа
ПД - КО или ПД - корпус в процессе решения обеспечиваются условия равенства нормальных составляющих скоростей и напряжений. Выделения границ ПД - воздух не производится.
Начальные условия конкретной задачи задаются распределением параметров р, р, w, vr и V в поле течения. Компоненты напряжений принимаются равными нулю.
Для получения более полной информации о процессах схлопывания КО или метания корпуса последние обычно маркируют реперными точками — маркерами или трассерами (позиция 8 на рис.1), в которых дополнительно вычисляют параметры текущего состояния среды Pm = {р, р, vr, V ...} по формуле
Pm =Е Р (к) Ак ]/(Дт Дг). (12)
к=1
Здесь к = 1 2, 3, 4; Ak- площади, определяемые текущим положением маркера m [2 - 5], а Pm(k) - параметры состояния среды в ячейках (к), окружающих маркер.
Система уравнений (1) - (11) интегрировалась модифицированным методом крупных частиц, подробно изложенным в следующем разделе (разд. 2), на программном комплексе КОЬОЦК [4], позволяющем проводить вычисления во всей области интегрирования на неподвижной (эйлеровой) сетке без предварительного анализа течения.
2. МЕТОД ЧИСЛЕННОГО РЕШЕНИЯ
а)
2, п+2 1, п+2 т+1, п+2
1, п+1 2, п+1 2, п+1 т+2, п+1
1,3+1
и Дг 1+1.1 т+2, /
А г
1, 2 2, 2 1, 2 т+1, т+2, 2
0 2, 1 1, 1 т+1, 1
б)
Рис. 2. Структура неподвижной (эйлеровой) разностной сетки: а) область интегрирования задачи (0 < г < Z, 0 < г < R); б) индексация ячеек разностной сетки («фиктивные» ячейки выделены серым цветом);
Г1, Г2, Г3, Г4 - границы расчетной области
Дискретизация задачи. Пусть в области интегрирования задачи (0 < z < Z, 0 < r < R) происходит плоское или осесимметричное неустановившееся движение неоднородной среды, описываемое уравнениями (1) - (10), при заданных начальных и граничных условиях. Область интегрирования разбивается на некоторое число прямоугольных ячеек со сторонами Ar = R/m и Az = Z/n, которые образуют неподвижную (эйлерову) разностную сетку (рис. 2, а). Значения целых чисел i = 2, 3, ... , m+1 и j = 2, 3, ... , n+1 обозначают центры ячеек разностной сетки (рис. 2, б), где m и n - число ячеек по осям r и z соответственно (для конкретного заряда это число обычно задают таким образом, чтобы по минимальной толщине оболочки оно было бы равно не менее 4 ... 5 целых ячеек, причем Ar/Az = 0.8 ... 1.2). В случае плоской симметрии отдельные ячейки представляют собой прямоугольники, в случае цилиндрической симметрии — тороиды. В точках (i, j) определяются средние значения параметров течения среды — компонент вектора массовой скорости, плотности, давления, концентрации, компонент напряжений.
Область интегрирования слева ограничена осью симметрии или жесткой стенкой (Г1, на рис. 2, а); снизу, справа, сверху — открытыми поверхностями (Г2, Г3 и Г4), через которые среда может вытекать или втекать. Чтобы не нарушать единообразия вычислений для граничных ячеек, вдоль всех упомянутых границ вводятся слои «фиктивных» ячеек (см. рис. 2, б), параметры которых определяют из смежных ячеек, как и в методе «частиц в ячейках» или в методе «крупных частиц» [2, 3]. Например, для оси симметрии или жесткой стенки (поверхность Г1) имеем:
)1, j =- )2 j; (~z )1, j = (~z )2, j;
(13)
\n+1 \n+1. \n+1 \n+1 . «+1 «+1 v '
(vr )n j=-(vr )2, j ; (vz )n J= (vz )n, j ; p1, j =p2, j ,
где параметры помеченные тильдой определяются в конце I-го этапа вычислений, а параметры с верхним индексом (n+1) — в конце III-го этапа [2, 3].
Для открытой поверхности Г2:
(~r )i, 1 = (~r )i, 2 , (~z )i, 1 = (~z )i, 2 ,
\n+1 Л n+1 \n+1 / , \n+1 «+1 «+1
(vr )n 1= (vr )i, 2 , (vz )n 1= (vz Л, 2 , Pi, 1 = P,, 2.
Для открытой поверхности Г3:
(vr )m+2, j = (~r ) m+1, j , (vz ) m+2, j (~z ) m +1, j (v ) n+1 = (v ) n+1 , (v ) n+1 = (v ) n+1 , Pn+1 = Pn+1 .
r m+2, j r m+1, j z m+2, j z m+1, j m+2, j m+1, j
(14)
(15)
Для открытой поверхности Г4:
(~г )г, п+2 = <Л )
Г /I, п+1 ■
(~г ),, п+2 = (~г )1
г/1, п+2 V г/1, п+1 э п+1 п+1 п+1
V) г, п+2 = V) г, п+1 V ) г, п+2 = V) г, п+1 г, п+2 г, п+1
Конечно-разностная схема. В методе «концентраций», который является модификацией метода «крупных частиц», в дополнение к расщеплению исходной системы уравнений по физическим процессам на три этапа (эйлеров, лагранжев и заключительный) проводится геометрическое расщепление уравнений по пространственным координатам. Реализация такого подхода, позволяет существенно упростить конечно-разностную схему физического расщепления двумерной задачи и избежать использования искусственной вязкости для размазывания скачков уплотнений.
В процессе отработки разностной схемы было рассмотрено несколько вариантов геометрического расщепления исходной системы уравнений (1) - (10) по пространственным координатам. Здесь излагается наиболее простой из них, в котором вместо соотношений (1) - (3) и (10) используются две подсистемы [2, 3]:
1 д Р + д (Р V ) + а Р V
2 д г
д г
_ 1 д (р w) д (р w уг ) р w Уг
0,--——- + -— + а —-г
1
2 д г
+ У„
д V,,
г 1
2 д г
д г р
д (- р) + + а(2Д + О )
а/., V агг агг!
г д г
д г а г
0,
1 д V.
2 д г
+ V..
д ^ д г
0
1 д р д (р Уг)
2 д г
+
д г
= 0
1 д (р w) д (р w Уг)
2 д г
д V.
1 д уг
—- + Уг —1
2 д г г д г
+
= 0,
1 д V д V 1
--- + vг —- = -
2 д г д г р
д (Оагг - р) + д Д
д г
а
+ —
г
= 0
(О агг )
(17)
(18)
д г д г
Последовательное интегрирование уравнений (17) и (18) с шагом т = Дгп/2 по каждому из направлений г и г в предположении неизменности компонент девиатора напряжений на текущем временном слое (слой п) позволяет получить значения состояния среды рп1, рп1,
vгn+1, V/1, wn+1 на новом временном слое гп1= гп + Дгп. При этом в обоих случаях устойчивость вычислений обеспечивается выбором шага интегрирования по времени
,.п+1
п+1_п
А гп =
шт
г, у
К
V
13
1, У
maX(| ^ ¡, I Vг |) + С"
1, У
(19)
г
где ¡, ] - индексы текущей ячейки; V¡,, С¿/, (уД/, (у)/ - объем, скорость звука материала и компоненты массовой скорости текущей в ячейке; Кг = 0.1 ... 0.35 - число Куранта.
На следующем этапе по вновь полученным характеристикам течения среды р ,
(уг )Г/!, (Уг)"+ рассчитываются значения производных (5 р/д 0П+и2, (д уг/д г)П+и2, (д уг/д 7)Г+1/2, (д уг/д г)П+12, (д уг/д z)П+1/2 и на основании соотношений (6), (7) определяются обновленные компоненты девиатора напряжений (Багг )г"++1, (Баи )г"++1, ()П+.
Далее излагается разностная схема для рассмотренного варианта расщепления исходной системы уравнений. Для этого будем полагать, что в текущий момент времени tn известны все величины потока двухкомпонентной среды в центрах ячеек эйлеровой сетки
РП,, к, у )",, у )П;, рП, , (Багг) I,, )П;, (Бт )П;. Получим соотношения для вычисления перечисленных параметров внутри и на границе расчетной области в момент времени ^+1= Г + Д^, где Д^ - шаг по времени, определяемый соотношением (19).
Рис. 3. Индексация границ ячейки на неподвижной разностной сетке
Этап I (эйлеров). На этом этапе вычислений среда предполагается «замороженной», т.е. отсутствует массообмен субстанцией между смежными ячейками разностной сетки. Изменение параметров течения среды происходит за счет напряжений, действующих в элементарном фиксированном объеме. В этом случае уравнения (17) и (18) принимают вид:
р = const, wp = const,
1 д (d_„„ - p) д d_,„ а р д r д z r
1 д vr д vr
--- + vr —-
2 д t д r
ZJD^.+а(2 От + DJ
д r д z
(20)
1 ду
2 д t
+ v„
д Vz
д z
д (Dc
д z
Р) д Dc
± / + с
д r
О„)
г
Заменим уравнения (20) разностными. Для этого определим конечно-разностные аналоги производных типа (д //д г)». и (д //д г)». следующим образом:
Л Л» -|
дГ] = л г /»)=— 5 / »- / ») - ^ /п - / »)]=
-/з^з - / - 5 - 5,)]= [5з(¿л,. - /;;)-- /;,»)]
= 1 [ /
/з }
д / ] »
д 2 у ;,}
(21)
2К
;,}
= Лг (/» )= ^ [54 /» - /» )- 52 (/2 » - /,» )] = ^ [/ »54 - /2 »52 ]
(22)
Здесь /} - значение некоторой сеточной функции в центре ячейки (;,}) на »-ом временном слое (рис. з); /", /г, /", /4" - значения произвольной сеточной функции на левой, нижней, правой и верхней границах ячейки (;,}) соответственно, т.е.
1111
г п _ ( 4 " 4 " ) 4 " _ I 4 " л. 4 п\ 4 " _ 14" I I 4 " _ 14" л. 4 "г
/1 -1,. /;,}', /2 _ 2 ,}- ;,} /з ~2^/;+1,1/;,./, /4 _ 2 „}+ ,}
У;у и 51, 52, 5з, 54 - соответственно объем и площади левой (1), нижней (2), правой (3) и верхней (4) граней (;,})-ой ячейки, определяемые формулами, приведенными в табл. 1.
С учетом соотношений (21) - (22) и величины шага интегрирования по времени т = Ы»/2 определяются промежуточные значения массовых скоростей течения среды (~г)»}. и
)"у в ячейке С
(~г );,} = (^г )» , +
А г
гл,} \ г /г, } р »
',}
Л [(а )]+Л [(В )".]+ — (2В + В )"
г гг ; , . г ; , . гг ; ,
;, .
(~);,} = (V ),\ +
А г
г;;,} \ г;;,} р»
;,}
Л г [(а гг ) » . ]+Л[( Бп ) ». ]+—(£>„ )
гг>;,}
;,}
Табл. 1
Формулы, определяющие объем и площадь граней фиксированной ячейки в плоском и
осесимметричном случаях
Параметр Вид симметрии
а=0 (плоская) а=1 (осевая)
II /Я — Дг [(/ -1) -12] АгАг
Я2 = Я/, 7-1/2 Дг (/ -1) АгАг
= Я/+1/2,7 Дг [(/ -1) +1/2] АгАг
Я 4 = ,, 7+1/2 Дг (/ -1) АгАг
К' Дг Дг (/ -1) Аг2 Аг
Этап II. (лагранжев). На этом этапе вычисляются эффекты переноса среды, учитывающие массообмен между отдельными ячейками эйлеровой координатной сетки. При этом предполагается, что перенос массы сплошной среды через границы ячеек за промежуток времени осуществляется с учетом промежуточных значений компонент Уг и вектора массовой скорости среды, рассчитанных на первом этапе.
: 0)
©
г-12 1+1/2
1-1 j
гф
1>3
■¡+и
©
Д 4
(~'гЬ+1/2,3 0 а)
¿0)
©
1-12 1+1/2
¡-и
гф
1>3
А
¡+и
Д
1+1/2,3
©
(!>),., , , , < О
б)
Рис. 4. К определению потока массы среды через границу между смежными ячейками (/,'') и (/'+1,7'): а) (~г )+1/2 . > 0 ; б) (уг )+1/2 . < 0; Д - донорская ячейка; А - акцепторная ячейка
В начале определим поток массы среды с одинаковой концентрацией (н =0 или н = 1) через границу между ячейками. Например, для правой границы произвольной ячейки (рис. 4) он определяется как
Ат
/'+1/2, '
Р/,' (У )/'+1/2,' ' Я/ +1/2,' А" пРН (У )/+1/2,' > 0;
г //+1/2, у
Р/+1,' (У )г+1/2,' * Я/+1/2,' пРН (У )г+1/2,' < 0
(23)
Здесь р. р;+1}- - плотности сред, содержащейся в смежных ячейках; 5;+1/2. -площадь правой границы ячейки (;, .), через которую производится массообмен между ячейками (;, .) и (; +1, }); (~г)м/2} = 12• [(~г);+1} + (~г);}]. Согласно выражению (2з) поток массы среды с концентрацией н = 1 (Дда1) через ту же границу ячейки рассчитывается по формуле (24)
(Ат1),
;+1/2,}
(Ат);+1/2,} • } при (~г );+1/2,} > 0;
~ (24)
(Ат);+1 /2,} • Н+1,} пРи <Л);+1/2,} <
Обратимся теперь случай, когда одна из смежных ячеек является смешанной (0 < н < 1), а другая — однородной (н = 1 или н = 0). Для примера будем рассматривать поток массы среды через правую границу смешанной ячейки (;,}) (см. рис. 4). При этом назовем ячейку, из которой происходит истечение среды, донорской (Д), а другую, в которую направлен поток, — акцепторной (А).
Как указывалось выше, в методе «концентраций» принимается, что из смешанной ячейки в однородную первоначально вытекает то вещество, которое находится в однородной ячейке, а обмен массами между смешанными ячейками происходит пропорционально объемным концентрациям веществ в акцепторной ячейке. С учетом сделанных допущений рассматриваемый поток массы среды определяется следующим образом:
А т+1/2,} = м 1 + м o,
м 1 =
мп
(25)
где
а ;+1, } • Л (Ат1) ;, } ПрИ (~г ); +1/2, } > 0,
а;,} • й(Ат1);+1}. при (~г);+1/2, < 0;
(1 - а ;+1,} ) • й (Ат0);,} при (~г );+1/2,} > 0, (1 - а;,} ) • й (Ат0) ;+1,} при (~г );+1/2,} <
Здесь м1, м0 - первая и вторая составляющие потока массы среды через правую границу смежной ячейки (;,}), соответствующие концентрации н = 1 или н = 0. При этом последовательность расчета величин м0 и м1 определяется значением концентрации вещества, содержащегося в однородной ячейке. Здесь а;. = (н;,}р;,})/(рн=1);,}, (1 - а. - объемные концентрации первой и второй компонент среды в (;,})-ой ячейке; й(Дт0)Ь, й(Дт1)Ь, Ь = [(;, }), (;+1,})] - потоки массы двухкомпонентной среды, определяемые текущими значениями плотностей (рн=0)Ь и (рн=1)Ь, а именно:
й (А Щ)Ь = (5 • ~);+12, } • (Рн=0)Ь А г П, й (А т)Ь = (5 • ~);+/2,} ■ (Рн=1)Ь А г ». (26)
Отметим, что соотношения (26) справедливы только в частном случае, когда в течение временного шага Дг" компоненты среды, находящиеся в смешанной ячейке, не
ЬИр:/ЛесЬпотац^и.гиМос/зз4177.Ь1:т1
1з
успевают полностью перетечь в акцепторную ячейку. С целью устранения указанного недостатка в общем случае необходимо предусмотреть возможность корректировки потоков перетекающей массы.
Рассмотрим вариант такой корректировки для потока М1 между смешанными ячейками (/, ]) и (/+1, ]) при условии (уг ) .+у2- > 0. Для этого обозначим массу вещества с
признаком н = 1, находящуюся в донорской ячейке (/,/), через (Ат^,-. Если (Ада1)г-7- стала меньше перетекающей массы М1, определяемой соотношением (25), значит, все вещество с признаком н = 1 вытекло из ячейки через правую границу. При этом время истечения составляет А^ = АС1 (Ат1/М1)гу. За промежуток времени А^2 = АС1 - А^ через площадь а/+1- Si+1/2,j в ячейку (/ +1, -) будет втекать второй компонент среды с признаком н = 0. С учетом сказанного составляющая М1 потока массы среды пересчитывается следующим образом:
М1 = (А тД. у + й (А то),у -ам.
=(А т1)/,у + й (А то)/,у 1 -
С л Л 1 А т1
М1 (А т1)
А, з з
а.+и-й (А т1) .,у
Сделав аналогичные выкладки для М1 при (уг ).+^2 - < 0 и для М0 при (уг )- > 0 или (уг ).+12 - < 0 , окончательно можно получить следующие соотношения:
М,
а н й (Ат1) Ь при а н й (Ат1) Ь < (Ат1) С (Ат1) Ь + й (Ат0) Ь -а, 1 (Ат1) 11 ^
Мп
а м й (Ат1) Ь при а н й (Ат1) Ь > (Ат1) Ь;
(1 - а^ ) й (Ат0 )ь ПрИ (1 - а^ ) й (Ат0 )Ь < (Ат0 )Ь
(Ат1) Ь С
С(Ат1)ь + й (Ат0)ь (1 -аN )
1
(1 -а N) й (Ат^ 1
пРи (1 - аN ) й (Ат0 ) 1 > (Ат0 ) 1
(27)
(28)
Значения нижних индексов Ь, N и константы £ в соотношениях (27) и (28) в зависимости от направления движения среды )+у2 - ^ 0 или )+у2 - < 0 приведены в
табл. 2.
Табл. 2.
Значения нижних индексов и константы в уравнениях (27) и (28)
(vr )i+12, j Z L N
> 0 1 i+1, j
< 0 -1 i+1, j
Давление в смешанной ячейке, содержащей контактный разрыв, определяется из решения системы нелинейных уравнений (11), представляющих собой условие аддитивности удельных объемов и равенство давлений находящихся в ней отдельных компонентов среды.
III Этап (заключительный). На этом этапе в пределах объемов ячеек эйлеровой координатной сетки определяются параметры состояния среды на временном слое n+1 с учетом переноса ее компонентов через границы ячеек:
p n+1 n , Am> -12, у -Amt +12, у . »+1 = Vt ■ X ] PL- + (Ami) t-1/2,j - (Ami) t +12, j . Vi, j Pu j + T/ p n+1 ' Wt,j v p n '
V. .р™ V .pl
(vk) I j = {? 1-12, j (~k) i-12, j Ami-12, ji+12, j (~k) i+12, j Ami+1/2, j +
+ (Vk) l j [pl V, j + (1 -Z i-12, j) Am.-12, j - (1 -Z i+12, j) Am.+i2,j ]}/p l+V, y.
Здесь k = r, z; Zl = 1, если поток массы среды через границу L = [(i - 1/2, j), (i + 1/2, j)], направлен вовнутрь ячейки (i, j), и ZL = 0 в противном случае.
Характерной особенностью рассмотренного алгоритма является тот факт, что деформируемая среда рассчитывается с учетом изменения шаровой компоненты напряжений на предыдущем полушаге. Последнее способствует увеличению аппроксимационной вязкости разностной схемы и подавляет счетные осцилляции, характерные для алгоритмов «частиц», без дополнительного введения искусственной вязкости.
Компоненты девиатора напряжений. Отличительной особенностью уравнений (6), определяющих изменение компонент девиатора напряжений в неподвижной системе координат по отношению к лагранжевой форме записи аналогичных уравнений, является наличие так называемых конвективных составляющих полных производных величин Darr, Darz, Dazz. Для получения их конечно-разностного аналога используется подход, предложенный Хагеман, Уолшем, Мейдером [5], в котором компоненты девиатора напряжений на временном слое 1+1 определяются по величинам (Darr )l j, (Darz )1 j, (Dazz )1 j, и
промежуточным значениям приращений (АОагг, (АОаг2(АВа2 )г"++1 , взятых с соответствующими весами. При этом последовательность вычислений следующая.
Во-первых, с учетом граничных условий (13) - (16) в каждой ячейке определяются производные вида
'д к ^^)£ у - V)Щ у.
V д г у
д V, л
д г
1,з
п+1
2 Аг
к )п+;+1 - к )п+-1
(V, )П+1 з - V )П+
п+1 2/1-1, ]
'1,з
2 Аг
Кд 2 У д 2
1,з
п+1
2 А2
(V,) П+1+1 - (V, )П+М
+1
2'и ]-
'1,з
2 А2
и далее
д кг
Уд Г У1,з Vд г У,,з
'д Ип+1 '
Vд г
+
д кг
Уд г У
Гд ^п+1 '
Vд г
+
дк1 Vд г У
1, з
1, з
д кг Vд 2 У и з 'д^ V д 2 У и з
(д ип+1 '
V д 2
+
V
д кг
Уд 2 У
Гд ^п+1 '
V д 2 У
+
п
д V V д 2 У
ь з
1, з
Во-вторых, вычисляются приращения девиатора напряжений:
X п+1/2
где
(АОт Г = Агп
(АОа„ Г = А Iп
(АОаа Г = Аг п О
Г 1 а \п+112 1 д р
3р д г
1,з
l, з
2О (д 1 д р
г +
V д г 3р д г
2О (д V2 1 д р
2 +
V д 2 3р д г
/д V дV > п+1/ 2
V д 2 д г уУ г, з
1 п+1 п
р,з - р1 , з
+ р п з А г п
+
ж (в„.
з
х п+1/2
- V У' ,
)
- ж(в - В )"
V агг а22 /г,
2
А г
( ~ п+1 _ п р -р
~п+1 . _ г Vр +р
'1, з
С П П Л п+1/2
' д V д V,
'1, з
2
Ж :
V д 2 д г
В-третьих, определяются значения составляющих девиатора напряжений на временном п+1 слое с учетом конечно-разностного аналога конвективных производных величин (Вагг )гп+\ уп+1, (Ва22 уп+1, взятых с соответствующими весами
§ г, з = (vr)А гп и 5 21, , = V)А гп.
^1, з
1, з
2->1, з
п
Рассмотрим конечно-разностные соотношения для расчета одной из компонент. Так как для других компонент соотношения аналогичны, то далее нижние индексы у составляющих девиатора напряжений не указаны.
1-1, 7
4 1
3 2
1-1И
Злу
2 3
1 4
а)
б-/./
./
■"гг.]
=С=!
в)
А./+/
ьз
б)
У 1+1,}
1 4
2 3
-1 1+1, И
г)
Рис. 5. К определению разностной записи конвективной производной: а) Ьгу < 0, Ьгу < 0; б) Ьг.у > 0, Ьгу > 0; в) Ьг.у < 0, Ьгу > 0; г) Ьу > 0, < 0
Если величины Ьг, и Ьгу, отсчитываемые от центра ячейки (., у), отрицательны (рис. 5, а), то имеем
4
пй+1. = О . ,
, у ^о., у
^ [(АО о) Г1 • АГ']к
Аг Аг
д=1
где
АГ+1 = = АГ, ++1 = (А г -1 5 г.,у |) (А г -1 5 г.,. |), (АОо )Г+1 = (АОо) Г+;
АГ+1 = А А = (А г -1 5 г, у|) | 5 г., - |, (АОо) 2+1 = = (АОо )Г+-1;
АП +1 = -1, у- -1 = | 5 Гу || 5 г. у|, (АОо )Г+1 = (АОо ^-1;
АГ +1 = Аг-+:- = | 5 г, у |(А г -1 5 г., у), (АО о) 4+1 = (АОо) Г+1, у.
При г = 2 (см. рис. 2) —
(AD, ) П+ ] = (AD, ) ] (AD, )n-+l^ ]-l = (AD, ) ];
при ] = 2 —
(AD, ) n+-l = (AD, ) ] (AD, )n+1]-l = (AD, ) ^ ;
при г = 2 и ] = 2 —
(AD, ) n+1 ] = (AD , )П+ ]-l = (AD, ) n++-l = (AD, ) ] Если ôr.] > 0 и ôz, > 0 (рис. 5, б), то
a: +1 = An,+ = (A r - S r,] ) (A z - S zt,] ), (AD,Г = (AD,)n+ ;
A2+1 = An,++l = (A r - S r.,] ) S zt,], (AD, )2+1 = (AD, )] ;
АП +1 = АП++1]+l = S r-, ] S z., ], (AD, )n+1 = (AD, )П+ ]+l ;
a:+1 = Aù = S r-,] (A z - S zг,] ), (AD, )4+1 = (AD, )П+].
При i = m+1 (см. рис. 2) —
при ] = n +1—
(AD, ) n+ll ] = (AD, )] (AD, ) n++l1 ]+1 = (AD, ]
(AD, ) П++1 = (AD, )] (AD, ) n+l ]+l = (AD, ) П++^ ] ;
при i = m+1 и ] = n +1—
(AD, ) ] = (AD, Х-, = (AD, ) n+ll ] = (AD, ) ]
Если ôrij < 0 и ôztj > 0 (рис.5, в), то
АП+1 = АП Г = (A r -1 S r-,] I) (A z - S z., ] ), (AD, )П+1 = (AD, Xl ; A2+1 = АП^ ++1 = (A r -1 S r-, ] I) S z., ], (AD, ) 2+1 = (AD, ) ] ;
АП +1 = АП-+1]+1 = IS r-, IS zt,], (AD,)П+1 = (AD,)n_+l1]+l;
АП+1 = Aû = I S r-,] I (A z - S z,] ), (AD,)П+1 = (AD,)П+].
При i = 2 (см. рис. 2) —
(AD, ) £1 ] = (AD, ) ] (AD, ) П-+1^ ]+l = (AD, )
при ] = n +1—
(AD, )] = (AD, ) ] (AD, )П-+1^ ]+l = (AD , ) ^ ;
при i = 2 и ] = n +1 —
(AD, ) П-11] = (AD, ) n++l = (AD, ) П-11]-1 = (AD, )]
Если 5гц > 0 и < 0 (рис 5, г), то
А" +1 = АП, 7+1: = (А Г -5 Г- 7 ) (А 2- | 5 2г,7 |), (АД.)"+1 = (АД.) "+1;
А2+1 = А", +-1 = (А Г -5 ги7) | 5 2г,7 |, (АД )"2+1 - (АД.) "+-1;
А+1- - А +1 - Аг +1,7 - 1 = 5 Г,7 | 5 ^ г,7 1 (АД. )"+1 = (АД.) "+117-1;
А4 +1 = Л" +1 Аг+1,7 = 5 Г,7 (А 2- | 5 2г,7 ), (АД.) 4+1 = =(ад. )"+ 7'
При I = т+1 (см. рис. 2) —
(АД.) Й у = (АД, )"+\ (АД.) "+; ^ = (АД.) ^
пРи7 = 2 —
(АД.) "+- = (АД. )"+, (АД.) "+; 7-1 = (АД.) ^;
при I = т+1 и при 7 = 2 —
(АД.) £-1 = (АД.) ^-1 = (АД.) Х = (АД.) "+1;
В-четвертых, реализация условия пластического течения обеспечивается выполнением процедуры приведения напряжений на круг текучести [2, 3]. Для этого в каждой ячейке рассчитываются значения
г+1 = 2(д2 + Д2 + Д2 + Д - Д )
•/г, ] \ .ГГ 2Г2 .22 .ГГ .22 /I
П+1
г, 7
В случае//+1 > (2/3)72 производится корректировка компонент девиатора напряжений по формулам
(д' )"+' = ^ .(д г1, (д' у+1 = ^ .(д )й+\ (д' у+1 = ^ -(д г1,
V .ГГ /г, 7 V .гг/г, 7 ' V .22/г, 7 V .22/г, 7 ' V .Г2/г, 7 V .Г2/г, 7 '
^ 7
^ 7
/г, 7
где скорректированные значения компонент девиатора напряжений определяются следующим выражением
^ = 7
1
2
3/,
"+1
г, 7
(Д ' )"+\ (Д ' )"+\ (Д ' )
V .гг /г,7 ' V .22 /г,7 ' V .Г2 /
\"+1 7 '
Реперные точки. Для получения более полной информации о процессах схлопыва-ния облицовки или метания корпуса последние обычно маркируют реперными точками (маркерами или трассерами), в которых по формуле (12) дополнительно вычисляют параметры текущего состояния среды, см. рис. 6.
1 7
А, 1..
/
Л4 А«
Маркер !
4
Рис. 6. К определению параметров течения среды в маркере
3. ВЕРИФИКАЦИЯ РАЗРАБОТАННОЙ методики расчета
До настоящего времени для большинства сложных нестационарных задач физики взрыва не доказаны математические теоремы существования и единственности получаемого решения, поэтому тестирование является важнейшим атрибутом проверки степени соответствия построенной математической модели и ее численного аналога реально протекающему физическому процессу. Для решения вопроса о достоверности получаемых результатов использовались различные приемы проверки качества расчета, изложенные в работах [2, 3]. А именно: проверка выполнения интегральных законов сохранения массы, импульса и энергии; оценка результатов решения одной и той же задачи по одному и тому же алгоритму, но с разным шагом пространственной разностной сетки; сравнение результатов расчета с точными решениями известных задач; сравнение результатов, полученных разными методами; использование экспериментальных данных и другие.
а)
б)
Рис. 7. Конструктивные схемы СФЗ - а) №1, б) №2: 1 - ВВ, 2 - корпус, 3 - ТИ, 4 - КО
а)
б)
Рис. 8. Осколочные цилиндрические макеты - а) №1, б) №2: 1 - ВВ; 2 - цилиндрическая оболочка; 3 - ТИ; 4, 5 -верхнее и нижнее донья соответственно
Сравнение с экспериментальными данными и альтернативными методами расчета. Особое внимание уделено сопоставлению результатов расчетов с данными рент-геноимпульсной съемки. Для этого использовались лабораторные СФЗ с внутренним диаметром корпуса 50 мм (макет №1) и 64.1 мм (макет №2), показанные на рис. 7, и осколочные цилиндрические макеты с диаметром шашки 20 мм, показанные на рис. 8. Кумулятивные заряды снаряжались литьевыми составами ТГ-50 (рвв = 1.65 г/см3, Пвв = 7.7 км/с, 0 = 4.61 МДж/ кг) и ТГ-40 (рвв = 1.68 г/см3, Ввв = 7.85 км/с, 0 = 4.61 МДж/ кг), осколочные макеты — прессованными шашками из флегматизированного гегсогена (рвв = 1.65 г/см , Пт = 8.1 км/с, 0 = 4.82 МДж/ кг) или флегматизированного октогена (рвв =
1.75 г/см3, Ввв =
8.6 км/с, 0 = 5.23 МДж/ кг). Кумулятивные облицовки СФЗ изготавливались из пластичной стали 11ЮА, физико-механические параметры которой принимались следующими:
плотность - р0 = 7.85 г/см ; модуль сдвига - G0 = 80 ГПа; динамический предел текуче*
сти - У0 = 0.5 ... 0.8ГПа, откольная прочность - а р = 1.65ГПа. Для корпусов кумулятив-
*
ных макетов применялась сталь 45Х (00 = 77.81 ГПа; У0 = 1.05 ... 1.1 ГПа, а р ~ 2.3ГПа),
а для цилиндрических оболочек осколочных макетов - сталь 20 (00 = 76.79 ГПа; У0 =
*
0.43 ГПа, а р = 1.65 ГПа). В конструктивном отношении осколочные макеты отличались между собой лишь формой верхнего дна.
а)
б)
в)
г)
20 мкс
25 мкс
33 мкс
Рис. 9. Сравнение экспериментальной и расчетной формы поражающего элемента, образованного взрывом СФЗ №1, на начальном этапе деформирования (20 мкс, 25 мкс и 33 мм): а) данные рентгеноимпульсной съемки, б) расчет по программе КОЬБИК, в, г) расчет по программе ЛИТОБУК; г) КО без отверстия в купольной части
Результаты отдельных выполненных тестов представлены на рис. 9 - 12. На рис. 9 показано сравнение расчетной формы ПЭ (б, в) с данными рентгеноимпульсной съемки (а)
на начальном этапе деформирования поражающего элемента (20 мкс, 25 мкс и 33 мкс), сформированного СФЗ №1, конструктивная схема которого схематично показана на рис. 7, а. При этом рис. 9, б соответствует расчету по программе КОЬБЦК [2 - 4], а рис. 9, в - АиТОБУК. Для сравнения на рис. 9, г показан расчет по программе АиТОБУК аналогичного заряда с облицовкой без отверстия в купольной части.
Отличительной особенностью СФЗ №1 является то обстоятельство, что используемая в нем кумулятивная облицовка имеет достаточно сложную и не традиционную форму. Ее наружная (обращенная к ВВ) поверхность представляет собой сферический сегмент с радиусом образующей 36 мм, а внутренняя (обращенная к преграде) поверхность очерчивается образующей в виде параболической кривой г = (39.5г) +5 относительно локальной системы координат, показанной на рис. 7, а. В купольной части КО имеется отверстие диаметром 4 мм. Как было показано в [10, 11] при взрыве такого заряда образуется удлиненный ПЭ, который разрушается при деформировании в процессе движения к преграде. Поэтому СФЗ №1, получивший впоследствии название заряд Улякова П.И. [11], удобно использовать для моделирования процесса разрушения УПЭ из стали 11ЮА. Как видно из представленных иллюстраций, полученные расчетным путем формы и кинематические характеристики ПЭ на начальной стадии взрыва (до ~35 мкс) при использовании всех представленных программ расчета имеют хорошее соответствие реально протекающему процессу. Неясностей физического толка для математического описания начальной стадии взрыва СФЗ не возникает.
а)
б)
30 мкс
62 мкс
93 мкс
33 мкс
64 мкс
84 мкс
Рис. 10. Сравнение расчетной формы ПЭ с данными рентгеноимпульсной съемки: а) заряд
№2 (с корпусом); б) заряд №2 (без корпуса)
На рис. 10 представлены результаты расчетов СФЗ №2 в сравнении с экспериментальными данными. Как уже отмечалось выше, в нем использовалась сегментная облицовка постоянной толщины 2.2 мм. Ее внутренняя и наружная поверхности представляли собой сферические сегменты с радиусами кривизны 54 мм и 56,2 мм соответственно.
Рис. 10 иллюстрирует степень влияния корпуса (поз. 2 на рис. 7) на геометрические характеристики формируемого взрывом ПЭ. Для этого показаны характерные стадии деформирования поражающих элементов, образованных СФЗ №2 с корпусом (рис. 10, а) и без корпуса (рис. 10, б). При этом сравниваются формы ПЭ, полученные расчетным путем при Y0 = 0.8 ГПа, с данными рентгеноимпульсной съемки в моменты времени 30 мкс, 62 мкс и 93 мкс для заряда с корпусом и в моменты времени 33 мкс, 64 мкс и 84 мкс для заряда без корпуса.
На рис. 11 представлены расчетные данных функционирования осколочного макета №1 с толщиной боковой стенки 1 мм, наполненного флегматизированным октогеном и полученные с помощью программ KOLDUN и Ls-DYNA в сравнении с данными импульсной рентгенографии. На рис. 12 представлено сравнение данных расчета функционирования осколочного макета №2 с толщиной боковой стенки 1 мм, наполненного флегмати-зированным гегсогеном по программе KOLDUN и программе AUTODYN, позволяющим рассчитывать процесс детонации ВВ и пространственный разлет оболочки разными методами, в частности, безсеточным методом сглаженных частиц SPH и эйлеровым методом FLIC.
Рис.11. Сравнение расчетных и экспериментальных данных функционирования осколочного макета № 1, наполненного окфолом, с толщиной стенки стальной оболочки 1 мм на момент времени 13 мкс а) расчет по программе КОЬБИК; б) расчет по программе ЬБ-
БУКЛ; в) эксперимент
В целом же по результатам выполненного тестирования можно сделать вывод о приемлемости результатов расчета как по формированию ПЭ, так и по метанию осесим-метричных оболочек произвольной конфигурации, получаемых с использованием разработанной численной методики. Не исключено также использование разработанной программы для тестирования других программ как отечественных, так и зарубежных. Кроме того, время расчета идентичного варианта по программе КОЬОЦК, реализующей изложенный алгоритм, примерно в два три раза меньше, чем у зарубежных аналогов.
Возможности расчетной методики. Программный комплекс «КОЛДУН» или «КОЬБЦК» [2 - 4, 12], отдельные модификации которого назывались также «КОЛУН» [13] и ГЕФЕСТ [1, 13 - 17], предназначен для решения задач, описывающих высокоскоростные двухфазные нестационарные двумерные течения упруговязкопластических, жидких и газообразных сред механики сплошной среды. Ранее он широко использовался для моделирования процессов взрыва в воздухе [18] и в воде [19 - 21]; решения задач кумуляции [1, 4, 22 - 24], в том числе связанных с функционированием ударного ядра [10, 11, 25]; моделирования взрывного метания осесимметричных оболочек [1, 12, 14, 15, 17] и пластин [13, 16, 17, 26]; высокоскоростного проникания [1 - 3, 24, 25].
а) б) в) г) д) е)
10 мкс 15 мкс
Рис 12. Сравнение данных расчета, полученных по программам KOLDUN и AUTODYN
(SPH и Euler), функционирования осколочного макета №2 с толщиной стенки 1 мм, наполненного А-91 на моменты времени 10 мкс и 15 мкс: а, г) расчет по программе KOLDUN; б, д) расчет по программе AUTODYN (SPH); в, е) расчет по программе по AUTODYN (Euler)
Ниже приводится краткий обзор наиболее интересных результатов, полученных на программном комплексе KOLDUN как самим автором, так и независимо от него.
В работах [1 - 3, 26] представлены результаты расчетных исследований критических условий образования кумулятивных струй, а также влияния типа симметрии (плоская, осевая) облицовки и сжимаемости ее материала на скорость КС для установившейся фазы истечения. Они проводилось в рамках идеализированных расчетных схем косого соударения металлических пластин и конических облицовок с жесткой стенкой или осью 77-30569/334177, №02 февраль 2012 г. http://technomag.edu.ru 24
симметрии в диапазоне скоростей соударения 0.5 ... 3 км/с. Показано, что результаты расчета скорости КС из облицовок с углами раствора больше 40° (соответствующих реальным на практике конструкциям), полученные с учетом сжимаемости материала и осесим-метричности процесса, не имеют существенных отличий от результатов гидродинамической теории кумуляции [1], созданной на основе теории соударения плоских струй несжимаемой жидкости. Существенные отличия имеет место для режимов схлопывания, близких к критическим условиям формирования КС.
В работах [1, 22] численно исследовались особенности формирования кумулятивного «ножа» удлиненного заряда с клиновидной выемкой, на основании которых предложена инженерная методика расчета функционирования данного вида КЗ от момента инициирования до окончания пробития.
В работах [1, 23] на основании анализа результатов численных расчетов формирования КС с полусферическими и сегментными облицовками высокого прогиба (отношением высоты внутреннего контура к внутреннему диаметру облицовки более 0.2) построена инженерная методика, позволяющая с приемлемой для практики точностью проводить сравнительный анализ эффективности действия таких зарядов. При этом показано, что образование струи в этом случае происходит через выворачивание кумулятивной облицовки, а распределение скорости вдоль струи близко к линейному распределению.
В работе [25] оценивались возможности использования КЗ, формирующих компактные и удлиненные ПЭ, для поражения подводных преград, имитирующих однокор-пусные и двухкорпусные подводные лодки. Получены результаты по возможному использованию высокоскоростных ПЭ для пробития подводных преград в условиях ограниченных габаритов воздушной полости перед КЗ, отсутствия или наличия воды при движении ПЭ между легким и прочным корпусами объектов, геометрических и кинематических параметров самих ПЭ. Показана низкая пробивная способность стальных компактных ПЭ, движущимися со скоростями до 2.5 км/с, по водной преграде, которая не превышает двух калибров заряда, формирующего компактный ПЭ.
Решения задач высокоскоростного проникания представлено в работах [1 - 3, 24, 25]. В частности, в [1 - 3, 24] представлены результаты математического моделирования процесса проникания элемента КС в виде стержня в воду. Они позволили с достаточной для практических целей точностью описать особенности движения образующейся в воде каверны в зависимости от геометрических (диаметр, длина), энергетических (скорость) и физико-механических (инерционных) характеристик взаимодействующих тел как до полного «срабатывания» внедряющегося элемента (гидродинамическая стадия проникания), так и после момента «срабатывания» (инерционная стадия проникания).
Подобного рода информация особенно необходима для построения адекватных инженерных моделей, описывающих проникание КС в жидкие преграды (см., например, [27]). Отличительной особенностью такого функционирования является большая глубина проникания КС, достигающая 20 ... 40 диаметров заряда (для сравнения укажем, что глубина пробития КЗ стальных преград, как правило, не превышает 6 ... 10 диаметров заряда). В результате этого практически вся КС проникает в воду в разорванном на отдельные элементы состоянии. Поэтому необходимо учитывать возможность увеличения глубины проникания струи за счет инерционного движения воды, имеющего место в период между временем окончания проникания текущего элемента струи и временем начала проникания последующего элемента струи.
Для примера на рис. 13 представлены результаты численных расчетов проникания медных и стальных компактных ПЭ различной формы со скоростью 2 км/с в монолитные полубесконечные преграды из свинца, стали [1] и в комбинированную преграду, представляющую собой чередующие слой воды и стали [25].
а)
б)
в)
г)
д)
е)
Рис. 13. К иллюстрации проникания медных (а - д) и стальных (е) компактных ПЭ различной формы со скоростью монолитные полубесконечные преграды из свинца (а), стали (б - д) и комбинированную преграду, представляющую собой чередующие слой воды и стали
В работе [28] исследовалась возможность использования разработанного алгоритма применительно к решению задачи проникания КС в идеально проводящую преграду с поперечным магнитным полем. Для этого исходная система уравнений (1) - (11) дополнялись соотношениями, учитывающими эволюцию магнитного поля в среде при условии идеальной проводимости материала:
Ж
чл р
В Эу Ву ЭУ Ж
+
р дх р ду Ж
В
л
ху
V р У
Вх дУу + ВУ дУу
р дх р ду
2 2 2
где Вх, Ву - компоненты тензора магнитной индукции, причем В = (Вх) + (Ву) , а уравнения движения среды корректировались с учетом действия объемных электромагнитных сил:
Тхх = 0.5(Вх2 - В2 )/ц0 ; Ту = 0.5(Ву2 - В2 0; Тху = (ВхВу )/ц0.
Здесь Тх, Ту, Тху - компоненты тензора магнитных напряжений; ц0 = 4п10-7 Гн/м - магнитная постоянная. При этом показана возможность компрессии и роста интенсивности созданного в проводящей преграде магнитного поля при взаимодействии с высокоскоростными ударниками в виде КС. На основании численного моделирования проникания плоской КС в идеально проводящую преграду с поперечным магнитным полем установлено, что следствием компресс поля, происходящей в тонком слое материала преграды на границе контакта с КС, может быть проявление эффектов, способных оказать влияние на процесс проникания. Рассмотрен один из возможных эффектов, заключающийся в схло-пывании образованной струей каверны в преграде.
а)__б)
10 мкс Я 40 О 320 1591 10 мкс Й 40 О 319 VI 2172
В) Г)
10 мкс Я 40 О 312 1363 10 мкс И 40 О 309 ЯГ 2150
Рис. 14. Компьютерная визуализация поля взрыва в воде цилиндрического заряда ВВ с внутренними и внешними воздушными полостями конической формы: а) КВ с углом раствора (2а) 60°; а) коническая внешняя воздушная полость с углом раствора (2Р) 90° при отсутствии КВ (2а = 180°); в) 2а = 60°, 2р = 90°; г) 2а = 90°, 2р = 90°
а)
5 мкс 8 25 35 С") о 45 □ 150 о 231 ШГ 1235 2734
б)
5 мкс ■ 25 8" 0 131 О 210 1212 3025
Рис. 15. Компьютерная визуализация поля взрыва в воде цилиндрического заряда ВВ при различных способах инициирования: а) центральное точечное инициирование одновременно с двух торцов, б) многоточечное кольцевое инициирование по краям одновременно
с двух торцов
В ряде случаев рассмотренная методика пригодна для расчетов фугасных макетов. Так, например, в работе [18] анализировались демпфирующие возможности жидкостных оболочек в плане ослабления параметров поля взрыва в воздухе. В работах [20, 21] с использованием разработанной численной методики исследовалось фугасное действие под водой цилиндрического заряда ВВ (ТГ-40). Показано, что реализация направленного действия подводного фугасного взрыва осесимметричного цилиндрического заряда ВВ может быть обеспечена различными способами и конструктивными решениями как в осевом, так и в радиальном направлениях. В частности, при использовании зарядов, имеющих на одном конце либо кумулятивную выемку (КВ) в самом заряде, либо прилегающую к заряду воздушную полость различной конфигурации, либо и то и другое вместе (рис. 14). При этом существенное превышение максимального давления на фронте ударной волны (1.4 ... 1.7 раза) и удельного импульса (1.2 ... 1.3 раза) можно ожидать только в ближней зоне взрыва, на дистанциях, не превышающей 5 ... 10 приведенных радиусов заряда.
Рассмотренные схемы многоточечного инициирования заряда могут быть рекомендованы для повышения эффективности действия (до 40%) цилиндрических зарядов ВВ в радиальном (боковом) направлении (рис. 15). В случае же необходимости повышения эффективности действия цилиндрического фугасного заряда в осевом направлении целесообразно формировать внутри заряда ВВ кумулятивную выемку, комбинируя ее с воздушной полостью снаружи заряда непосредственно перед его торцом. Вместе с тем, указанные технические решения могут обеспечить достаточно ощутимый эффект только в ближней зоне взрыва, не превышающей дистанции 5 ... 10 приведенных радиусов заряда.
Повышение эффективности действия подводного взрыва зарядов ВВ ограниченной массы на преграды, расположенные на дистанциях свыше 10 приведенных радиусов может быть связано с применением КЗ, имеющих металлическую облицовку кумулятивной
выемки (для осевого направления), переходом на более мощные ВВ и совершенствованием схем многоточечного инициирования (для радиального направления).
РЭРС N01. NSWC
Рис. 16. Моделирование функционирования стандартного осколочного цилиндрического
макета №12 (RSFC)
В работах [13, 16, 17, 26] смоделирован процесс метания круглой металлической пластины с торцевой поверхности заряда ВВ при трех различных схемах инициирования. Рассмотрены волновые процессы в пластине и ПД. Установлены законы формирования углов разлета пластины. Рассмотрено влияние принятых величин предела текучести на результаты компьютерного моделирования процессов осесимметричного метания пластин, облицовок и оболочек. Показано, что это влияние является существенным и может приводить к изменению, как количественных характеристик процесса, так и его физической картины.
В работе [15] проведено моделирование процесса взрыва стандартного осколочного цилиндрического макета № 12 или RSFC (Russian Standard Fragmenting Cylinder) [1], наполненного составом A-IX-20 (рвв =
1.79 г/см , Двв = 7.96 км/с), см. рис. 16. Результаты расчета сравнивались с экспериментальными данными по углу разлета осколков и распределению массы осколков в различных угловых зонах. Представленные в ней результаты интересны также и тем, что получены на достаточно подробных разностных сетках, содержащих 102 тысячи и 625 тысяч ячеек, при этом размер пространственной разностной сетки составлял 0.65 мм и 0.14 мм соответственно. Полученные расчетные распределения масс и скоростей по двух градусным угловым зонам сравнивались с экспериментальными данными. При этом показано, что, начиная с некоторого предела, уменьшение размера разностной сетки уже не вносит заметных изменений в получаемые характеристики.
В работах [1, 12, 14, 15, 17, 26] представлены обширные расчетные исследования закономерностей расширения осесимметричных оболочек осколочно-фугасных снарядов под действием ПД. Рассмотрены волновые процессы, протекающие в ПД, стенках снарядов и придонной части корпуса. Установлены общие закономерности набора скоростей и формирования углов разлета оболочки. Разработана методология количественного определения характеристик фрагментации осколочных боеприпасов естественного и заданного дробления с использованием расчетных и экспериментальных данных [1]. При этом континуальная стадия расчета заканчивается процедурой построения законов распределения масс метаемой оболочки, как сплошной субстанции, и ее начальных скоростей по угловым зонам разлета в меридиональном направлении. На следующем этапе полученному угловому распределению с использованием эмпирических соотношений ставится в соответствие осколочный поток. При этом степень соответствия подобранного аналитического закона распределения реальному распределению осколков определяют обычно с применением критерия согласия Пирсона. Наконец, на заключительном этапе полученному распределению осколочного поля по меридиональному углу ставиться в соответствие годограф скоростей, производится оценка формы осколков и их эффективность действия по конкретной преграде.
Таким образом, разработана методика расчетного определения параметров взрывных полей, формируемых взрывными устройствами разного назначения, основанная на численном решении упругопластических задач механики сплошной среды в двумерной постановке. При этом к определяемым с достаточной для практики точностью характеристикам для КЗ можно отнести геометрические, массовые и кинематические параметры КС и ПЭ, осколочных макетов - временные распределения масс и скоростей оболочек по меридиональному углу разлета, для фугасных зарядов - характер затухания ударных волн в окружающей среде. Продемонстрированы возможности разработанной методики для моделирования быстропротекающих процессов в сравнении с экспериментальными данными и результатами расчетов по универсальным программным комплексам АКБУБ и АИТО-БУК разных версий. Представлены наиболее интересные результаты, полученные с ее использованием.
ЛИТЕРАТУРА
1. Физика взрыва: В 2 т. / В.И. Колпаков [и др.]; Под ред. Л.П. Орленко. Изд. 3-е, испр. М.: ФИЗМАТЛИТ, 2004. Т.2. 656 с.
2. Численные методы в задачах физики взрыва и удара: Учебник для втузов / В.И. Колпаков [и др.]; Науч. ред. В.В. Селиванов. 2-е изд., испр. М.: Изд-во МГТУ, 2000. 516 с.
3. Численные методы в задачах физики быстропротекающих процессов Учебник для втузов / В.И. Колпаков [и др.]; Науч. ред. В.В. Селиванов. М.: Изд-во МГТУ, 2006. 520 с.
4. Колпаков В.И., Ладов С.В., Рубцов А.А Математическое моделирование функционирования кумулятивных зарядов: Метод. указания. М.: Изд-во МГТУ им. Н.Э. Баумана, 1998. 38 с.
5. Мейдер Ч. Численное моделирование детонации. М.: Мир, 1985.384 с.
6. Уилкинс М. Расчет упруго-пластических течений // Вычислительные методы в гидродинамике. М.: Мир, 1967. С.212-263.
7. Ях К. Численный анализ задач классической и обратной кумуляции // ПМТФ. 1987. №2. С.123-129.
8. Высокоскоростное взаимодействие тел / В.М. Фомин [и др.] Новосибирск: Издательство СО РАН, 1999. 600 с.
9. Грязнов Е.Ф., Колпаков В.И., Уткин А.В. Экспериментальное исследование волновой стадии расширения стальной цилиндрической оболочки при осевой детонации заряда ВВ // Труды международной конференции «VII Харитоновские тематические научные чтения». Саров: РФЯЦ-ВНИИЭФ, 2005. С.565-570.
10. Оценка влияния технологических факторов на кинематические параметры удлиненного поражающего элемента кумулятивного заряда / В.И. Колпаков [и др.] // Труды международной конференции «IX Харитоновские тематические научные чтения». Саров: РФЯЦ-ВНИИЭФ, 2007. С.585-590.
11. Колпаков В.И., Кружков О.А., Шикунов Н.В. Математическое моделирование взрывного формирования УПЭ КЗ (тест П.И. Улякова) // Материалы VI Всерос. науч. конф. Томск: Том. гос. ун-т, 2008. С.253-254.
12. Одинцов В.А., Сидоренко Ю.М., Туборезов В.С. Моделирование процесса взрыва осколочно-фугасного снаряда с помощью двумерного гидрокода // Оборонная техника. 2000. № 1-2. С.49-55.
13. Моделирование процесса метания осколочной пластины с помощью двумерного гидрокода / В.А. Одинцов [и др.] // Оборонная техника. 2001. № 1-2. С.8-12.
14. Одинцов В.А., Сидоренко Ю.М. Угловое улавливание осколков // Оборонная техника. 2001. № 1-2. С.13-16.
15. Одинцов В.А., Сидоренко Ю.М. Моделирование процесса взрыва стандартного осколочного цилиндра с различной степенью детализации // Оборонная техника. 2001. №1-2. С.17-20.
16. Влияние положения точки инициирования на характеристики осколочного осевого потока / В.А. Одинцов [и др.] // Оборонная техника. 2002. №1-2. С.53-58.
17. Влияние величины предела текучести на результаты компьютерного моделирования процессов метания пластин, облицовок и оболочек / В.А. Одинцов [и др.] // Труды международной конференции «V Харитоновские тематические научные чтения». Саров: РФЯЦ-ВНИИЭФ, 2003. С.68-73.
18. Численная оценка эффективности действия жидкостных локализаторов взрыва в двухмерной постановке / В.И. Колпаков [и др.] // Двойные технологии. 2000. №2. С.5-10.
19. Колпаков В.И., Орленко Л.П., Рубцов А.А. Математическое моделирование направленного подводного взрыва осесимметричного заряда ВВ // Оборонная техника. 1995. №4. С. 25-28.
20. Колпаков В.И., Ладов С.В. Оценка степени влияния «направленности» подводного взрыва цилиндрического заряда ВВ при различных схемах инициирования // Труды международной конференции «V Харитоновские тематические научные чтения». Саров: РФЯЦ-ВНИИЭФ, 2003. С. 482-486.
21. Колпаков В.И., Ладов С.В. Численный анализ конструктивных схем зарядов, обеспечивающих направленное фугасное действие подводного взрыва // Оборонная техника. 2003. №3-4. С.49-55.
22. Колпаков В.И., Ладов С.В., Федоров С.В. Расчет формирования кумулятивных «ножа» удлиненного заряда с клиновидной выемкой // Оборонная техника. 1995. №1. С.24-29.
23. Колпаков В.И., Ладов С.В., Федоров С.В. Инженерная методика расчета действия КЗ с полусферическими и сегментными облицовками // Оборонная техника. 1999. №1-2. С.39-45.
24. Орленко Л.П., Бабкин А.В., Колпаков В.И. Задачи прикладной газодинамики: Результаты численного решения. М.: Изд-во МВТУ. 1988. 88 с.
25. Использование компактных поражающих элементов для пробития подводных преград / В.И Колпаков [и др.] // Оборонная техника. 2004. №1-2. С.67-71.
26. Компьютерное моделирование процессов осевого метания пластин, облицовок и оболочек при различных значениях предела текучести их материалов / В.А. Одинцов [и др.] // Оборонная техника. 2003. №3-4. С.61-69.
27. Колпаков В.И., Ладов С.В., Орленко Л.П. Методика расчета глубины проникания кумулятивной струи в воду // Оборонная техника. 2004. № 11. С.60-64.
28. Федоров С.В., Колпаков В.И., Бабкин А.В. Проникание плоской кумулятивной струи в идеально проводящую преграду с поперечным магнитным полем // Вестник МГТУ. Серия «Естественные науки». 2000. №2 (5). С.80-92.
electronic scientific and technical periodical
SCIENCE and EDUCATION
_EL № KS 77 -3()56'J..VaU421100025. ISSN 1994-jMOg_
Mathematical simulation of the explosive devices' performance 77-30569/334177
# 02, February 2012 Kolpakov V.I.
Bauman Moscow State Technical University Kolpakov@sm.bmstu.ru
Physico-mathematical statement of the problem of performance of explosive devices of different purposes, the fundamentals of numerical method for solving of this problem and description of algorithm for its implementation were considered in this article. Various examples of calculations were compared with experimental data; the possibilities of developed numerical technique were demonstrated.
Publications with keywords: method of calculation, numerical modelling, high explosive, pene-trator, model of explosive devices, explosive charge, shaped charge, the cumulative shell, shaped-charge jet, explosively formed projectiles, parameters field blast
Publications with words: method of calculation, numerical modelling, high explosive, penetra-tor, model of explosive devices, explosive charge, shaped charge, the cumulative shell, shaped-charge jet, explosively formed projectiles, parameters field blast
Reference
1. V.I. Kolpakov, et al., in: L.P. Orlenko (Ed.), Physics of explosion, 2 Volums, The vol. 2, Moscow, FIZMATLIT, 2004, 656 p.
2. V.I. Kolpakov, et al., in: V.V. Selivanov (Ed.), Numerical methods in problems of physics of explosion and shock, Moscow, Izd-vo MGTU - BMSTU Press, 2000, 516 p.
3. V.I. Kolpakov, et al., in: V.V. Selivanov (Ed.), Numerical methods in problems of physics of fast processes, Moscow, Izd-vo MGTU - BMSTU Press, 2006, 520 p.
4. Kolpakov V.I., Ladov S.V., Rubtsov A.A., Mathematical modeling of the functioning of shaped charges, Izd-vo MGTU im. N.E. Baumana - MGTU Press, 1998, 38 p.
5. Meider Ch, Numerical simulation of detonation, Moscow, Mir, 1985, 384 p.
6. Uilkins M, Calculation of elastic-plastic flows, Computational Methods in Fluid Dynamics, Moscow, Mir, 1967, pp. 212-263.
7. Iakh K., Numerical analysis of the problems of classical and inverse cumulative, Prikladnaja mehanika i tehnicheskaja fizika 2 (1987) 123-129.
8. V.M. Fomin, et al., High-speed interaction of bodies, Novosibirsk, Izdatel'stvo SO RAN, 1999, 600 p.
9.
10
11
12
13
14
15
16
17
18
19
20
21
22
Griaznov E.F., Kolpakov V.I., Utkin A.V., Experimental study of wave stage of expansion of steel cylindrical shells under axial detonation of the explosive charge, in: Proc. of the international conference "VII Kharitonov thematic scientific reading", Sarov, RFIaTs-VNIIEF, 2005, pp. 565-570.
V.I. Kolpakov, et al., Assessing the impact of technological factors on the kinematic parameters of an elongated striking element of the cumulative charge, in: Proc. of the international conference "IX Kharitonov thematic scientific reading", Sarov, RFIaTs-VNIIEF, 2007, pp. 585-590.
Kolpakov V.I., Kruzhkov O.A., Shikunov N.V., Mathematical modelling of explosive formation of an elongated striking element of the cumulative charge (test of PI. y^akob), in: Proceedings of VI All-Russian Scientific Conference, Tomsk, Tom. gos. un-t - TGU, 2008, pp. 253-254.
Odintsov V.A., Sidorenko Iu.M., Tuborezov V.S., Modeling of the process of explosion of high-explosive fragmentation shell through a two-dimensional hydrocode, Oboronnaia tekhnika 1-2 (2000) 49-55.
V.A. Odintsov, et al., Modeling of the process of throwing fragmentation plate with the help of two-dimensional hydrocode, Oboronnaia tekhnika 1-2 (2001) 8-12.
Odintsov V.A., Sidorenko Iu.M., Angular trapping of splinters, Oboronnaia tekhnika 1-2 (2001) 13-16.
Odintsov V.A., Sidorenko Iu.M., Modeling of the process of explosion of standard fragmentation cylinder, with varying degrees of detail, Oboronnaia tekhnika, 1-2 (2001) 17-20.
V.A. Odintsov, et al., The influence of point of initiation on the characteristics of fragmentation axial flow, Oboronnaia tekhnika 1-2 (2002) 53-58.
V.A. Odintsov, et al., The influence of size of limit of fluidity at the results of computer simulation of processes of throwing of plates, linings and shells, in: Proc. of the international conference "V Kharitonov thematic scientific reading", Sarov, RFIaTs-VNIIEF, 2003, pp. 68-73.
V.I. Kolpakov, et al., Numerical evaluation of the effectiveness of the liquid localizers explosion in the two-dimensional formulation, Dvoinye tekhnologii, 2 (2000) 5-10.
Kolpakov B.I., Orlenko L.P., Rubtsov A.A., Mathematical modeling of directional underwater explosion of an axisymmetric explosive charge, Oboronnaia tekhnika 4 (1995) 25-28.
Kolpakov V.I., Ladov S.V., Assessment of influence of "direction" of underwater explosion of cylindrical explosive charge at the various schemes of initiation, in: Proc. of the international conference "V Kharitonov thematic scientific reading", Sarov, RFIaTs-VNIIEF, 2003, pp. 482-486.
Kolpakov V.I., Ladov S.V., Numerical analysis of the structural schemes of the charges, providing directional high-explosive operation of underwater explosion, Oboronnaia tekhnika 3-4 (2003) 49-55.
Kolpakov V.I., Ladov S.V., Fedorov S.V, Calculation of the formation of cumulative "knife" of elongated charge with wedge-shaped groove, Oboronnaia tekhnika, 1 (1995) 24-29.
23. Kolpakov V.I., Ladov S.V., Fedorov S.V., Engineering method of calculation of the actions of the cumulative charge with hemispheric and segment facings, Oboronnaia tekhnika1-2 (1999) 39-45.
24. Orlenko L.P., Babkin A.V., Kolpakov V.I., Problems of applied gas dynamics: the Results of numerical solution, Moscow, Izd-vo MVTU - BMSTU Press. 1988, 88 p.
25. V.I Kolpakov, et al., The use of compact damaging elements for penetration of underwater barrier, Oboronnaia tekhnika 1-2 (2004) 67-71.
26. V.A. Odintsov, et al., Computer simulation of the processes of the axial throwing plates, linings and shells under various values of the limit of fluidity of their materials, Oboronnaia tekhnika, 3-4 (2003) 61-69.
27. Kolpakov V.I., Ladov S.V., Orlenko L.P., Methods of calculating of the depth of penetration of cumulative jet into the water, Oboronnaia tekhnika 11 (2004) 60-64.
28. Fedorov S.V., Kolpakov V.I., Babkin A.V., The penetration of the flat cumulative jet in a perfectly conductive barrier with transverse magnetic field, Vestnik MGTU. Ser. Estestven-nye nauki - Bulletin of BMSTU. Ser. Natural science 5 (2) (2000) 80-92.
