НАУКА и ОБРАЗОВАНИЕ

УДК 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 – реперные точки (маркеры)

Так как для обеих схем постановка задачи примерно идентична, рассмотрим ее на примере КЗ. Для этого будем полагать, что в начальный момент времени (t = 0) в точке 3 (точка инициирования) осуществляется подрыв заряда ВВ с начальной плотностью ρ_вв и теплотой взрывчатого превращения Q. От точки инициирования начинает распространяться фронт ДВ (кривая 5) со скоростью D_вв с образованием ПД (зона 6). С течением времени ДВ начинает отражаться от поверхностей корпуса (позиция 2) и КО (позиция 4), на которые действуют давления порядка 20 … 60 ГПа. Его величина зависит от свойств ВВ, угла подхода фронта ДВ к поверхности облицовки, материала и толщины облицовки. Под действием ПД кумулятивная облицовка начинает обжиматься с образованием КС или ПЭ. При этом для получения общих закономерностей или особенностей их формирования для конкретного КЗ, обусловленных формой облицовки, геометрией заряда, формой и месторасположением линзы и ТИ, физико-механическими свойствами используемого состава ВВ или материалов облицовки и корпуса, целесообразно использовать модель сжимаемой идеальной упругопластической среды с баротропным уравнением состояния [1 – 4]. Вязкими свойствами материала облицовки в процессе формирования КС или ПЭ можно пренебречь. Использование баротропного уравнения состояния позволяет избежать интегрирования уравнения энергии в системе соотношений, описывающих поведение взаимодействующих сред. Причем если при формировании КС можно ограничиться газодинамической моделью, то для определения формы и кинематических характеристик ПЭ, особенно на поздних стадиях движения, является важным учет проявления упругопластических свойств материала. Для зоны ПД течение считается изэнтропическим, так как даже отражение фронта ДВ от абсолютно жесткой поверхности не приводит к заметному росту энтропии. Распространение детонации рассматривается вне общей системы уравнений и задает границу области, охваченную течением. С учетом сделанных допущений система уравнений, описывающая двумерное течение в переменных Эйлера, имеет вид:

Здесь ρ – плотность; p – давление; t – текущее время; v_r, v_z – компоненты вектора скорости по осям 0r и 0z выбранной системы координат (рис. 1); σ_rr, σ_zz – нормальные и σ_rz – касательные напряжения; D_σrr, D_σrz, D_σzz – компоненты девиатора напряжений; D(…) / Dt – производная в смысле Яуманна, учитывающая поправки к составляющим девиатора напряжений, обусловленные поворотом фиксированного элемента среды как целого [2, 3, 5, 6]; G, Y – текущие значения модуля сдвига и динамического предела текучести материала среды, принимаемых в областях разгрузки в виде функциональных зависимостей от их начальных значений G₀,Y₀ и безразмерной плотности γ = ρ /ρ₀, причем при G₀ = 0 и Y₀ = 0 исходная система уравнений естественным образом трансформируется в уравнения идеального газа, описывающие поведение ПД; α – коэффициент симметрии (для плоского случая α = 0, для осесимметричного случая α = 1).

В приведенной системе уравнений соотношения (1) – (3) представляют собой, соответственно, законы сохранения массы и импульса; (4) – баротропные уравнения состояние взаимодействующих сред (КО, корпуса, линзы, ВВ и ПД), конкретизируемые ниже. Соотношения (5) связывают компоненты тензора полных напряжений с шаровой и девиаторной составляющими; (6) – закон Гука в дифференциальной форме; (7) – условие пластического течения Мизеса σ_i ≤ Y (σ_i – интенсивность напряжений). Если условие (7) нарушается, т.е. (f > 2Y²/3) и материал деформируемой среды находится в состоянии пластического течения, то компоненты D_σrr, D_σrz, D_σzzумножаются на множитель Y(2/3f)^0.5 (процедура приведения напряжений к кругу текучести [6]).

Для описания разрушения металлических элементов конструкции зарядов под действием ПД использовалась комбинация критериев откольной прочности p = – σ^*_p(σ^*_p – откольная прочность) и предельных пластических деформаций ε_р ≤ ε^*_р (ε^*_р – предельная пластическая деформация). Первый из них использовался на этапе нагружения элементов конструкции ДВ, второй – при последующем пластическом деформировании. Причем при r/r₀ < 1 и p = – σ^*_p – величина давления ограничивалась в соответствии с условиями [7, 8]

