ESTIMATE OF PROJECTILE INITIAL VELOCITY AS A SOLUTION OF A TWO-POINT BOUNDARY VALUE PROBLEM

The estimating of a projectile initial velocity is formulated as a two-point boundary value problem. To solve it, the data of a Doppler Radar or the results of solving the Cauchy problem can be used. The projectile initial velocity v0 estimation process is based on the numerical solution of a system of ordinary differential equations and the bisection method. The iterative calculating process is interrupted when a predetermined accuracy of a projectile initial velocity and a predetermined value of the width of velocity's search interval is reached. In the article, the block diagram of the algorithm for the projectile initial velocity process is developed. The Mathcad program code for mathematical modeling and a computer simulation of the projectile initial velocity estimation process for a 57mm armor-piercing projectile of ZIS-2 anti-tank gun 1943 model is given.


Introduction
A large number of practical problems of external ballistics of barreled artillery systems require information about a projectile initial velocity (PIV), briefly denoted as 0 v .This type of assessment is based on experimental data containing information about a projectile velocity 1 at trajectory points with known time or at trajectory points with known horizontal ranges.Typically, methods and devices that implement 0 v estimates use empirical mathematical dependencies based on assumptions about the nature of a projectile deceleration process (Alon, 1989), (Zaytsev et al, 2004) or general information about a law of frontal air resistance (Semenov, 2013).
The aim of this paper is to substantiate the PIV estimating process for tubed artillery based on the numerical solution of a two-point BVP.The PIV 2 estimating is based on: -multiple repeated solution of the system of ordinary differential equations (ODEs) of the longitudinal movement of an artillery projectile using a method of numerical integration; -comparison of the received result with the right boundary condition.Unlike previously known solutions, the technique assumes the use of the actual law of air drag, ballistic coefficient (BC), atmospheric parameters and actual angle of departure.If necessary, this data can be corrected before the start of computational or physical experiments.

