CC BY
119
http://dx.doi.org/10.5937/vojtehg66-15328
Introduction
As a physical phenomenon, a tubed artillery shot belongs to the class of fast processes. One of the ways of studying them is based on
CO
0
CO
I
SB
01 !± cp
MATHEMATICAL MODELING AND COMPUTER SIMULATION OF A BASIC PROBLEM OF TUBE ARTILLERY EXTERNAL BALLISTICS BY MEANS OF THE MATHCAD SOFTWARE
Vadim L. Khaikov j
independent researcher, Krasnodar, Russian Federation, e-mail: wadimhaikow@inbox.ru, ORCID iD: http://orcid.org/0000-0003-1433-3562
w
S3
TS
FIELD: Computational Ballistics I
ARTICLE TYPE: Original Scientific Paper 2
ARTICLE LANGUAGE: English
Summary:
The paper presents mathematical modeling and computer simulation (MM&CS) in the area of numerical solving of the basic problem of external ballistics for tubed artillery. Five-stage MM&CS scheme for conducting a ballistic simulation is developed. It is shown that a formal mathematical procedure allowing to solve the basic problem of external ballistics is a numerical solution of the Cauchy problem for a system of ballistic "s differential equations. The trajectory of a projectile flight for the 57-mm ^ ZIS-2 anti-tank cannon is estimated. A solving algorithm and the Mathcad program code are given. The numerical solution for a system of four first order ballistic differential equations is a five-dimensional space.The o possibility of visual presentation for a numerical solution was proposed in u the form of a square matrix. The boundaries of each subspace are determined. A procedure based on spline functions is developed for checking the correctness of the numerical solution. As a result of such verification, the effects of a light increase in the error at the edges of the integration interval are observed. A comparison of the numerical solution of the basic ballistics problem is conducted by means of "soft" and "stiff" solver-functions.The trajectory parameters estimated by "soft" and "stiff" methods are the same up to the fifth decimal place.
Key words: external ballistics, projectile trajectory, mathematical modeling, computer simulation, Cauchy task, numerical solution, Mathcad.
cp
I
CD CD
"o
>
03
o CM
of
UJ
a.
Z) O
o <
o
X
o
LU
H
^
OH <
H
<
CD >o
X LU H
O
O >
the method of mathematical modeling and computer simulation which replaces the real physical process by its formalized mathematical description. A problem is thus dealt with a significant amount of calculation, but it can be partially solved by using computer technology and automating calculations due to a program environment. As such, the Mathcad software was selected.
In view of the insignificant number of printed publications related to the technology of conducting external ballistic calculations in Mathcad1, as well as their non-systemic nature, the author wanted to contribute to filling this gap. To this goal, basic program modules have been developed for making it possible to carry out the MM&CS process for an external-ballistic trajectories assessment and to realize the calculations of flight path parameters.
In this part, the paper is a continuation of (Khaikov, 2018), where a method for estimating a projectile initial velocity as a two-point boundary value problem was proposed.
The aim of this paper is to solve the following subtasks:
- the development of a program module for the Cauchy problem solution due to various numerical methods for "soft" systems of ordinary differential equations (ODEs);
- the proposal of a method for the correctness verification of the obtained numerical solution;
- the application of numerical methods for the solution of the classical basic problem of external ballistics as probably "stiff" ODEs.
All the subtasks must be determined in the Mathcad software.
MM&CS in the area of external ballistics
MM&CS (mathematical modeling and computational experiment2) in the ballistics of tubed weapons is regarded here as a research methodology and technology, based on applied mathematics, mechanics, programming and computer technologies for:
- a deeper understanding of a shot as a physical phenomenon;
- developing new and improving existing designs of small firearms and artillery tubed weapons;
- developing new and improving existing ammunition;
- developing measuring instruments for experimental ballistics goals.
1 Publications from Eastern Europe were checked. The wide capabilities of the Mathcad software are used in many engineering activities. A preliminary study carried out by the author showed that there is no manual for a Mathcad application in analyzing internal and external ballistics of tubed weapons.
2 Here the term "computational experiment" is used as a synonym for "computer simulation".
With regard to ballistics problems, the process of MM&CS can be conditionally divided into five stages (Figure 1). Their names and essence are listed in Table 1.
Figure 1 - A generalized presentation of the MM&CS process Рис. 1 - Обобщенное представление процесса матмоделирования Слика 1 - Уопштена представа процеса MM&KS
Table 1 - Name of the stage and the essence of the performed actions Таблица 1 - Характеристика этапов математического моделирования Табела 1 - Називи фаза и кратки описи предузетих активности
№ Name of the stage The essence of the performed actions
1 2 3
1 Development of a mathematical ballistic model The construction of a mathematical model as a formal mathematical description of a physical phenomenon or object.
2 Development of a calculation algorithm Selection of numerical methods for the calculation and development of a computational algorithm.
3 Choice of the programming environment and the development of a program code Implementation of the developed algorithm in the selected software for engineering calculation. An important element in the development of program code is a definition of the existence area for the initial database (DB) and the development of management tools.
со о со
I
00 сч !± с
Í
ф
£
(Л
S3 I
ТВ .Q
Ф Ф
Ф .Q
3
su
_Q 2 ü
о
го
.Q ГО
t5 с о
:d
о.
тз %
О)
тз о
75 о
£ ■<3
CD CD
O >
O CM
DC LLJ
OC ZD
o
o <
o 2:
X
o
LU
I— >-
QC <
< -J
CD >0
X LU I— O z
—3
O >
№ Name of the stage The essence of the performed actions
1 2 3
4 Calculation performance Carrying out the calculations and the postprocessing of the resulting computing data.
5 Analysis of the results Analysis of the calculation results, comparison (if possible) with the data of physical experiments or with conclusions of other computational simulations.
In the fifth stage, researchers can go back to the previous stages and make the necessary changes (for example, they can change the program-code or refine the initial data). The ways to perform feedback are shown by a dashed line in Figure 1.
Further on, the paper will describe a sequence of stages for assessing unguided projectile external-ballistic trajectories.
Ballistic trajectories and their elements as a result of numerical integrating of a system of ODEs throughthe Mathcad solver-functions
The following will be demonstrated: 1) how to calculate a trajectory and evaluate its parameters on the basis of a numerical solution, 2) how to verify the correctness of the obtained numerical solution, and 3) how to solve the basic ballistics (sometimes the word "main" is used) problem of external ballistics by various numerical methods.
The realization of all five stages of the MM&CS process will be shown below by solving problem (1).
Numerical calculation of an external-ballistic flight path and its parameters
MM&CS Stage 1. Mathematical model of a projectile external-ballistic flight path.
The longitudinal motion of an artillery projectile 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). This type of mathematical expressions belongs to the type of the Point-Mass Trajectory Model:
— = v cos(O) dt
dy = v sin(o)
dt
dO g (a)
~dt =~ ~v COS(O) ddV = -( + g ■ sin(O))
(1)
CO
0
CO
I
SB
01
cp
(D £
in
here3 x - the abscissa of the trajectory (m); y - the ordinate of the trajectory (m); 0 - the angle of a velocity vector relative to the base of the $
trajectory (degrees or radians); v - the instantaneous projectile velocity |
2
(m/s); t - the time of flight (s); g - the acceleration of gravity (m/s2); ts J - the acceleration of the air resistance force (m/s2).
In order to assess an external-ballistic trajectory, we will perform an example for the flight of the 57-mm Armor Piercing Tracer Solid Projectile 53-BR-271SP (weight 3.14 kg) of the anti-tank gun, ZIS-2 1943 model4. The ZIS-2 gun is a so-called direct fire weapon and provides shooting with a flat trajectory.
MM&CS Stage 2. Projectile flight path calculation algorithm.
A calculation algorithm is a simple chain of performing computational Jj operations. It is presented in Figure 2. The chain consists of eight £ elements. In boxes 1, 2, and 6, the values of the constants are determined. Boxes 3 and 4 define individual functions, and box 7 sets a functional matrix-column. The value of a ballistic coefficient is calculated in box 5. -§
The numerical solution of the Cauchy problem is done in box 8. The difference between box 8 and boxes 1-7 is that box 8 is called the built-in Mathcad solver-function. Data from boxes 1-8 is necessary for the solver-function for its correct functioning. In the block diagram, Figure 2, it is
sn
_Q 2 O
'¡B
.Q
15 O
3 A more complete description of a system of ODEs and its elements is given in (Khaikov, 2018).
4 Despite the almost identical initial data, the difference between (Khaikov, 2018) and this paper is in the following: the direct external-ballistical task is solved here, and the inverse task fullfiled in the previous work. The direct task (basic ballistic problem) requires the 5 solution of the Cauchy task, and the inverse task deals with the solution of the boundary value problem.
>
implied that after box 8 the results (errors) are output in a tabular or graphical form.
CD CD
О >
О CM
DC Ш
DC ZD О
о <
о z
X
о ш
I—
>-
DC <
(Л <
-J
CD >o
X Ш I—
о z
о >
START
Ж
Determination of the boundaries Determination of the number of
of the integration interval —► points in the integration interval
N
Determination of the air drag law -Ь Determination of the law of
F(v) resistance H(y)
A
J 4 +
+
Calculations of the value of ь The determination of the initial
the ballistic coefficient conditions for a system of ODEs
The definition of the right-hand
sides of the system of ODEs
l7-
The solution of the Cauchy problem
l8-
EN □
Figure 2 - An algorithm for calculating a projectile flight path Рис. 2 - Алгоритм расчета траектории полета снаряда Слика 2 - Алгоритам за израчунаваше путане лета проjектила
MM&CS Stage 3. Program-code for a numerical solution.
In accordance with Figure 1, a process for program-code preparing is divided into two parts: the preparation of a database and the development of computational Mathcad procedures. The structure of the initial DB is given in Table 2.
Table 2 - Characteristics of the initial database Таблица 2 - Характеристика исходных данных Табела 2 - Карактеристике почетне базе података
I. Departure point conditions
1 X coordinate of the departure point (X, m) 0
2 Y coordinate of the departure point (V, m) 0
3 Angle of departure (Q0, angle minute) 10
4 Projectile initial velocity (Vo, m/s) 990
II. Integration interval
5 Time interval (t, s) tbeg.. tend
5.1 Beginning of the time-interval (tbeg, s) 0
5.2 End of the time-interval (tend, s) 1.1
6 Number of points (dimensionless) 1000
7 Integration step (At, milliseconds) 1.1
III. Functions, variables (constants) of the right parts of a system of ODEs
8 Air drag law (see (Khaikov, 2018)) quadratic
9 The law of density variation with altitude (H(y)) Eberling5
10 Coefficient that depends on a shape of a projectile, i (dimensionless) 0.483
11 Caliber of a projectile (main diameter) (d, m) 0.057
12 Projectile weight (g,kg) 3.14
13 Air density at the departure point (Pn, kg/m3) 1.206
14 Acceleration of gravity at the departure point (g, m/s2) 9.81
When using system (1), 12 constants (variables) and 2 functions are necessary for a numerical solution of the Cauchy problem.
The Mathcad-code program is commented below. Using formula 1 and Table 1, one starts programming the elements of the Cauchy problem.
Let us determine the caliber for an artillery projectile (57 mm = 0.057 m), the weight (3.14 kg) and the value of the i coefficient (according to the quadratic drag law):
CO
0
CO
I
00
01 !± CP
S
Ф £
(Л
S3 I
TS .Q
"S
ф ф
ф -Q
■3
SU
-Q 2 ü
о
то
-Q ТО
Ъ
С
0
1
3
CP
тз
s
О) С
TD тз о
ТВ о
d:=0.057 q:=3.14 i:=4.83
(2)
By Eberling law, H(y) is determined as: H(y)=10
_л /~>-0.000046(y)
CD CD
O >
O CM
DC LLJ
DC ZD
o
o <
o 2:
X
o
LU
I— >-
Q1 <
< -J
CD >0
x
LU I— O z
—3
O >
An angle of departure (in radians) is calculated as a set of angular degrees, minutes and seconds according to the formula:
Gradus:=0
Min:=10
Sec:=0
t Min Sec 60 :=-1 Gradus h---h -
180 1
60 3600,
= 2.909 • 10"
(3)
At a point of departure, the value of the acceleration coefficient of gravity is determined as 9.18 m/s2. Further on, it is necessary to determine the time interval of integration, its boundaries and the total number of points:
tbeg=0 tend=1.1 npoints-' = 1000 (4)
As a law of air resistance, we will use the quadratic dependence of the type:
F,ap (v )= Bs
(5)
where Bsi = 0.00002694. This simplification is made in order to avoid the complicated mathematical formulas that describe the law of air resistance.
The value of the ballistic coefficient is estimated from formula (3) (Germershausen, 1982, p.159). It is equal to 0.5 m2/kg. The acceleration
of the air resistance force J will be characterized by form6
J = 0.5 • ))• (0.00002694v2)
(6)
The initial conditions (for (1)) are determined as a matrix-column y, which will contain their known numerical values:
y :=
( y0 = x(0)^ y 1 = y (0)
y 2 = 0(0)
y3 = v(0).
3
2
v
6 Here it should be noted that the formula for calculating the ballistic coefficient and the method of its introduction into the J-equation for different countries (their national standards) may be not the same. In order to avoid ambiguous understanding, it must be agreed in advance.
In view of the fact that the initial velocity of the 57 mm armor-piercing solid tracer projectile 53-BR-271SP is 990 m/s (Burlov et al, 2006, p.425), the matrix-column y will look like:
CO
o
CO
I
SB
CM !± Cp
f
y :=
0 0
A
2.909 • 10" 990
(7)
V y
The matrix-column D is the right-hand part of the system of ODEs (1). It includes variables with the following notation: the abscissa of the trajectory x-y0; the ordinate of the trajectory y-yi, the angle of the inclination of the tangent 0 -y2, and the instantaneous projectile velocity
v -y3:
y3 • C0S(y2 ) y3 • sin(y2 )
D(t, y) :=
46
0.5 •
10
1000000
(l
---COS(y2 )
y3
A
2694
10000000
y32 1 + g sin(y2)
(8)
Calling the solver-function rkfixed for numerical solution, (1) is performed by the command:
<D £
in
S3 I
TS .Q
TB
<D <D
<D .Q
3
sn
_Q £
O
'¡B
TO .Q CD
t3 c
0
1
:d
cp
Result1:= rkfixediy, tbeg, tend, nPomt, D).
For solving differential equations, the solver-function rkfixed implement the fourth order Runge-Kutta methods with a fixed step.
The variable Resultl is the matrix containing the results of the numerical solution of the system of ODEs (1). In this case, the matrix has a size of 5 by 1001 elements and contains 5005 numbers. Five rows of the Resultl matrix are the independent argument time (t) and the elements of the matrix y (or D). The first column (argument time) counts 1001 numbers. This number includes the initial conditions plus the numbers of the points in the integration interval.
<D T3 O
15 o
iS
3
CD CD
O >
O CM
DC LLJ
DC ZD
o
o <
o 2:
X
o
LU
I— >-
Q1 <
A special feature of the basic problem of external ballistics is that the system of differential equations (1) contains the quantity J, which does not have an exact analytic expression. The fact is that J is a empirical function which can be described as
J = CH (y )F (v), (9)
where C - a ballistic coefficient estimated by the Siacci mathematical expression (Germershausen, 1982, p.159), m2/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 the instantaneous projectile velocity.
The functions H(y) and F(v) can be described by means of:
- an approximate analytical function;
- a piecewise function;
- a spline function.
< -J
CD >0
x
LU I— O z
—3
O >
For example, the 1943 air drag law7 can be described by a piecewise function (Konovalov, 1979, p.84):
(m ) =
0.157 0.033M + 0.133 3.9M2 - 6.419M + 2.8025831 1.5M -1.176
0.1 < M < 0.73 0.73 < M < 0.82 0.82 < M < 0.91 0.91 < M < 1.00
A
- 1.6M2 + 3.7632M - 1.828716 1.00 < M < 1.18
0.384sin(1.85M-1 ) 0.29M-1 + 0.172 - 0.011M + 0.301
1.18 < M < 1.62 1.62 < M < 3.06 3.06 < M < 3.53
J
which will lead to the complication of the program-code. The drag coefficient Cx43 is a function of the Mach number (M).
7 It means the law of air resistance that was developed in the USSR during World War II. It has a short name: the 1943 law.
By using the drag functions Cx(M) the quantity J will be calculated as J = 4.732 • 10"4 CH (y)alM 2Cx (M), (10)
where a0 - the speed of sound constant (about 340 m/s).
MM&CS Stage 4. Calculation.
An example of the numerical solution of ODEs (1) as a Resultl matrix is shown in Figure 4. Note that the Mathcad software independently enumerates the columns and rows. The vertical column and the horizontal row with numbers and gray background should not be confused with numerical data.
The first row contains the initial conditions for (1). The leftmost column is the independent argument (in this case, the time of flight). The values of the first column are limited by the boundaries of the integration interval.
Resultl =
0 1 2 3 4
0 0 0 0 2.909-10-3 990
1 1.110"3 1.089 3.162'10-3 2.Э9Э-10-3 599.955
2 2.2 10-3 2.178 6.311 '10-3 2.997'10"3 999.709
3 3.3-10"3 3.266 9.44a-10"3 2.976-10-3 999.564
4 4.4-10-3 4.355 0.013 2.965-10-3 999.419
5 5.5'10-3 5.443 0.016 2.954-10-3 999.272
Figure 3 - An example of data for a numerical solution of ODEs (1) Рис. 3 - Пример матрицы численного решения системы ОДУ (1) Слика 3 - Пример података за нумеричко решете ОДJ
со
0 со
I
00
01 !± CP
S
ф
£
(Л
S3 I
тз
.Q
ТВ
ф ф
ф .Q
3
_Ф
.Q £
О
ГО .Q ГО
Ъ
С
0
1
3
The time step of the Table is constant and it is determined as (in milliseconds)
Step =
^end tbegin npo rnts + 1
1.1 - 0 1001
= 1.09ms.
In order to use the found values independently from the matrix Resultl will create 6 variables; their purpose are explained in Table 3. If necessary, the variables can be combined into a new matrix.
CP
T3
s
O) С
TD тз о
TS о
T: = Result^00 0Ra: = Result1(3)
X: = Result1(1) Oct: = 180m1 Or*,
Y: = Result1(2) V: = Result1(4)
i 'Ш
CD CD
О >
О CM
ОС ш
ос
ZD О
О _|
< с;
X О ш
I— >-
ОС <
Usually, the output of tabular data is accompanied by graphs. Their coordinates type and scale must be convenient for further analyses.
Table 3 - Characteristics of variables of a numerical solution of the system of ODEs (1) Таблица 3 - Характеристика переменных численного решение системы ОДУ (1) Табела 3 - Карактеристике промен^ивих нумеричког решена система ОДJ(1)
№/№ Designation Specification
1 T A time of flight (t, s)
2 X A horizontal range (X, m)
3 Y A height of a trajectory (Y, m)
4 вка An angle of the velocity vector relative to the base of a trajectory in radians
5 вог An angle of the velocity vector relative to the base of a trajectory in degrees
6 V An instantaneous projectile velocity (m/s)
(Л <
-j
О ■О
X ш I— О
О >
MM&CS Stage 5. Analysis of results.
It is known that the analysis of any two quantities can be expressed analytically, in a tabular form and graphically. The numerical solution of the system of ODEs (1) enables ananalysis of tabular and graphic data. The graphs of functions performed in the Cartesian coordinates for the system of ODEs (1) are shown in Figure 4.
Figure 4 shows the complete set of dependencies for the five quantities which characterize the system of ODEs (1). This image is constructed in the form of a square matrix. The quantities indicated in the upper horizontal line (X, Y,6Gr, V, T) of Figure 4 are plotted respectively on the ordinate axis of each graph. The quantities showed in the left vertical column of this figure are plotted respectively on the abscissa axis. The total number of dependencies is 25. Five graphs along the matrix diagonal are functions of the magnitude itself. Thus, there are 20 dependencies. From this number, 10 are direct dependencies and 10 are reverse ones. The functions Y-X, V-X, and V-T are the most widespread ones in practice.
Sometimes for the value of Y as a height of the trajectory and as an argument of function, there are 2 points (ordinates) on the graphs (for one Y-value). Such ambiguity characterizes the ascending and descending parts of the trajectory.
A starting point in the form of a square is added to each graph. This is done to explain the dynamics of the projectile flight and it displays the initial state from the point of view of the physical process.
The domain of the numerical solution for the system of ODEs (1) is the five-dimensional space (5D) of X Y, 6Gr, V and T quantities. The boundaries of each of the quantities are shown in Table 4.
If the quantities X Y, 6Gr and V are the parametric equations from the argument T, then, in the three-dimensional space, they can be expressed by two types of 3D-graphs: X Y, 6Gr X, Y and V. A parametric equation and the three-dimensional coordinate space give an opportunity to create a better graphic image.
CO
0
CO
I
00
01 !± cp
Ф £ 45
(Л
S3 I
TS .Q
Ф Ф
Ф .Q
3
_ф
_Q 2 О
ГО .Q ГО
t3 с
0
1
3
cp
Ф T3
о
15 о
Figure 4 - Complete set of 2D dependencies between the values of the system of ODEs
(1)
Рис. 4 - Полное множество 2D зависимостей между величинами системы ОДУ (1) Слика 4 - Комплет 2D зависности измену вредности система ОДJ (1)
!
CD CD
О >
О CM
ОС ш
ос
ZD О
О _|
< с;
X О ш
н >-
ОС <
(Л <
-J
О >о
X ш н о
о >
Table 4 - Characterization of the boundaries of the quantities (X, Y, 6Gr, V, T) Таблица 4 - Характеристика границ величин (X, Y, 6Gr,V, T) Табела 4 - Карактеризаци^а граница величина X, Y, QGr, V, T
№ Physical quantities Min Max
1 The horizontal range (X, m) 0 1000
2 The height of the trajectory (Y, m) 0 1.1
3 The angle of the velocity vector relative to the horizontal plane, in degrees -0.5 2.909-10-3
4 The instantaneous projectile velocity (m/s) 820 990
5 The time of flight (t, s) 0 1.1
Verification of the correctness for external-ballistics trajectory parameters obtained as result of a numerical solution
In view of the fact that the existing solver-functions are nothing more than "black boxes", and the information about their implementation algorithm is characterized by incompleteness, it may be necessary to verify the correctness of the obtained solution by the substitution method.
At the first level of this process, we create a matrixline of discrete time values:
tin ' tbeg ' tbeg +
tend tbeg 1001
.. t.
end
At the second level, we perform the spline interpolation of the four obtained quantities (X, Y, 6Gr, V) that are the result of the numerical solution of the system of ODEs, and express them as a spline function with the argument tin. Thus, we define the four new functions: Xin(tin), Yin(tin), 0,n(tm) and Vin(tin).
X;
in(tin) Yin(tin ) ein(tin) Vin(tin )
= interp (lspline (T, X), T, X, tin) = interp (lspline (T, Y), T, Y, tin) = interp (lspline (T, e), T, e, tin) = interp (lspline (T, V), T, V, tin)
In this case, a linear spline function is used, but M allows the use of quadratic and cubic spline functions.
For the first differential equation (1), namely
dx(t)
dt
= v(t )cos(6>(t))
express the left part and right one in the form of the spline functions. For the left-hand side, we perform numerical differentiation and obtain the function LX(tin)
lXU:= dyW^n))
The right-hand side of RX(tin) is obtained by multiplying
RxOrn) := Vin(v)-cos (SinOin)) .
One visualizes the results (both the right and the left part) in the form of a graph in the Cartesian coordinates
Figure 5 - Comparison of the right and left part of the first differential equation of (1) Рис. 5 - Сравнение правой и левой частей первого уравнения системы ОДУ (1)
Слика 5 - Поре^еше левог и десног дела прве диференцщалне ]едначине (1)
The graph of the right and left part of the first differential equation of system (1) is presented as a decreasing function of time. Figure 5 shows its form for the integration interval 0-1.1 s. In order to find a magnitude of the relative error of their difference, one will look for (in percent):
1001 Lx(tin) - %У|
Delta
X1
(tin) -
CO
0
CO
I
00
01 !± cp
Ф £
(Л
S3 I
TS .Q
Ф Ф
Ф .Q
3
_ф
_Q 2 О
ГО .Q ГО
Ъ
С
0
1
3
cp
Ф T3
о
ТВ о
Lx'
(tin)
I
To improve the image, the value of DeitaX(tin) is summed with the constant a = 1-1Q'10.
CD CD
О >
О CM
DC Ш
DC ZD О
о <
о
Z X О ш
i>-
DC <
(Л <
-J
О >о
X ш н
о
О >
The graph in Figure 6 shows that at the ends of the integration interval there are "edge" effects and the error of numerical solving increases.
Figure 6 - Visual presentation of the relative error of the right and left parts of the first differential equation expressed in percent (integration interval 0-1.1 s.)
Рис. 6 - Визуальное изображение относительной погрешности правой и левой частей первого уравнения, выраженной в процентах (интервал интегрирования
0-1.1 с.)
Слика 6 - Визуелни приказ релативне грешке десног и левог дела прве диференцц'алне ]едначине изражене у процентима (интеграциони интервал 0 - 1,1s)
In a similar way, the left and right parts for the second, third and fourth differential equations of system (1) are compared (Figure 7).
There is a coincidence of the function graphs. The right-hand sides of the second and third equations of the system of ODEs (1) are monotonically decreasing time functions. Their shape is close to linear dependence. In contrast, the right-hand side of the fourth equation is an increasing function. By the adopted initial conditions, their domain is negative.
A visual presentation of the errors of the right-hand and left-hand sides of equations 2-4 of the system of ODEs (1) is shown in Figure 7.
A visual presentation of the relative error of the right and left parts for the second, third and fourth differential equation expressed in percent (integration interval 0-1.1 s.) is shown in Figure 8.
Figure 7 - A proof of the numerical solution for the second, third and fourth differential
equation
Рис. 7 - Проверка численного решения второго, третьего и четвертого дифференциальных уравнений Слика 7 - Доказ нумеричког решена друге, теПе и четврте диференци^алне
jедначине
со
0 со
I
00
01 !± Cp
ф
£
45
(Л
S3 I
ТВ .Q
ТВ
Figure 8 - A visual representation of the relative error of the right and left parts of
differential equations
Рис. 8 - Визуальное представление относительной погрешности правой и левой частей дифференциальных уравнений Слика 8 - Визуелни приказ релативне грешке десног и левог дела диференцц'алних jедначина
ф ф
ф
.Q
3
_Ф
.Q £
О
ГО .Q ГО
Ъ
С
0
1
cp
The conclusion for paragraph 2.2 is that the verification of the numerical solution can be carried out using spline functions with the help of which the approximation of the numerical solution is performed and the right parts of the systems of differential equations are calculated. The left-hand parts of the system of ODEs are calculated by numerical differentiation. At the final stage, we compared the right and left parts of the chosen interval of integration. For system (1), the error in the numerical solution does not exceed 110-7%.
ф
T3
о
TS о
i
CD CD
О >
О CM
DC LLJ
DC ZD О
о <
с; z
X
о ш
i>-
DC <
Numerical solving of the basic externai-baiiistics problem by numerical methods for "soft" and "stiff" systems of ODEs
From a practical point of view, it is interesting to consider the following problems:
I. Comparison of the calculation results obtained through various numerical methods.
II. Determination of the "stiff" properties of the system of ODEs characterizing the external ballistician problems.
Regarding the first problem: doing research for comparing the usefulness of various numerical methods by solving a system of ballistic differential equations will possibly allow an increase in a number of methods that can be conducted for numerical solving procedures. During the research, the level of errors, computational costs and the average solution time will be established.
The Mathcad solver-functions for solving the ODEs for "soft" and "stiff" systems are collected in Table 5. There are 7 solver-functions.
(Л <
-J
CD >o
X Ш I—
о
N
о >
Table 5 - Functions-solvers of ODEs for "soft" and "stiff" systems Таблица 5 - Функции-решатели ОДУ для «мягких» и «жестких» систем Табела 5 - Функци^е за решаваше ОДJ за „меке" и „круте" системе
Functions-solvers of "soft" systems of ODEs Functions-solvers of "stiff" systems of ODEs
Adams
rkfixed Radau
Rkadapt Stiffr
Bulstoer Stiffb
The commands for calling a solver-function and the required variables are shown below:
Result1:= Adams(y, tbeg, tend, n^ D, 10r9) Result2:= rkfixed(y, tbeg, tend, npon, D) Result3:= Rkadapt(y, tbeg, tend, ripoint, D) Result4:= Bulstoer(y, tbeg, tend, D)
The computational experiments conducted by the author showed that the calculations of system (1) performed by "soft" ODEs such as rkfixed, Rkadapt and Bulstoer give identical results up to the 5th digit after the decimal point (there were no experiments with the solver-function Adams).
Regarding the second problem: a parallel solution of the same ballistic problem by several numerical methods will make it possible to detect the appearance of probable "stiff" properties of the system of ODEs.
System (1) was solved with numerical methods for solving "stiff" ODEs: for instance, the Stiffb solver-function which implements the Bulirsch-Stoer method and the Stiffr solver-function - the Rosenbrock method.
To evaluate a trajectory using "stiff" methods, it is necessary to know the matrix Jkthat has a size n*(n+1): n rows and n+1 columns. The first column Jk contains the partial derivatives dF/dt, while the remaining columns and rows represent the Jacobi matrix dF/dy.
The structure of the matrix Jk applied to (1) has the following form
J (t, y ) =
dF0 dF0 dF0 dF0 dF0
dt dy1 dy3
dF1 dF1 dF1 dF1 dF1
dt ^0 dy1 dy3
dF2 dF2 dF2 dF2 dF2
dt dy(> dy1 dy3
dF3 dF3 dF3 dF3 dF3
dt dy0 dy1 ^2 dy3
For D(t,y) (8), the matrix Jk(t,y) of size 4*5 (4 rows and 5 columns) has the kind
cos(y 2 )
sin(y 2 )
gcos(y 2)
J ( y ) =
0 0 0 - y3 sin(y2)
0 0 0 - y3 COs(y2 )
0 0 0 g sin(y2 )
y3
0 0 0 g COs(y 2 )
y
2
0.0000538Cy3
100.000046 yj
(11)
The element of the matrix a has the form
CO
0
CO
I
SB
01 !± CP
i
(D £ 45 in
S3 I
TS .Q
TS
(D (D
(D .Q
3
sn
.Q 2 o
'¡S
.Q CD
t3 c
0
1
:d
CP
T3 £ O) C
TD T3
o
TS o
a =
dF3 = 1.237 • 10-9 ln(l0)Cy
10'
|0.000046y1
(12)
i TS
CD CD
0
01
(Л
>
The result of testing solutions (1) by the methods of solving a "stiff" system of ODEs Stiffb, Stiffr gives a similar result as "soft" methods. For such solutions, there are 5 digits after the decimal point.
£ Conclusion
In this article, it is proposed to consider that a numerical solution of the basic problem of external ballistics can be considered as a five-stage w process of MM&CS. ^ It is shown that a formal mathematical procedure allowing to solve
o the basic problem of external ballistics is a numerical solution of the Cauchy problem for a system of ballistic differential equations.
0 The flight trajectory of a projectile for the 57-mm ZIS-2 anti-tank
1 cannon is estimated and visualized. A solving algorithm and the Matcad
o
uj program-code are given.
The numerical solution for a system of four first order ballistic differential equations is a five-dimensional space. The possibility of a 2D visual representation for a numerical solution was proposed in a form of a square matrix. The boundaries of each subspace are determined.
A procedure based on spline functions is developed for checking the correctness of the numerical solution. For verification purposes, spline fj functions of three kinds can be used: a linear spline, a quadratic spline ^ and a cubic spline. As a result of such verification, the effects of a light >o increase in the error at the edges of the integration interval are observed.
A comparison of the numerical solution of the basic ballistics TE problem is conducted by means of "soft" and "stiff" solver-functions. The o article considers a possibility of solving a ballistic problem using five ¡3 types of the Mathcad solver-functions. The trajectory parameters estimated by "soft" and "stiff" methods are the same up to the fifth decimal place.
References
Burlov, V.V. et al. 2006. Ballistics of Tubed Artillery Systems. Moscow: Mashinostroenie (in Russian). (In the original: Бурлов, В.В. и др. 2006. Баллистика ствольных систем. Москва: Машиностроение).
Germershausen, R. 1982. Handbook on weaponry. Düsseldorf: Rheinmetall GmbH.
Khaikov, V.L. 2018. Estimate of projectile initial velocity as a solution of a two-point boundary value problem. Vojnotehnicki glasnik/Military Technical Courier, 66(1), pp.9-27. Available at: http://dx.doi.org/10.5937/vojtehg66-15097.
Konovalov, A.A., Nikolayev, Yu.V. 1979. Vneshnyaya ballistika. Moscow: 8
00 Ol !± CP
Tsentral'nyy Nauchno-issledovatel'skiy Institut Informatsii (in Russian). (In the original: Коновалов, А.А., Николаев, Ю.В. 1979. Внешняя баллистика. Москва: Центральный научно-исследовательский институт информации).
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ И ВЫЧИСЛИТЕЛЬНЫЙ ЭКСПЕРИМЕНТ ПРИ ОЦЕНКЕ ВНЕШНЕЙ БАЛЛИСТИКИ СТВОЛЬНОЙ АРТИЛЛЕРИИ СРЕДСТВАМИ ПРОГРАММЫ MATHCAD
ф £
Вадим Л. Хайков Ъ
независимый исследователь, г. Краснодар, Российская Федерация ш
ОБЛАСТЬ: вычислительная баллистика
ВИД СТАТЬИ: оригинальная научная статья $
ЯЗЫК СТАТЬИ: английский !
2
Резюме:
В статье представлены примененный метод математического моделирования и компьютерный эксперимент при решении задач оценки внешней баллистики ствольных систем. Разработана ~ общая схема проведения баллистического моделирования, состоящая из пяти уровней. Формальной математической з процедурой, позволяющей решить основную задачу внешней баллистики является решение задачи Коши системы :§ дифференциальных уравнений. Для 57-мм орудия ЗИС-2 оценена
ш
2
траектория полёта её снаряда. Разработан алгоритм решения § задачи и программный Ма№сасС код. и
Численное решение системы четырех баллистических ^ дифференциальных уравнений первого порядка представляет собой пятимерное пространство. Показана возможность ~ визуального представления решения в виде квадратной матрицы и определены границы подпространств. Разработан механизм |[ проверки численного решения системы дифференциальных уравнений с использованием сплайн функций. В результате проверки обнаружены эффекты незначительного увеличения погрешности на краях интервала интегрирования. Проведено сравнение решения основной задачи баллистики при помощи «мягких» и «жестких» функций-решателей. По результатам оценок «мягкими» и «жесткими» методами параметры траектории одинаковы до пятого знака после запятой.
Ключевые слова: внешняя баллистика, внешнебаллистическая траектория, математическое моделирование, вычислительный эксперимент, задача Коши, численное решение, программа МаШаС.
ш
о
™ ПРИМЕНА ПРОГРАМА MATHCAD У МАТЕМАТИЧКОМ
Ф МОДЕЛОВА^У И РАЧУНАРСКОJ СИМУЛАЦШИ ОСНОВНОГ
ПРОБЛЕМА СПО^НЕ БАЛИСТИКЕ ЦЕВНЕ АРТИ^ЕРШЕ
Вадим Л. XajKOB
° независни истраживач, Краснодар, Руска Федерац^а
ОБЛАСТ: балистика, рачунарске науке ВРСТА ЧЛАНКА: оригинални научни чланак
0
01
ос иЕЗИК ЧЛАНКА: енглески ш
^ Сажетак: о
0 У раду су представъени математичко моделоваше и < комп}утерска симулаци}а (MM&KS) у области нумеричког ^ решаваша основног проблема споъне балистике цевне
1 артиъерще. Разви]ена ¡е шема MM&KS у пет фаза за изво^еше ш балистичке симулаци]е. Показано jе да формални математички
поступак ко\и омогуЬава решете основног проблема споъне ос балистике представка нумеричко решава^е Кош^евог
проблема за систем балистичких диференцц'алних 1'едначина. Испитана jе путала лета проjектила противтенковског топа ЗИС-2 калибра 57мм. Представъени су алгоритам и код програма Mathcad ко\и су коришЬени за решава^е. Нумеричко й решете система ко\и се састоjи од четири балистичке
О 2
диференцщалне ¡едначине првог реда ¡есте тродимензионални простор. Предложена jе могуЬност да се нумеричко решеше о представи визуелно у облику квадратне матрице. Одре^ене су и
границе сваког потпростора. Разви]ен }е поступак заснован на ш „spline" функци]ама ради провере исправности нумеричког
о решеша. Приликом дате провере, примеЬени су ефекти лаког
g повеЬаша грешке на кра}евима интеграционог интервала.
> Нумеричко решеше основног балистичког проблема упоре^ено }е
помоЬу „меких" и „крутих" функци]а за решаваше. Параметри путаша испитани „меким" и „крутим" методама исти су до пете децимале.
Къучне речи: спорна балистика, путаша про]ектила, математичко моделираше, комп}утерска симулаци}а, Коши]ев задатак, нумеричко решеше.
Датум приема чланка / Дата получения работы / Paper received on: 12.10.2017. Датум достав^а^а исправки рукописа / Дата получения исправленной версии работы / Manuscript corrections submitted on: 08.11.2017. Датум коначног прихвата^а чланка за об]ав^ива^е / Дата окончательного согласования работы / Paper accepted for publishing on: 10.11.2017.
© 2018 Аутор. 06jaBMO Во}нотехнички гласник / Vojnotehnickiglasnik / MilitaryTechnicalCourier (www.vtg.mod.gov.rs, втг.мо.упр.срб). Ово je чланак отвореног приступа и дистрибуира се у складу са Creative Commons licencom (http://creativecommons.org/licenses/by/3.0/rs/).
© 2018 Автор. Опубликовано в «Военно-технический вестник / Vojnotehnicki glasnik / Military Technical Courier» (www.vtg.mod.gov.rs, втг.мо.упр.срб). Данная статья в открытом доступе и распространяется в соответствии с лицензией «Creative Commons» (http://creativecommons.org/licenses/by/3.0/rs/).
© 2018 The Author. Published by Vojnotehnicki glasnik / Military Technical Courier (www.vtg.mod.gov.rs, втг.мо.упр.срб). This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/3.0/rs/).
©CD