где

γ₁, γ₂ – постоянные материала, числовые значения которых, например, для стали по данным работ [9] составляют – γ₁ = 0.95, γ₂ = 0.90 для условий нагружения пластин иγ₁ = 0.975, γ₂ = 0.95 –для цилиндрических оболочек.

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

,

где ρ₀ – начальная плотность; А_b, n_b – эмпирические константы [1], причем при n_b ≠ 1 это уравнение называется ударной адиабатой Тэта, а при n_b = 1 оно вырождается в линейную баротропную зависимость, в которой A_b = K₀, где K₀ – модуль объемного сжатия.

Для продуктов детонации в качестве уравнения состояния рекомендуется использовать изоэнтропу в форме степенного двучлена

Коэффициенты изоэнтропы В, С и n (γ = 1.2…1.4) определяются по параметрам в точке Чепмена – Жуге [1 – 3]

 , ,

где   – давление, плотность и энергия на фронте ДВ, k – показатель адиабаты (для конденсированных ВВ k ≈ 3):

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

,

в котором коэффициент С_вв определяется из условия

.

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

Граничные условия для рассматриваемых задач в рамках идеализированных расчетных схем (рис. 1), задаются на участках поверхностей взаимодействующих сред. Так, например, поверхности КО и корпуса, контактирующих с воздухом, и поверхность продуктов детонации, истекающих из корпуса, допустимо рассматривать как свободные от действия внешних поверхностных сил. То есть, пренебрегая силами атмосферного давления – σ_ijn^j = 0 (pn^j = p_nj = 0 для ПД), где n^j – вектор единичной нормали. На поверхностях КО, линзы и корпуса, контактирующих с ПД, накладываются ограничения на скорости движения индивидуальных точек в соответствии с условием непроницаемости

а также на напряженное состояние, реализующееся в этих точках в соответствии с третьим законом Ньютона

При формулировке граничных условий на оси симметрии (ось 0z) необходимо учитывать, что при r = 0 частицы среды движутся только в осевом направлении (v_r = 0), а осевые ускорения этих частиц должны быть ограничены. Из уравнений движения (2) – (3) следует, что это может быть реализовано только при отсутствии касательных напряжений на оси симметрии (σ_rz = 0). Давление, плотность и энергия во фронте ДВ задаются равными параметрам в точке Чепмена – Жуге.

Выделение контактных разрывов типа ПД – корпус или ПД – КО осуществлялось методом «концентраций» [2, 3], разработанного В.В. Кореньковым. Для этого в дополнение к основным характеристикам течения (ρ, p, v_r, v_z, D_σrr, D_σrz, D_σzz) для неоднородной системы ПД (ВВ) – КО (корпус) определяют массовую концентрацию веществ следующим образом:

Здесь M₁ – масса ВВ, ПД или воздуха в ячейке, а M – масса ячейки. Таким образом, для однородных ячеек, содержащих только первую компоненту w = 1, для однородных ячеек со второй компонентой w = 0 и, наконец, для смешанных ячеек 0 < w < 1.

Локализация контактных границ осуществляется из анализа текущего распределения концентрации взаимодействующих веществ, которое определяется законом сохранения концентрации

Давление в смешанной ячейке вычисляется из условия аддитивности удельных объемов и равенства давлений содержащихся в ней компонент

где ρ, ρ_{w = 0}, ρ_{w = 1} – текущая плотность смешанной ячейки и плотности содержащихся в ней компонент; p₀(ρ_{w = 0}) – ударная адиабата материала корпуса; p₁(ρ_{w = 1}) – изэнтропа ПД или уравнение состояния ВВ. По значениям поля концентраций определяется вид контактной границы и на его основе рассчитывается поток из смешанной ячейки. На разрывах типа ПД – КО или ПД – корпус в процессе решения обеспечиваются условия равенства нормальных составляющих скоростей и напряжений. Выделения границ ПД – воздух не производится.

Начальные условия конкретной задачи задаются распределением параметров ρ, p, w, v_r и v_z в поле течения. Компоненты напряжений принимаются равными нулю.

