Моделирование формирования порохового заряда, состоящего из зерненых пороховых элементов

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

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

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

где ρ _j - плотность газа для части заряда с ин­дексом j;

m - пористость смеси;

S - переменная площадь поперечного се­чения каморы и канала ствола;

t - время;

v, w - скорость движения газа и твердой фазы в канале ствола соответственно; х - координата;

ρ = ρ₁ +ρ₂ - суммарная плотность смеси продуктов горения;

τ _w1, τ _w2 - сила взаимодействия, приходяща­яся на единицу объема, вызванная разностью скоростей между фазами для первой и второй частей заряда;

G₁, G₂ - газоприход с поверхности пороха в единице объема за единицу времени;

П_с - периметр ствола;

q_c - тепловой поток, направленный на по­верхность канала ствола;

τ_c - сила трения газа о поверхность ствола, приходящаяся на единицу площади;

ε - внутренняя энергия единицы массы пороховых газов;

Q₁, Q₂ - теплотворная способность (потен­циал) пороха;

р - давление;

α - коволюм;

θ = k — 1, k = c_p / c_v - показатель адиабаты; c_p, c_v - теплоемкость газа при постоянном давлении и объеме;

a_j - счетная концентрация зерненых поро­ховых элементов;

δ - плотность материала пороха;

ψ_j - относительная доля сгоревшего по­роха для первой и второй частей заряда;

Λ₀, S₀ - начальный объем и поверхность порохового зерна;

σ(ψ) - отношение текущей поверхности горения к первоначальной;

U_k - линейная скорость горения пороха;

æ_р - коэффициент формы порохового эле­мента;

e₁ - половина начальной толщины горяще­го свода.

Индекс j в уравнениях обозначает номер части заряда для комбинированных зарядов.

Начальные и граничные условия имеют следующий вид [1, 2].

1. Начальные условия (t = 0):

где р₀ - давление форсирования;

ω _j - масса j-го полузаряда;

W_j - объем каморы, занятый j-м полуза- рядом;

m₀ - начальная пористость заряда;

W_км - объем каморы орудия;

U_{1 j} - скорость горения пороха при атмо­сферном давлении (единичная скорость горе­ния);

v _j - показатель степени в уравнении для скорости горения;

f_j - сила пороха.

2. Граничные условия.

При х = 0 (у дна каморы)

v = w = 0;

при х = х_сн (у дна снаряда)

v = v_сн

где х_сн - координата дна снаряда в данный момент времени;

v_сн - скорость снаряда в данный момент времени.

При расчете внутрибаллистических па­раметров в уравнениях газодинамического подхода используется такая величина, как по­ристость - m. Пористость - это объем пустот в единице объема, занятого пороховыми эле­ментами. Таким образом, пороховые элементы будут занимать объем 1 - m. По мере выгорания пороха пористость будет увеличиваться.

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

где Δ - плотность заряжания;

δ - плотность пороха.

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

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

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

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

Цилиндр заполняется пороховыми эле­ментами. Каждый пороховой элемент возникает на верхнем срезе цилиндра в произвольной точ­ке оси OX₀ - точка возникновения элемента со­ответствует расположению его центра тяжести в начальный момент времени. Затем он падает вниз до тех пор, пока не достигает дна цилин­дра или других пороховых элементов, располо­женных под ним. После этого на верхнем срезе цилиндра возникает новый пороховой элемент, который имеет произвольную ориентацию в пространстве, т. е. его ось симметрии может составлять любой угол α с осью OX₀.

После того как падающий пороховой эле­мент достигает слоя нижележащих элементов, он останавливается или продолжает движение в зависимости от положения центра тяжести элемента относительно точки опоры. При рас­чете движения порохового элемента кроме ос­новной системы координат OX₀Y₀ для каждого падающего элемента строится еще одна систе­ма координат Oxyz (рис. 2, ось z на рисунке не показана), связанная с пороховым элементом. За начало координат принята точка O - центр тяжести элемента.

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

qa_x = G sin α- F_тp;

где q - масса порохового элемента;

a_x, a_y - проекции ускорения элемента на оси х и у;

α - угол поворота порохового элемента относительно оси OX₀;

I_z - момент инерции порохового элемента относительно оси z;

ε - угловое ускорение;

g - ускорение свободного падения;

f - коэффициент трения.

Если проекции ускорений представить в виде вторых производных координаты по вре­мени, можно получить следующую систему дифференциальных уравнений:

где х и у - координаты точки опоры в системе координат Oxyz; t - время.

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

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

1. Переход к свободному падению эле­мента (рис. 3, а).

В процессе движения порохового эле­мента он может потерять контакт с другими пороховыми элементами, это возможно при выполнении условия