An instantaneous projectile velocity (IPV) as a function of time
It is known that the longitudinal motion of an artillery projectile (stabilized by rotation) in the Earth's atmosphere can be described by the system of ODEs with an independent argument time (t) (Burlov et al, 2006, p.249).Systems of ODEs (1) are obtained by assuming the absence of the influence of lifting and lateral forces.There is no wind effect. 1 The article considers only the longitudinal motion of the artillery projectile.Projectile movements such as rotational, nutational and precessional movements are not taken into account.
2 Under PIV or v0 is meant a fictitious velocity for TOF t=0.Muzzle velocity is the actual velocity of the projectile that left a gun barrel.Information about PIV is used to calculate tabular (graphical) firing tables or for developing firing control equipment. here x -the abscissa of the trajectory, m; y -the ordinate of the trajectory, m; θ -angle of velocity vector relative to the horizon plane, degree or radian; v -IPV, m/s; t -TOF, s; g -acceleration of gravity, m/s 2 ; J -acceleration of the air resistance force, m/s 2 : where C -BC is estimated by the expression of Siacci (Germershausen, 1982, p.159), m 2 /kg; H(y) -a function characterizing the dependence of air density vs altitude; F(v) -a function characterizing the dependence of the air resistance vs an IPV.The ballistic coefficient is calculated as where i -value of the coefficient taking into account the shape of the projectile (dimensionless); from a formal point of view, this is a fitting coefficient; d -caliber of projectile (main diameter), m; q -projectile weight, kg.
On the basis of the foregoing, the Cauchy problem, namely the problem with initial conditions (IC), is formulated.The totality of its elements is shown on Figure 1.The acronym ADR means an Artillery Doppler Radar, which is intended for a projectile velocity measurement.
Since the system of ODEs (Figure 1, box 6) does not have "stiff" properties, it is therefore recommended to implement a numerical solution by dint of the fourth-order Runge-Kutta method or Bulirsch-Stoer method.As a result of the solution, there are five trajectory parameters: TOF (t), horizontal range (x), height of trajectory (y), angle of velocity vector (see ( 1)) (θ), IPV (v).
From a mathematical point of view, the difference between PIV and IPV is as follows: PIV is a constant; IVP is the interval from the PIV value to the so-called collision velocity (that is the velocity of the projectile, which has an impact on a target).
As an example, we can estimate the external ballistics of a 57mm 53-BR-271SP Armor-Piercing Tracer Solid projectile with a weight of 3.14 kg, fired from an anti-tank gun ZIS-2 M1943, whose gun barrel had a length of 73 caliber (Burlov et al, 2006, p.425).This artillery cannon was widely used in the World War II on the Soviet-German front.The serial production of the ZIS-2 was carried out in the period from 1941 to 1949.
The left graphical window (Figure 2) shows the dynamics of the deceleration process as a change of the IPV value vs TOF3 .The values of the variables t, x, y, θ, v were determined every 1.1 milliseconds of the projectile's flight, which corresponds to 1000 points inside the time interval.For each curve, the value of BC was constant for the whole domain of functions V=f(t).
The curves of Figure 2 are monotonically decreasing functions of TOF.In (Semenov, 2013) they are called "trajectories" of velocities.In the time interval (0-1.1 s), the graphs in form are close to linear functions.The dashed line t = 0.5 s (time after leaving the gun muzzle), crossing the velocity-time curves, forms three points of intersection.Their position is determined by the time interval t and the value of the BC.The higher the ballistic coefficient, the more significant the reduction of IPV in the chosen time interval by the constant angle of departure.Thus, an advanced knowledge of the BC value and time delay will allow to restore 0 v .
In the Cartesian coordinates Y-X in (the right graphic window of Figure 2), the shape of the trajectories is shown.Numerals 0.6, 0.5, 0.4 indicate a discrete decreasing of the BC value.
The growth of the BC causes a decrease in the horizontal range and a decrease in the value of the summit of a projectile flight path.The displacement of these parameters causes significant transformations of the descending branch of the trajectory in comparison with the ascending branch.
In order to simplify the calculations and reduce the computation time, Siacci air drag law was used in the form of the F-type curve (Germershausen, 1982, p.160 with the velocity range 0-1200 m/s replaced by a quadratic dependence, which has the form: where the value of the coefficient si B is 0.0002696 (dimensionless).
The maximum degree of the power by variable v by Siacci air drag law corresponds to 10, but formula (4) allows its reduction to a power of 2. The result of this approximation is shown in Figure 3.It should be noted that the comparison of the trajectories in the Y-X coordinates calculated according to Siacci air drag law and using formula (3) (only the law of the frontal air resistance was changed, the remaining parameters retained the previous values) leads to a decrease of a horizontal range about 10% regarding the result for Siacci.However, the replacement of (3) by ( 4) does not affect an accuracy of the solution of the BVP, since the conditions on the boundary are also calculated using (4).The gain from replacement consists of mathematical expressions for compactness and saving of calculation time, which is important for repeatedly recurring calculations at the testing stage and for checking the algorithm correctness.Time measurement for the numerical solution for the system of ODEs (1) using the fourth-order Runge-Kutta method fulfilled on the PC Lenovo Ideapad Z510 and Mathcad 15 software gave the following results: 17 ms for Siacci air drag law contra 14 ms for the quadratic law.The gain in calculating time was 3 ms, that is 17%.
If the BC is constant and PIV is a variable value, then the family of curves will have the form shown in Figure 4.The invariable ballistic coefficient (C = 0.5 m 2 /kg) and a discrete increasing of PIV ( 0 v = 980; 990; 1000 m/s) lead to the velocity change in the velocity-time coordinates in the form of three monotonically decreasing curves (the left graphic window of Figure 4).As before, the three curves characterize the trajectories of a projectile flight by the angle of departure of 19 angular minutes.Just like before, the presented graphs are close to linear dependences.However, the difference between Figure 4 and Figure 2 is the parallelism of the graphs of functions relative to each other, which predetermines the absence of intersection points.That is why the approximation of the functions (Figure 4) by the polynomial of the first degree will give the same tangents of the line slope (the condition of parallelism).
The right graphic window of Figure 4 shows the changes in the trajectories of the projectile in the Cartesian Y-X coordinates.The PIV increases discretely (v 0 = 980, 990, 1000 m/s) but the value of BC and the angle of the projectile departure remain constant.
If a velocity value at a given time is determined through the ADR measurement or by means of the Cauchy problem solution and if the angle of a projectile departure and the coordinates of the departure point are known, then in this case one can formulate a two-point BVP.

The two-point boundary value problem
The essence of the two-point BVP is that, in contrast to the Cauchy problem, the boundary conditions are not strictly defined on the left (as ICs), but also on the right boundary.For example, the value of velocity condition (IPV t=k ) restricts the result on the right boundary.Thus, the search for a solution to the systems of ODEs (1) that satisfies the boundary conditions for the time interval t=k (see Table 1) is organized as a procedure consisting of the following items: 1. Solution of the Cauchy problem described in the previous section; 2. Comparison of the result with a given condition on the right boundary; 3. Making a decision about changing the ICs in the ODEs with the subsequent return to point 1 or interrupting the calculations.The value of the projectile velocity at the right boundary of the interval is determined by a theoretical calculation or a physical experiment.In the first case, the Cauchy problem (1) is solved according to the beforehand known ICs.The estimation of 0 v by an experimental determination is carried out by using a ballistic chronograph or ADR.
The result of solving the boundary problem is obtaining the velocity of the projectile on the left border of the interval, herewith the fourth element of the matrix of ICs, v 0 (see Figure 1, box 5) is determined.
It should be noted that sometimes the boundary problem is called "bringing the projectile velocity to the muzzle".This idiomatic expression is valid only from a mathematical point of view.In view of the fact that the physics of the projectile motion process in a gun barrel and in the air by the shot is very complex, the muzzle velocity   m v and PIV   0 v are not equal to each other due to the presence of the aftereffect."PIV is a fictitious velocity that a projectile would have at the point of departure, so that in the absence of a period of aftereffect, its velocity outside the aftereffect line segment would be equal to the actual velocity.Such a velocity is 1-2% greater than the projectile muzzle velocity" (Burlov et al, 2006, p.115).
Thus, the fictitious value of 0 v exists as a phenomenon of a simplified mathematical description of a shot as a physical process, but plays an independent role in the preparation of tabular or graphical firing tables.Consequently, it is necessary to know it when the two-point boundary problem needs to be solved on the basis of experimental data.
If the condition on the right boundary is regulated by the IPV at the point t = 0.5 s, then the solution will be found in the iteration at which . (5) In this expression, the decrement   0 5 .

, v t v c
correspond the IPV at a point with a time t=0.5 s by PIV= 0 v .The velocity 0 v is the solution of the Cauchy problem.This is reminiscent of the subscript "с".The subtrahend is the user-defined right boundary condition.Expression ( 5) is a nonlinear equation of the argument 0 v .In fact, the right-hand side of ( 5) is an admissible error of the difference denoted as f  .Under these conditions, (5) is transformed to the form Since solution ( 6) is carried out in the interval of PIV 02 01 v v  , it is rational to impose an additional condition on its limiting width (after which further reduction of the interval is not rational).This restriction makes the iteration interrupt mechanism more flexible.
Nonlinear equation ( 6) is solved numerically by the method of bisection (division of a line segment in half).The algorithm of action is shown in Figure 5.
Box "1" is intended for entering the necessary constants.The block diagram assumes that the right boundary condition for velocity is determined from outside, for example, using the ADR or another method.Box "2" sets the limits of the search interval for initial speeds.It is assumed here that condition (5) has different signs on the left and right selected boundaries, and the left-hand side of function ( 6) is continuous.The task of box "3" is to divide the line segment of the search in half.The solution of the Cauchy problem is carried out in box "4".A comparison of the calculated value and the data of the ADR is performed in box "5".Boxes "6" and "7" serve to interrupt the iterative process of estimating the value of the PIV.The calculations are terminated upon reaching the predetermined accuracy f  of determining the initial speed and the pre-set value of the width of the search interval v  .In boxes "8-11", the process of controlling the change in the boundaries of the interval of the search for the PIV when it is divided in half is carried out.The number is 0.5 s present in the description of blocks "1.1" and "5" shows that the solution is sought for the right boundary condition of 0.5 s (the velocity is fixed at 0.5 s from the start of the flight).
The computer simulation of a BVP solution will be illustrated due to the Mathcad engineering math software and the sbval function (Solodov & Ochkov, 2004, p.101).Such commented program code is presented below.
At the first stage, it is necessary to determine the angle of departure The starting value z is entered to the program-code as a matrix that contains one element z 0 .
The matrix-column D is the right-hand parts of the system of ODEs (1).It includes variables with the following notation: -the abscissa of the trajectory  The fourth element of the matrix D stores information about the BC, about the function H(y) and about the law of an air drag.When using an air drag law other than quadratic, it is necessary to modify this element of the D-matrix.
The matrix-column L contains the ICs in the form of their known numerical values and unknown variables.In our case, these are 0 The parameter of the left side of expression t1 corresponds to the conditions on the left boundary of the interval.Here 0 z is specified as an element of the matrix z.
The matrix-column score stores the difference between the initial condition given at time-point t2 and the value of the sought solution at that point.The dimension of the matrix L is equal to the dimension of the matrix z.The parameter t2 denotes an employment of conditions from the right boundary The value of the velocity is 976.945m/s and it is the speed of the projectile at the right boundary of the time interval 0-0.1 s. and corresponds to the solution of the Cauchy problem with the following parameters and ICs: the parameters: C=0.5 m 2 /kg; law of air resistance -formula (3); change in air density -Eberling's formula; time interval of integration 0-1.1 s. with 1000 calculation points; ICs: x 0 = 0; y 0 = 0; θ 0 = 19 min.; Calling the solver function sbval for solution of the two-point BVP is performed by the command: S:= sbval (z, t1, t2, D, L, score).
The matrix-column of the solutions S has the same dimensions with the score and z matrices.As a result, we obtain S = 989.987m/s.The error in the numerical solution of the BVP was 0.0013% (the actual value of the PIV was 990 m/s).
If instead of the instantaneous velocity v in the ODEs, the velocity projections on the OX, OY axis are used, then the system of ODEs ( 1) is transformed into expression (8)4 (Chujko, 1976, p.38, 42) where u is the projection of IPV on the abscissa axis, m/s; ω is the projection of IPV on the ordinate axis, m/s.
Formulating the instantaneous velocity v through its projections on the axis OX, OY, we modify the first three equations ( 8) and obtain (9): Similar actions can be performed with respect to the quadratic law of air resistance (4):

  
. The obtained value of the IPV agrees with the value which was found in the previous example (the condition of expression ( 7)).The value of the relative error is 0.0001%.
Solving the BVP for the system of ODEs ( 9

0
(zero in subscript is the symbolical index) in radians (u=19 angular minutes), the ballistic coefficient and the acceleration of gravity as a will enter the starting value to find the required PIV (m/s), the starting and ending values of the time interval (s) from the virtual matrix Y previous example, we find the values of the parameters u,  in 0.1 s after a 57mm projectile has left the gun barrel.The following result is a solution of the Cauchy problem for the system of ODEs (9): u = 976.938m/s;  = 3.856 m/s.Since the horizontal and vertical projections of the velocity v in the Cartesian coordinate system are related to the velocity by the relation

.
) at u = 976.938m/s;  =3.856 m/s (right boundary conditions) we obtain the velocity projections on the OX, OY axis at time t = 0 (unknown IC) equal to u 0 = 989.975m/s; 0  = 4.903 m/s, from which it follows that the PIV 0 v is equal to:The errors in the numerical solution of the BVP for the system of ODEs (1) and the system of ODEs (9) are the same and equal to 0.0013%, so the uniqueness of the solution of the two-point boundary problem takes place.From the point of view of practice, it is of interest to determine the dependence of the error in estimating the PIV PIV  as a function of error in determining the boundary condition RB PIV   .The results of such an assessment are shown in Figure 6.