Для получения более полной информации о процессах схлопывания КО или метания корпуса последние обычно маркируют реперными точками — маркерами или трассерами (позиция 8 на рис.1), в которых дополнительно вычисляют параметры текущего состояния среды P_m = {ρ, p, v_r, v_z …} по формуле

Здесь k = 1 2, 3, 4; A_k – площади, определяемые текущим положением маркера m [2 – 5], а P_m(k) – параметры состояния среды в ячейках (k), окружающих маркер.

Система уравнений (1) – (11) интегрировалась модифицированным методом крупных частиц, подробно изложенным в следующем разделе (разд. 2), на программном комплексе KOLDUN [4], позволяющем проводить вычисления во всей области интегрирования на неподвижной (эйлеровой) сетке без предварительного анализа течения.

2. МЕТОД ЧИСЛЕННОГО РЕШЕНИЯ

Рис. 2. Структура неподвижной (эйлеровой) разностной сетки:

а) область интегрирования задачи (0 ≤ z ≤ Z, 0 ≤ r ≤ R); б) индексация ячеек разностной сетки («фиктивные» ячейки выделены серым цветом); Г₁, Г₂, Г₃, Г₄ – границы расчетной области

Дискретизация задачи. Пусть в области интегрирования задачи (0 ≤ z ≤ Z, 0 ≤ r ≤ R) происходит плоское или осесимметричное неустановившееся движение неоднородной среды, описываемое уравнениями (1) – (10), при заданных начальных и граничных условиях. Область интегрирования разбивается на некоторое число прямоугольных ячеек со сторонами Δr = R/m и Δz = Z/n, которые образуют неподвижную (эйлерову) разностную сетку (рис. 2, а). Значения целых чисел i = 2, 3, … , m+1 и j = 2, 3, … , n+1 обозначают центры ячеек разностной сетки (рис. 2, б), где m и n – число ячеек по осям r и z соответственно (для конкретного заряда это число обычно задают таким образом, чтобы по минимальной толщине оболочки оно было бы равно не менее 4 … 5 целых ячеек, причем Δr/Δz = 0.8 … 1.2). В случае плоской симметрии отдельные ячейки представляют собой прямоугольники, в случае цилиндрической симметрии — тороиды. В точках (i, j) определяются средние значения параметров течения среды — компонент вектора массовой скорости, плотности, давления, концентрации, компонент напряжений.

Область интегрирования слева ограничена осью симметрии или жесткой стенкой (Г₁, на рис. 2, а); снизу, справа, сверху — открытыми поверхностями (Г₂, Г₃ и Г₄), через которые среда может вытекать или втекать. Чтобы не нарушать единообразия вычислений для граничных ячеек, вдоль всех упомянутых границ вводятся слои «фиктивных» ячеек (см. рис. 2, б), параметры которых определяют из смежных ячеек, как и в методе «частиц в ячейках» или в методе «крупных частиц» [2, 3]. Например, для оси симметрии или жесткой стенки (поверхность Г₁) имеем:

где параметры помеченные тильдой определяются в конце I-го этапа вычислений, а параметры с верхним индексом (n+1) — в конце III-го этапа [2, 3].

Для открытой поверхности Г₂:

Для открытой поверхности Г₃:

Для открытой поверхности Г₄:

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

В процессе отработки разностной схемы было рассмотрено несколько вариантов геометрического расщепления исходной системы уравнений (1) – (10) по пространственным координатам. Здесь излагается наиболее простой из них, в котором вместо соотношений (1) – (3) и (10) используются две подсистемы [2, 3]:

Последовательное интегрирование уравнений (17) и (18) с шагом τ = Δtⁿ/2 по каждому из направлений r и z в предположении неизменности компонент девиатора напряжений на текущем временном слое (слой n) позволяет получить значения состояния среды ρⁿ⁺¹, pⁿ⁺¹, v_rⁿ⁺¹, v_zⁿ⁺¹, wⁿ⁺¹ на новом временном слое tⁿ⁺¹= tⁿ + Δtⁿ. При этом в обоих случаях устойчивость вычислений обеспечивается выбором шага интегрирования по времени