где a = l или a = d в зависимости от располо­жения в пространстве падающего элемента.

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

2. Прекращение движения элемента (рис. 3, б).

В этом случае система уравнений изменя­ется с учетом появления второй точки контакта:

Если точки соприкосновения расположе­ны на двух взаимно перпендикулярных сторо­нах порохового элемента, то поступательное движение прекратится. При этом поворот от­носительно оси z может продолжаться в зави­симости от соотношения между величинами Δr, h и α. Можно доказать, что прекращение вращательного движения элемента произойдет в случае выполнения условия:

Если условие выполняется, то расчет дви­жения данного элемента прекращается, генери­руется новый пороховой элемент. Если условие не выполняется, то продолжается расчет пово­рота данного элемента относительно оси z.

3. Прекращение вращательного движения при продолжении поступательного (рис. 3, в).

Если точки соприкосновения расположе­ны на одной стороне порохового элемента и линия действия силы тяжести проходит между точками опоры, то поступательное движение элемента вдоль оси х будет продолжаться, а вра­щение элемента прекратится. Для доказатель­ства этого утверждения достаточно составить два уравнения моментов относительно первой и второй точек опоры; решая их совместно, можно получить значение ε = 0. Для этого слу­чая справедлива следующая система уравнений:

В случае равенства нулю углового уско­рения угловая скорость поворота вокруг оси z остается постоянной величиной: ω = const. Так как к моменту появления второй точки опоры угловая скорость больше нуля, то по­ворот элемента будет продолжаться (рис. 3, г) и произойдет потеря контакта с первой точкой опоры. На пороховой элемент действует сила тяжести G, вызывающая замедление вращения элемента, поэтому угловая скорость ω посте­пенно уменьшится до нуля, а затем начнется поворот элемента в противоположном направ­лении с угловой скоростью ω', и он вернется в исходное положение. В некоторых случаях, при больших значениях угла α, угловая ско­рость вращения может остаться положитель­ной в момент достижения значения α= 90°, в этом случае происходит переворот элемента.

Для проверки переворота элемента ис­пользуются следующие уравнения:

где ω₀ - значение угловой скорости ω в мо­мент отрыва от второй точки опоры;

ε_с - среднее значение углового ускорения при повороте элемента на угол α_π -α₀;

t_п - промежуток времени, за который про­изойдет поворот на угол α_π - α₀;

α_п = 90° - значение угла α в момент пере­ворота элемента;

α₀ - значение угла α в момент отрыва по­рохового элемента от второй точки опоры.

Решая полученные уравнения, можно найти значение угловой скорости ω в момент переворота элемента т. е. когда α = α_π = 90°. При этом если угловая скорость ω окажется больше нуля ω>0, то произойдет переворот порохового элемента; если ω ≤ 0, то перево­рот невозможен, т. е. пороховой элемент вер­нется в исходное положение с двумя точками опоры.

4. Изменение положения точки опоры (см. рис. 3, д).

Если точки соприкосновения располо­жены на одной стороне порохового элемента и линия действия силы тяжести не проходит между точками опоры, то поступательное и вращательное движения элемента продолжат­ся, но зона контакта с нижним слоем элементов переместится из точки C₁ в точку C₂. После по­явлении второй точки контакта под действием силы тяжести угол α будет продолжать увели­чиваться, в результате элементы в зоне точки C₁ перестанут соприкасаться. В этом случае единственной точкой контакта между порохо­выми элементами будет точка C₂, а для расчета движения элемента можно использовать систе­му уравнений (1), приняв в качестве начальных координат точки контакта координаты C₂.

Задача решалась в двумерной постанов­ке с помощью численного моделирования на ЭВМ. При этом рассматриваемое сечение раз­бивалось на отдельные элементы с размерами Δ х и Δу. Таким образом, был сформирован двумерный массив элементов (рис. 4).

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

После того как координата центра тяже­сти и ориентация порохового элемента опреде­лены, начинается расчет процесса его падения. При этом всем элементам массива, в которых в данный момент он располагается, присваива­ется значение, равное 2. Падение продолжает­ся до тех пор, пока под пороховым элементом остается свободное пространство.

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

При моделировании движения порохово­го элемента уравнения основной системы ре­шаются численно с помощью метода Эйлера. В данном случае используется регулярная вре­менная сетка. Она образуется линиями t_n = n% (n = 0, 1, 2, ...), где τ - шаг по времени. Значения рассчитываемых параметров (v_x, ω, х, α) определяются в целых точках по времени т . Таким образом, можно вычислить линейную скорость v_x и угловую скорость по­ворота ω в любой момент времени, а также значения координаты х и угла α:

При расчете по уравнениям для каждого момента времени определяются координаты центра тяжести порохового элемента и других точек в основной системе координат ОХ₀ Y₀:

