Simulation of the formation of powder charge consisting of granulated powder elements

| S pa ce R es ea rc h an d R oc ke t S ci en ce | UDC 623.5 Khmelnikov E. A., Zavodova T. E., Smagin K. V., Dubinina S. F. Simulation of the formation of powder charge consisting of granulated powder elements The paper focuses on simulation of the formation of powder charge consisting of granulated powder elements. The purpose of the study is to develop a model which makes it possible to take into account the uneven distribution of powder elements along the length of the charge when calculating the intraballistic characteristics according to the gas-dynamic approach.

Determination of the propellant gas pressure, munition departure velocity and other intraballistic characteristics is of great importance in artillery armament design.
Today, the gas-dynamic approach is widely used for solving this problem. One of the welldeveloped models simulating the powder burning process, the motion of propellant gases and powder elements based on this approach is described in [1,2]. In the existing model, the distribution of powder elements along the charge length at the initial time is assumed uniform, which fails to comply with real processes observed when firing. That is why we developed a model to form a powder charge consisting of granulated powder elements. Such a model will allow to take into account non-uniform distribution of powder elements along the charge length.
A system of equations applied for the gas-dynamic approach consists of the equation of continuity, equation of motion, equation of energy, equation of powder charge burning, as well as additional equations required to calculate interphase interaction forces and other parameters. Therefore, the system of equations is represented as follows [1][2][3][4]: where ρ j -gas density for charge part with index j ; m -mixture porosity; S -variable area of cross-section of chamber and bore; t -time; v w , -gas and solid phase motion velocities in bore, respectively; х -coordinate; ρ ρ ρ = + 1 2 -total density of combustion products mixture; τ τ w w 1 2 , -interaction force per unit of volume caused by the difference of velocities between phases for the first and second parts of charge; | Space Research and Rocket Science | G G 1 2 , -gas supply from powder surface in unit of volume per unit of time; П с -barrel perimeter; q c -heat flow directed to bore surface; τ c -force of gas friction against barrel surface per unit of area; ε -internal energy of propellant gas unit of mass; -heating value (potential) of powder; р -pressure; a -covolume; , -gas specific heat at constant pressure and volume; a j -total particulate count of granulated powder elements; δ -powder material density; ψ j -relative fraction of burnt powder for the first and second parts of charge; Λ 0 0 , S -initial volume and surface of powder grain; σ ψ ( ) -current burning surface to initial burning surface ratio; u k -linear velocity of powder burning; ae p -powder element shape factor; e 1 -half of initial web thickness. Index j in the equations designates the number of charge part for combined charges.
Initial and boundary conditions are represented as follows [1,2]. 1. Initial conditions ( where р 0 -shot-start pressure; ω j -weight of the j-th half charge; W j -chamber volume occupied by the j-th half charge; m 0 -initial porosity of charge; W км -gun chamber volume; u j 1 -powder burning rate at atmospheric pressure (unit rate of burning); v j -exponent in equation for burning rate; f j -force of powder. 2. Boundary conditions.
where х сн -coordinate of projectile bottom at that time; v сн -projectile velocity at that time.
To calculate intraballistic parameters, the equations of gas-dynamic approach include such parameter as porositym. Porosity is the void volume in the unit of volume occupied by powder elements. Thus, powder elements will occupy volume 1− m . As powder burns out, porosity will increase.
In most computation methods that involve the gas-dynamic approach, the initial porosity value is determined based on confinement data: where Δ -confinement; δ -powder density. In these equations for calculating the initial porosity value, it is assumed that this value is constant for the whole volume of the chamber, i.e. powder elements are uniformly distributed over the whole volume behind the shot. Thus, the conditions of powder burning and movement of propellant gases at the beginning of computation will be the same for the entire length of the chamber. However, the assumption regarding uniform distribution of the powder elements over the volume behind the shot does not correspond to real processes observed during firing.
Therefore, if porosity at the initial time is considered variable depending on the position of the volume element to be observed inside the chamber, we can build a new model with non-uniform distribution of powder elements. The proposed model of charge formation is developed with account for the following factors: • when filling the shell with granulated powder, each powder element takes a random spatial position and orientation; • after a powder element falls down on the lower powder layer, it can continue to move along the side surfaces of the formed volume of powder (i.e. "falling" of charge).
When simulating the powder charge formation, the volume behind the shot being filled with powder is represented as a vertical hollow cylinder of the same volume. In case of 2D simulation, only one vertical cross-section of the cylinder passing through its centreline is considered. In this case, the ОХ Y 0 0 coordinate system is used (Fig. 1). The cylinder is filled with powder elements.
Each powder element appears on the upper section of the cylinder at a random point of ОХ 0 -axis -the point where the element appears corresponds to the position of its centre of gravity at the initial time. Then the powder element falls down until it reaches the cylinder bottom or other powder elements located underneath. Then a new powder element appears on the upper section of the cylinder. Such an element will be randomly oriented in space, i.e. its centreline may be located at any angle a relative to ОХ 0 -axis.
After a falling powder element reaches the layer of elements located below, it will stop or continue to move, depending on the position of the element's centre of gravity relative to the support point. When calculating the motion of a powder element, in addition to the basic coordinate system ОХ Y 0 0 , another coordinate system Оxyz connected with the powder element is built for each falling element (Fig. 2, z-axis not shown). Point О , i.e. the element's centre of gravity, is assumed the origin of coordinates.
After contacting other elements, the element under consideration is exposed to the gravity, support reaction, and friction force. For calculations, we use the equations of forward motion in the projections on х -axis and y -axis, as well as the equation of rotation around z-axis, plus additional equations to determine the said forces: where q -powder element weight; a x , a y -projections of element acceleration on х -axis and y-axis; a -angle of turning of powder element relative to ОХ 0 -axis; I z -moment of inertia of powder element rela tive to z-axis; ε -angular acceleration; g -gravity acceleration; f -friction coefficient.
where x and y -coordinates of support point in Оxyz coordinate system; t -time. By solving the obtained equations, we can determine linear and angular accelerations of the element at any time, as well as its motion speed along х-axis, angular speed of rotation around z-axis, and coordinates of the support point in the coordinate system associated with the element. Since the support point coordinates in the basic coordinate system are known, we can calculate the coordinates of the element's centre of gravity and (with the a angle value known) completely determine the position of a falling powder element relative to other elements at any time.
Solving the equations allows to determine the time when the loads applied to the elements change; therefore, the equations no longer apply and computation is interrupted. Possible criteria for computation interruption: 1. The element starts a free fall (Fig. 3, a). While moving, a powder element can lose its contact with other powder elements; it is possible if the following expression is true where a l = or a d = depending on the spatial position of a falling element.
In this case, computation by equations (1) is interrupted, because a falling element is exposed to gravity force G only. Computation of further motion of the powder element is similar to computation of the element's free fall when it first appears on the upper section of the cylinder simulating volume behind the shot.
2. Element motion interruption (Fig. 3, b). In this case, the system of equations changes depending on the appearance of the second contact point: If acceleration projections are represented as second derivatives of the time coordinate, the following system of differential equations can be obtained: Ox, Oy -coordinate system axes; О -element's centre of gravity; G -gravity; N -support reaction force; F тр -friction force; h -distance between axis Oy and line of support reaction force; r -distance between axis Ox and line of friction force; a -powder element turning angle If contact points are located on two mutually perpendicular sides of the powder element, forward motion will be interrupted. Turning around z-axis may continue depending on the relation between values ∆r , h , and a. We can prove that the rotational motion of the element will be interrupted if the following expression is true: If the expression is true, computation of the element motion is interrupted, and a new powder element is generated. If the expression is false, computation of element rotation around z-axis will continue.
3. Interruption of rotational motion with forward motion in progress (Fig. 3, c).
If contact points are located on one side of the powder element and the gravitational vertical passes between the support points, forward motion of the element along x-axis will continue, while element rotation will be interrupted. To prove this assumption, it is sufficient to form two equations of moments relative to the first and second support points; solving them together, we can obtain value ε = 0. In this case, the following system of equations will apply:  If angular acceleration is equal to zero, the angular speed of turning around z-axis will remain constant: ω = const . As by the time of appearance of the second support point the angular speed is over zero, therefore, the element will continue rotating (Fig. 3, d) and the contact with the first support point will be lost. The powder element is exposed to gravity force G, that decelerates the element rotation; thus, angular speed ω will gradually drop down to zero, then the element will start to turn in the opposite direction at angular speed ′ ω until it reaches the initial position. Sometimes, at high values of angle a , the angular speed of rotation may remain positive at the time when it reaches α = 90º; in this case, the element will be turned over.
The following equations are used to check the element turnover: where ω 0 -value of angular speed ω at the time of separation from the second support point; ε с -average value of angular acceleration when the element turns through angle α α п − 0 ; t п -time interval needed for the element to turn through angle α α п − 0 ; α п =°90 -value of angle a when the element turns over; a 0 -value of angle a when the powder element is separated from the second support point.
Solving the obtained equations allows to determine the value of angular speed ω when the element turns over, i.e when α α = =°п 90 . If angular speed ω is greater than zero ( ω > 0), the powder element will turn over; if ω ≤ 0, turnover is not possible, i.e. the powder element will return to the initial position between both support points.
If contact points are located on one side of the powder element and the gravitational vertical does not pass between the support points, forward and rotational motion of the element will continue, but the zone of contact with the lower layer of elements will move from point С 1 to point С 2 . After the second contact point appears, angle a will increase under gravity; as a result, the elements in the zone of point С 1 will no longer contact with each other. In this case, point С 2 will be the only contact point between the powder elements and to calculate the element motion, we can use system of equations (1), taking the contact points of coordinate С 2 as the origin coordinates.
The problem was solved in a 2D environment using computer-aided numerical simulation. The cross-section under consideration was divided into individual elements with their sizes equal to ∆ х and ∆у. Thus, a 2D-element array was formed ( Fig. 4).
At the initial time, value 0 was assigned to all the elements of the array -this value corresponds to an empty unfilled element of volume. Then a powder elements appeared in an element of the array. Its position (angle a ) was randomly determined. In order to determine the coordinate of the centre of gravity and its position, the standard function for random number generation was used.
After coordinate of the centre of gravity and the orientation of the powder element are determined, computation of the element's falling process is initiated. The value equal to 2 is assigned to all the array elements where the powder element is located at that time. Falling continues until there is free space under the powder element.
This computation within the simulation model is implemented as a check of array values, which are located under the area occupied by the powder element. If at least one of the array elements located under the powder element has the value equal to 1, the motion of the powder ele ment will be interrupted and the process of the element's falling computation will be over. Value 1 corresponds to those array elements where the stationary powder element is located.
When simulating the motion of the powder element, the equations of the basic system are numerically calculated using the Euler method. In For computations based on the equations, coordinates of the powder element's centre of gravity and coordinates of other points in the basic coordinate system ОХ Y 0 0 are determined for any specific time: where Х i , Y i -coordinates of the point with number i in coordinate system ОХ Y 0 0 . The calculation diagram is shown in Fig. 4. After the coordinates of points 3, 4, 5, and 6 are determined, the contour of a moving powder element at time t n is plotted. Х 0 , Y 0 -basic coordinate system axes; x, y -local coordinate system axes; l ц -length of segment 1-2; γ -angle between ОY 0 -axis and segment 1-2; h -length of segment 1-7; r -length of segment 2-7; β -angle between segments 1-7 and 1-2; a -length of falling powder element along x-axis; a -falling element turning angle If over 50 % of the area of any array element is located within the plotted contour 4-6-3-5, it is assumed that this array element is occupied by a moving powder element, and the value equal to 2 is assigned to this element. Thus, a particular area consisting of array elements will be formed in accordance with the moving powder element. The smaller the sizes of array elements (i.e. the smaller the value ∆х and ∆у), the more precisely the plotted area will match the moving powder element contour.
After the falling powder element position is determined, a check for conformity with the criteria of computation interruption is carried out. If none of the criteria are met, the computation will continue and the position of the powder element at the next time point t n+1 will be determined. If any criterion is met, the computation procedure is changed in accordance with the rules described above. After a complete stop of the powder element, value 1 is assigned to those array elements where the powder element is located at that time, and the process of powder element motion computation will be finished while a new powder element will be generated.
After total weight of the powder elements in the cylinder reaches a value equal to the nominal weight of the powder element, the porosity value will be determined for all horizontal cross-sections of the cylinder, as well as the confinement for the whole volume of the cylinder.
During porosity computation, the number of empty array elements (i.e. the elements with value 0) and the number of filled elements (the elements with value 1) is determined in each horizontal cross-section. Then the porosity value is calculated for a given horizontal cross-section: where m K -porosity value for horizontal crosssection with coordinate К; i K -number of filled array elements for cross-section with coordinate K .
Based on the developed model, we conducted a series of computations to form a powder charge for the AK-230 gun. Fig. 5 shows filling of the cylinder with powder elements. The figure corresponds to one of the computations and shows a part of the volume behind the shot (48 mm in length) at the beginning (Fig. 5, a) and at the end (Fig. 5, b) of its filling. Table 1 contains the results of 20 computations completed with the input data corresponding to design parameters of the AK-230 30 mm gun. According to the analysis of obtained data, porosity values for various horizontal crosssections of the cylinder may substantially differ from one another. The average porosity of the charge also changes within certain limits -it depends on the level of filling of the shell with powder elements (along ОY 0 -axis). Fig. 6 shows graphs illustrating porosity distribution along the chamber length. The first graph corresponds to computation 3 (see Table 1) with a drastic change of porosity in the chamber bottom area; the second graph corresponds to computation 8 with the most uniform estimated porosity distribution along the chamber length.
During simulation of powder formation, we observed a considerable difference in porosity for adjacent cross-sections. For example, in computation No. 3 for cross-sections with coordinates In order to determine how the porosity distribution along the charge length affects propellant gas pressure values when firing, computations of intraballistic parameters based on the gas-dynamic approach were conducted for all formed charges. Fig. 7 shows the graphs of propellant gas pressure variations in the munition bottom areap x с ( ), in chamber bottom areap x к ( ), as well as munition velocity variation graphs V. Moreover, the figure shows the curves of propellant gas pressure variation in case of constant porosity along the charge length. All the curves shown in the graphs are dependencies of propellant gas pressures on munition bottom coordinate x. Zero value on x -axis corresponds to the gun chamber bottom.
According to obtained data, we have found out that the behaviour of porosity distribution along the charge length has a considerable effect on the powder charge burning process and, therefore, on the pressure of propellant gases observed during the firing process. For charges with the similar average porosity, but with different behaviour of porosity distribution, maximum pressures of propellant gases may differ to a great extent. Thus, for example, in computations 7 and 8, at close values of average porosity, the difference between maximum pressures of propellant gases was about 70 MPa.
The observed phenomenon can be explained as follows: some charges have areas with a drastic porosity variation in a relatively short segment of the charge length. This leads to a local pressure rise (or drop) in these areas during the powder burning process.
The calculated maximum propellant gas pressure values are given in Table 2.
For calculation with constant porosity, the maximum propellant gas pressure values and munition velocity values during firing were as follows: The experimental values of maximum propellant gas pressures during firing for the same gun lie in the range of 307.034 to 362.012 MPa (the data is obtained when testing the artillery system at FSOE "NTIIM") -one of the experimental graphs is shown in Fig. 7.
Therefore, the calculation of intraballistic parameters with account for non-uniform distribution of the powder elements along the charge length leads to dramatic changes in the results; Note: P max м -maximum propellant gas pressure through the whole period of computation for the entire space behind the shot; P см -maximum propellant gas pressure through the whole period of computation near munition bottom; P км -maximum propellant gas pressure through the whole period of computation near chamber bottom; V -munition departure velocity.

| Space Research and Rocket Science |
such a model of charge formation allows to make a way better description of the real powder charge burning process. This model can be used to calculate initial values of porosity in gas-dynamic approach equations. Bibliography