где i, j – индексы текущей ячейки; V_i,j, C_i,jⁿ, (v_r)_i,jⁿ, (v_z)_i,jⁿ – объем, скорость звука материала и компоненты массовой скорости текущей в ячейке; K_r = 0.1 … 0.35 – число Куранта.

На следующем этапе по вновь полученным характеристикам течения среды    рассчитываются значения производных      и на основании соотношений (6), (7) определяются обновленные компоненты девиатора напряжений   

Далее излагается разностная схема для рассмотренного варианта расщепления исходной системы уравнений. Для этого будем полагать, что в текущий момент времени tⁿ известны все величины потока двухкомпонентной среды в центрах ячеек эйлеровой сетки        Получим соотношения для вычисления перечисленных параметров внутри и на границе расчетной области в момент времени tⁿ⁺¹= tⁿ + Δtⁿ, где Δtⁿ – шаг по времени, определяемый соотношением (19).

Рис. 3. Индексация границ ячейки на неподвижной разностной сетке

Этап I (эйлеров). На этом этапе вычислений среда предполагается «замороженной», т.е. отсутствует массообмен субстанцией между смежными ячейками разностной сетки. Изменение параметров течения среды происходит за счет напряжений, действующих в элементарном фиксированном объеме. В этом случае уравнения (17) и (18) принимают вид:

Заменим уравнения (20) разностными. Для этого определим конечно-разностные аналоги производных типа  и  следующим образом:

Здесь f_i,jⁿ – значение некоторой сеточной функции в центре ячейки (i, j) на n-ом временном слое (рис. 3); f₁ⁿ, f₂ⁿ, f₃ⁿ, f₄ⁿ – значения произвольной сеточной функции на левой, нижней, правой и верхней границах ячейки (i, j) соответственно, т.е.

V_i,j и S₁, S₂, S₃, S₄ – соответственно объем и площади левой (1), нижней (2), правой (3) и верхней (4) граней (i, j)-ой ячейки, определяемые формулами, приведенными в табл. 1.

С учетом соотношений (21) – (22) и величины шага интегрирования по времени τ = Δtⁿ/2 определяются промежуточные значения массовых скоростей течения среды  и  в ячейке (i, j):

Формулы, определяющие объем и площадь граней фиксированной ячейки в плоском и осесимметричном случаях

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

Рис. 4. К определению потока массы среды через границу между смежными ячейками (i, j) и (i+1, j): а) ; б) ; Д – донорская ячейка; А – акцепторная ячейка

В начале определим поток массы среды с одинаковой концентрацией (w = 0 или w = 1) через границу между ячейками. Например, для правой границы произвольной ячейки (рис. 4) он определяется как

Здесь ρ_i,j, ρ_i+1,j – плотности сред, содержащейся в смежных ячейках; S_i+1/2,j –площадь правой границы ячейки (i, j), через которую производится массообмен между ячейками (i, j) и (i +1, j); . Согласно выражению (23) поток массы среды с концентрацией w = 1 (Δm₁) через ту же границу ячейки рассчитывается по формуле (24)

Обратимся теперь случай, когда одна из смежных ячеек является смешанной (0 < w < 1), а другая — однородной (w = 1 или w = 0). Для примера будем рассматривать поток массы среды через правую границу смешанной ячейки (i, j) (см. рис. 4). При этом назовем ячейку, из которой происходит истечение среды, донорской (Д), а другую, в которую направлен поток, — акцепторной (А).

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

где

Здесь M₁, M₀ – первая и вторая составляющие потока массы среды через правую границу смежной ячейки (i, j), соответствующие концентрации w = 1 или w = 0. При этом последовательность расчета величин M₀ и M₁ определяется значением концентрации вещества, содержащегося в однородной ячейке. Здесь α_i,j = (w_i,jρ_i,j)/(ρ_w=1)_i,j, (1 – α_i,j) – объемные концентрации первой и второй компонент среды в (i, j)-ой ячейке; d(Δm₀)_L, d(Δm₁)_L, L = [(i, j), (i+1, j)] – потоки массы двухкомпонентной среды, определяемые текущими значениями плотностей (ρ_w=0)_L и (ρ_w=1)_L, а именно:

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