где Х_i, Y_i - координаты точки с номером i в системе координат ОХ₀ Y₀.

Схема расчета представлена на рис. 4.

После определения координат точек 3, 4, 5 и 6 строится контур движущегося порохового элемента в момент времени t_n .

Если более 50 % площади какого-либо элемента массива находится внутри постро­енного контура 4-6-3-5, то считается, что

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

После определения положения падаю­щего порохового элемента проводится про­верка выполнения условий прекращения рас­чета. Если ни одно из условий не выполняется, расчет продолжается и определяется положе­ние порохового элемента в следующий мо­мент времени t_n+1. Если какое-либо условие выполняется, порядок расчета изменяется в соответствии с правилами, описанными выше. После полной остановки порохового элемента элементам массива, в которых в данный мо­мент он находится, присваивается значение 1 и процесс расчета перемещений порохового элемента заканчивается, генерируется новый пороховой элемент.

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

При расчете пористости в каждом из гори­зонтальных сечений определяется число пустых элементов массива (т. е. имеющих значение 0) и число заполненных элементов (имеющих значе­ние 1). Затем вычисляется значение пористости для данного горизонтального сечения:

где m_K - значение пористости для горизон­тального сечения, имеющего координату К;

i_K - число заполненных элементов массива для сечения с координатой K.

На основе составленной модели был про­веден ряд расчетов по формированию поро­хового заряда для пушки АК-230. На рис. 5 показано заполнение цилиндра пороховыми элементами. Рисунок соответствует одному из расчетов и изображает часть заснарядного про­странства (длиной 48 мм) в начале (рис. 5, а) и в конце (рис. 5, б) его заполнения.

В табл. 1 приведены результаты 20 расче­тов, выполненных при начальных данных, со­ответствующих конструктивным параметрам 30-милиметровой пушки АК-230.

Анализ полученных данных показывает, что значения пористости для различных горизон­тальных сечений цилиндра могут значительно отличаться друг от друга. Средняя пористость заряда также изменяется в некоторых преде­лах - она зависит от уровня заполнения гильзы пороховыми элементами (вдоль оси OY₀).

На рис. 6 представлены графики распре­деления пористости по длине каморы. Первый график соответствует расчету 3 (см. табл. 1) с резким изменением пористости в районе дна каморы, второй - расчету 8 с наиболее равно­мерным из всех расчетов распределением по­ристости по длине заряда.

При моделировании формирования за­ряда наблюдалась значительная разница в по­ристости для соседних сечений. Например, в расчете № 3 для сечений с координатами от 26 мм до 32 мм значение пористости изменя­ется от 0,2217 до 0,7322. Такое неравномерное распределение пористости может приводить к значительным перепадам давления пороховых газов для соседних сечений каморы при расче­те внутрибаллистических параметров.

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

Графики изменения давлений пороховых газов в районе дна боеприпаса - р_с (х), в райо­не дна каморы - р_к (х) , а также графики измене­ния скорости боеприпаса V показаны на рис. 7. Кроме того, приведены кривые изменения давления пороховых газов в случае постоян­ной пористости по длине заряда. Все кривые на графиках представляют собой зависимости давлений пороховых газов от координаты дна боеприпаса х. Нулевое значение по оси х со­ответствует дну каморы орудия.

На основании полученных данных уста­новлено, что характер распределения пори­стости по длине заряда значительно влияет на процесс горения порохового заряда и, соот­ветственно, на давление пороховых газов, на­блюдаемых в процессе выстрела. Для зарядов с одинаковой средней пористостью, но разным характером распределения пористости, макси­мальные давления пороховых газов могут от­личаться весьма значительно. Так, например, в расчетах 7 и 8 при близких значениях средней пористости разница между максимальными давлениями пороховых газов составляла око­ло 70 МПа.

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

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

Для расчета при условии постоянной по­ристости значения максимальных давлений по­роховых газов и скорости боеприпаса при вы­лете составили: P_{max м} = 377 8, МПа, P_см = 369 9, МПа, P_км = 353 3, МПа, V = 946 5602, м/с.

Экспериментальные значения макси­мальных давлений пороховых газов при вы­стреле для этого же орудия находятся в пре­делах от 307,034 до 362,012 МПа (данные получены при испытаниях артиллерийской системы на полигоне ФКП «НТИИМ») - один из экспериментальных графиков показан на рис. 7.

Таким образом, учет неравномерного распределения пороховых элементов по дли­не заряда при расчете внутрибаллистических параметров приводит к серьезным изменениям в результатах; подобная модель формирования заряда позволяет значительно лучше описать процесс горения реального порохового заряда, она может быть использована для вычисления начальных значений пористости в уравнениях газодинамического подхода.