Рассмотрим вариант такой корректировки для потока M₁ между смешанными ячейками (i, j) и (i+1, j) при условии  Для этого обозначим массу вещества с признаком w = 1, находящуюся в донорской ячейке (i, j), через (Δm₁)_i,j. Если (Δm₁)_i,j стала меньше перетекающей массы M₁, определяемой соотношением (25), значит, все вещество с признаком w = 1 вытекло из ячейки через правую границу. При этом время истечения составляет Δt₁ = Δtⁿ (Δm₁/M₁)_i,j. За промежуток времени Δt₂ = Δtⁿ – Δt₁  через площадь α_i+1,j S_i+1/2,j в ячейку (i +1, j) будет втекать второй компонент среды с признаком w = 0. С учетом сказанного составляющая M₁ потока массы среды пересчитывается следующим образом:

Сделав аналогичные выкладки для M₁ при  и для M₀ при  или , окончательно можно получить следующие соотношения:

Значения нижних индексов L, N и константы ζ в соотношениях (27) и (28) в зависимости от направления движения среды  или  приведены в табл. 2.

Значения нижних индексов и константы в уравнениях (27) и (28)

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

III Этап (заключительный). На этом этапе в пределах объемов ячеек эйлеровой координатной сетки определяются параметры состояния среды на временном слое n+1 с учетом переноса ее компонентов через границы ячеек:

Здесь k = r, z; ζ_L = 1, если поток массы среды через границу L = [(i – 1/2, j), (i + 1/2, j)], направлен вовнутрь ячейки (i, j), и ζ_L = 0 в противном случае.

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

Компоненты девиатора напряжений. Отличительной особенностью уравнений (6), определяющих изменение компонент девиатора напряжений в неподвижной системе координат по отношению к лагранжевой форме записи аналогичных уравнений, является наличие так называемых конвективных составляющих полных производных величин   . Для получения их конечно-разностного аналога используется подход, предложенный Хагеман, Уолшем, Мейдером [5], в котором компоненты девиатора напряжений на временном слое n+1 определяются по величинам    и промежуточным значениям приращений   , взятых с соответствующими весами. При этом последовательность вычислений следующая.

Во-первых, с учетом граничных условий (13) – (16) в каждой ячейке определяются производные вида

и далее

Во-вторых, вычисляются приращения девиатора напряжений:

;

;

;

где

В-третьих, определяются значения составляющих девиатора напряжений на временном n+1 слое с учетом конечно-разностного аналога конвективных производных величин   взятых с соответствующими весами

  и   

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

Рис. 5. К определению разностной записи конвективной производной:

a) δr_i,j < 0, δz_i,j < 0; б) δr_i,j ≥ 0, δz_i,j ≥ 0; в) δr_i,j < 0, δz_i,j ≥ 0; г) δ_i,j ≥ 0, δz_i,j < 0

Если величины δr_i,j и δz_,j, отсчитываемые от центра ячейки (i, j), отрицательны (рис. 5, а), то имеем

где

При i = 2 (см. рис. 2) —

при j = 2 —

при i = 2 и j = 2 —

Если δr_i,j ≥ 0 и δz_i,j ≥ 0 (рис. 5, б), то

При i = m+1 (см. рис. 2) —

, ;

при j = n +1—

при i = m+1 и j = n +1—

Если δr_i,j < 0 и δz_i,j ≥ 0 (рис.5, в), то

При i = 2 (см. рис. 2) —

при j = n +1—

при i = 2 и j = n +1 —

Если δr_i,j ≥ 0 и δz_i,j < 0 (рис 5, г), то

При i = m+1 (см. рис. 2) —

при j = 2 —

при i = m+1 и при j = 2 —

В-четвертых, реализация условия пластического течения обеспечивается выполнением процедуры приведения напряжений на круг текучести [2, 3]. Для этого в каждой ячейке рассчитываются значения

В случае f_i,jⁿ⁺¹ > (2/3)Y² производится корректировка компонент девиатора напряжений по формулам

где скорректированные значения компонент девиатора напряжений определяются следующим выражением

.

Реперные точки. Для получения более полной информации о процессах схлопывания облицовки или метания корпуса последние обычно маркируют реперными точками (маркерами или трассерами), в которых по формуле (12) дополнительно вычисляют параметры текущего состояния среды, см. рис. 6.

Рис. 6. К определению параметров течения среды в маркере

3. ВЕРИФИКАЦИЯ РАЗРАБОТАННОЙ МЕТОДИКИ РАСЧЕТА

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

Сравнение с экспериментальными данными и альтернативными методами расчета. Особое внимание уделено сопоставлению результатов расчетов с данными рентгеноимпульсной съемки. Для этого использовались лабораторные СФЗ с внутренним диаметром корпуса 50 мм (макет №1) и 64.1 мм (макет №2), показанные на рис. 7, и осколочные цилиндрические макеты с диаметром шашки 20 мм, показанные на рис. 8. Кумулятивные заряды снаряжались литьевыми составами ТГ-50 (ρ_вв = 1.65 г/см³, D_вв = 7.7 км/с, Q = 4.61 МДж/ кг) и ТГ-40 (ρ_вв = 1.68 г/см³, D_вв = 7.85 км/с, Q = 4.61 МДж/ кг), осколочные макеты — прессованными шашками из флегматизированного гегсогена (ρ_вв = 1.65 г/см³, D_вв = 8.1 км/с, Q = 4.82 МДж/ кг) или флегматизированного октогена (ρ_вв = 1.75 г/см³, D_вв = 8.6 км/с, Q = 5.23 МДж/ кг). Кумулятивные облицовки СФЗ изготавливались из пластичной стали 11ЮА, физико-механические параметры которой принимались следующими: плотность – ρ₀ = 7.85 г/см³; модуль сдвига – G₀ = 80 ГПа; динамический предел текучести – Y₀ = 0.5 … 0.8 ГПа, откольная прочность – σ^*_p = 1.65 ГПа. Для корпусов кумулятивных макетов применялась сталь 45Х (G₀ = 77…81 ГПа; Y₀ = 1.05 … 1.1 ГПа, σ^*_p ≈ 2.3 ГПа), а для цилиндрических оболочек осколочных макетов – сталь 20 (G₀ = 76…79 ГПа; Y₀ = 0.43 ГПа, σ^*_p = 1.65 ГПа). В конструктивном отношении осколочные макеты отличались между собой лишь формой верхнего дна.

Рис. 9. Cравнение экспериментальной и расчетной формы поражающего элемента, образованного взрывом СФЗ №1, на начальном этапе деформирования (20 мкс, 25 мкс и 33 мм): а) данные рентгеноимпульсной съемки, б) расчет по программе KOLDUN, в, г) расчет по программе AUTODYN; г) КО без отверстия в купольной части

Результаты отдельных выполненных тестов представлены на рис. 9 – 12. На рис. 9 показано сравнение расчетной формы ПЭ (б, в) с данными рентгеноимпульсной съемки (а) на начальном этапе деформирования поражающего элемента (20 мкс, 25 мкс и 33 мкс), сформированного СФЗ №1, конструктивная схема которого схематично показана на рис. 7, а. При этом рис. 9, б соответствует расчету по программе KOLDUN [2 – 4], а рис. 9, в – AUTODYN. Для сравнения на рис. 9, г показан расчет по программе AUTODYN аналогичного заряда с облицовкой без отверстия в купольной части.

Отличительной особенностью СФЗ №1 является то обстоятельство, что используемая в нем кумулятивная облицовка имеет достаточно сложную и не традиционную форму. Ее наружная (обращенная к ВВ) поверхность представляет собой сферический сегмент с радиусом образующей 36 мм, а внутренняя (обращенная к преграде) поверхность очерчивается образующей в виде параболической кривой r = (39.5z)^0.5 +5 относительно локальной системы координат, показанной на рис. 7, а. В купольной части КО имеется отверстие диаметром 4 мм. Как было показано в [10, 11] при взрыве такого заряда образуется удлиненный ПЭ, который разрушается при деформировании в процессе движения к преграде. Поэтому СФЗ №1, получивший впоследствии название заряд Улякова П.И. [11], удобно использовать для моделирования процесса разрушения УПЭ из стали 11ЮА. Как видно из представленных иллюстраций, полученные расчетным путем формы и кинематические характеристики ПЭ на начальной стадии взрыва (до ~35 мкс) при использовании всех представленных программ расчета имеют хорошее соответствие реально протекающему процессу. Неясностей физического толка для математического описания начальной стадии взрыва СФЗ не возникает.

Рис. 10. Сравнение расчетной формы ПЭ с данными рентгеноимпульсной съемки: а) заряд №2 (с корпусом); б) заряд №2 (без корпуса)

На рис. 10 представлены результаты расчетов СФЗ №2 в сравнении с экспериментальными данными. Как уже отмечалось выше, в нем использовалась сегментная облицовка постоянной толщины 2.2 мм. Ее внутренняя и наружная поверхности представляли собой сферические сегменты с радиусами кривизны 54 мм и 56,2 мм соответственно.

Рис. 10 иллюстрирует степень влияния корпуса (поз. 2 на рис. 7) на геометрические характеристики формируемого взрывом ПЭ. Для этого показаны характерные стадии деформирования поражающих элементов, образованных СФЗ №2 с корпусом (рис. 10, а) и без корпуса (рис. 10, б). При этом сравниваются формы ПЭ, полученные расчетным путем при Y₀ = 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 мкс а) расчет по программе KOLDUN; б) расчет по программе LS-DYNA; в) эксперимент

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

Возможности расчетной методики. Программный комплекс «КОЛДУН» или «KOLDUN» [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].

Рис 12. Сравнение данных расчета, полученных по программам KOLDUN и AUTODYN (SPH и Euler), функционирования осколочного макета №2 с толщиной стенки 1 мм, наполненного А-91 на моменты времени 10 мкс и 15 мкс: а, г) расчет по программе KOLDUN; б, д) расчет по программе AUTODYN (SPH); в, е) расчет по программе по AUTODYN (Euler)

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

В работах [1 – 3, 26] представлены результаты расчетных исследований критических условий образования кумулятивных струй, а также влияния типа симметрии (плоская, осевая) облицовки и сжимаемости ее материала на скорость КС для установившейся фазы истечения. Они проводилось в рамках идеализированных расчетных схем косого соударения металлических пластин и конических облицовок с жесткой стенкой или осью симметрии в диапазоне скоростей соударения 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) дополнялись соотношениями, учитывающими эволюцию магнитного поля в среде при условии идеальной проводимости материала:

; ,

где B_x, B_y – компоненты тензора магнитной индукции, причем B² = (B_x)² + (B_y)², а уравнения движения среды корректировались с учетом действия объемных электромагнитных сил:

;    ;   .

Здесь Т_х, Т_у, Т_xу – компоненты тензора магнитных напряжений; μ₀ = 4π·10^–7 Гн/м – магнитная постоянная. При этом показана возможность компрессии и роста интенсивности созданного в проводящей преграде магнитного поля при взаимодействии с высокоскоростными ударниками в виде КС. На основании численного моделирования проникания плоской КС в идеально проводящую преграду с поперечным магнитным полем установлено, что следствием компресс поля, происходящей в тонком слое материала преграды на границе контакта с КС, может быть проявление эффектов, способных оказать влияние на процесс проникания. Рассмотрен один из возможных эффектов, заключающийся в схлопывании образованной струей каверны в преграде.

Рис. 14. Компьютерная визуализация поля взрыва в воде цилиндрического заряда ВВ с внутренними и внешними воздушными полостями конической формы: а) КВ с углом раствора (2α) 60°; а) коническая внешняя воздушная полость с углом раствора (2β) 90° при отсутствии КВ (2α = 180°); в) 2α = 60°, 2β = 90°; г) 2α = 90°, 2β = 90°

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

В ряде случаев рассмотренная методика пригодна для расчетов фугасных макетов. Так, например, в работе [18] анализировались демпфирующие возможности жидкостных оболочек в плане ослабления параметров поля взрыва в воздухе. В работах [20, 21] с использованием разработанной численной методики исследовалось фугасное действие под водой цилиндрического заряда ВВ (ТГ-40). Показано, что реализация направленного действия подводного фугасного взрыва осесимметричного цилиндрического заряда ВВ может быть обеспечена различными способами и конструктивными решениями как в осевом, так и в радиальном направлениях. В частности, при использовании зарядов, имеющих на одном конце либо кумулятивную выемку (КВ) в самом заряде, либо прилегающую к заряду воздушную полость различной конфигурации, либо и то и другое вместе (рис. 14). При этом существенное превышение максимального давления на фронте ударной волны (1.4 … 1.7 раза) и удельного импульса (1.2 … 1.3 раза) можно ожидать только в ближней зоне взрыва, на дистанциях, не превышающей 5 … 10 приведенных радиусов заряда.

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

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

Рис. 16. Моделирование функционирования стандартного осколочного цилиндрического макета №12 (RSFC)

В работах [13, 16, 17, 26] смоделирован процесс метания круглой металлической пластины с торцевой поверхности заряда ВВ при трех различных схемах инициирования. Рассмотрены волновые процессы в пластине и ПД. Установлены законы формирования углов разлета пластины. Рассмотрено влияние принятых величин предела текучести на результаты компьютерного моделирования процессов осесимметричного метания пластин, облицовок и оболочек. Показано, что это влияние является существенным и может приводить к изменению, как количественных характеристик процесса, так и его физической картины.

В работе [15] проведено моделирование процесса взрыва стандартного осколочного цилиндрического макета № 12 или RSFC (RussianStandardFragmentingCylinder) [1], наполненного составом А-IX-20 (ρ_вв = 1.79 г/см³, D_вв = 7.96 км/с), см. рис. 16. Результаты расчета сравнивались с экспериментальными данными по углу разлета осколков и распределению массы осколков в различных угловых зонах. Представленные в ней результаты интересны также и тем, что получены на достаточно подробных разностных сетках, содержащих 102 тысячи и 625 тысяч ячеек, при этом размер пространственной разностной сетки составлял 0.65 мм и 0.14 мм соответственно. Полученные расчетные распределения масс и скоростей по двух градусным угловым зонам сравнивались с экспериментальными данными. При этом показано, что, начиная с некоторого предела, уменьшение размера разностной сетки уже не вносит заметных изменений в получаемые характеристики.

В работах [1, 12, 14, 15, 17, 26] представлены обширные расчетные исследования закономерностей расширения осесимметричных оболочек осколочно-фугасных снарядов под действием ПД. Рассмотрены волновые процессы, протекающие в ПД, стенках снарядов и придонной части корпуса. Установлены общие закономерности набора скоростей и формирования углов разлета оболочки. Разработана методология количественного определения характеристик фрагментации осколочных боеприпасов естественного и заданного  дробления с использованием расчетных и экспериментальных данных [1]. При этом континуальная стадия расчета заканчивается процедурой построения законов распределения масс метаемой оболочки, как сплошной субстанции, и ее начальных скоростей по угловым зонам разлета в меридиональном направлении. На следующем этапе полученному угловому распределению с использованием эмпирических соотношений ставится в соответствие осколочный поток. При этом степень соответствия подобранного аналитического закона распределения реальному распределению осколков определяют обычно с применением критерия согласия Пирсона. Наконец, на заключительном этапе полученному распределению осколочного поля по меридиональному углу ставиться в соответствие годограф скоростей, производится оценка формы осколков и их эффективность действия по конкретной преграде.

Таким образом, разработана методика расчетного определения параметров взрывных полей, формируемых взрывными устройствами разного назначения, основанная на численном решении упругопластических задач механики сплошной среды в двумерной постановке. При этом к определяемым с достаточной для практики точностью характеристикам для КЗ можно отнести геометрические, массовые и кинематические параметры КС и ПЭ, осколочных макетов – временные распределения масс и скоростей оболочек по меридиональному углу разлета, для фугасных зарядов – характер затухания ударных волн в окружающей среде. Продемонстрированы возможности разработанной методики для моделирования быстропротекающих процессов в сравнении с экспериментальными данными и результатами расчетов по универсальным программным комплексам ANSYS и AUTODYN разных версий. Представлены наиболее интересные результаты, полученные с ее использованием.

ЛИТЕРАТУРА

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.        Оценка влияния технологических факторов на кинематические параметры удлиненного поражающего элемента кумулятивного заряда / В.И. Колпаков [и др.] // Труды международной конференции «IХ Харитоновские тематические научные чтения». Саров: РФЯЦ-ВНИИЭФ, 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.        Колпаков B.И., Орленко Л.П., Рубцов А.А. Математическое моделирование направленного подводного взрыва осесимметричного заряда ВВ // Оборонная техника. 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.