CC BY
82
MSC 80A20, 65D25
DOI: 10.14529/mmp150106
CONSTRUCTION OF APPROXIMATE MATHEMATICAL MODELS ON RESULTS OF NUMERICAL EXPERIMENTS
V.A. Tenenev, Kalashnikov Izhevsk State Technical University, Izhevsk, Russian Federation, tenenev@istu.ru,
I.G. Rusyak, Kalashnikov Izhevsk State Technical University, Izhevsk, Russian Federation, primat@istu.ru,
V. G. Sufiyanov, Kalashnikov Izhevsk State Technical University, Izhevsk, Russian Federation, vsufiy@mail.ru,
M.A. Ermolaev, Kalashnikov Izhevsk State Technical University, Izhevsk, Russian Federation, michael6@yandex.ru,
D.G. Nefedov, Kalashnikov Izhevsk State Technical University, Izhevsk, Russian Federation, denisnefedov1@yandex.ru
A mathematical model of an artillery shot is represented as a system of non-stationary one- and two-dimensional differential equations of the multiphase gas dynamics and heat transfer. Conjunction Euler-Lagrange method is used for the numerical solution of gas-dynamic equations. The initial mathematical model is approximated by a system of ordinary differential equations using a vector of correction functions. Correction functions are found from solutions of multiobjective optimal control problem. Multiobjective optimization is carried out using a hybrid genetic algorithm. The resulting model is adequate and allows doing more processing series of calculations the main process parameters (projectile velocity and maximum pressure) depending on the input parameters. Comparative analysis of different approximators (linear multiple regression, support vector machines, multi-layer neural network, radial network, the method of fuzzy decision trees) showed that an acceptable accuracy 0,4-0,5 % is provided by only non-linear approximation methods, such as multi-layer and radial neural networks. Constructed approximate models are not require much computing time and can be implemented in the control systems.
Keywords: mathematical model of the shot; multi-phase gas dynamics; approximate models; multi-objective optimization.
Introduction
The present mathematical modeling level allows us to construct the models of very complex processes and to make a detailed numerical experiment. As a rule, the numerical implementation of such mathematical models involves time-consuming. In some cases, calculation results of processes characteristics are required in real time. This can be the movement control of any objects or any process. Having information on the results of numerical or natural experiment we can construct mathematical models with lower dimension, which reproduce the characteristics of these processes by small time consuming. We consider the example of mathematical models dimension reduction. This example is related with numerical modeling of artillery shot, which is characterized by the diversity of occurring physical, chemical and gas-dynamic phenomena. In various barrel systems tests the operative analysis of a natural experiment is very important in order to improve the accuracy characteristics and identification of using calculation methods.
1. Mathematical Model of the Shot
The movement of burning gunpowder elements along the barrel at the shot is a classic case of heterogeneous system movement, due to processes of heat and mass exchange and friction between the phases [1].
The main assumptions of the model:
1) distances, at which the flow parameters vary considerably, much more than distances between the particles and size of the particles;
2) various phases are present simultaneously at all points of space, at the same time each phase occupies a portion of the mixture;
3) calculation of each phase movement can be performed independently from the mixture under the condition that appropriated account of the interaction between the different phases is provided;
4) viscosity and thermal conductivity is significant only in the phases of interaction processes;
5) particles have the same size on average, collisions, or interaction between the particles, can be neglected;
6) phases movement is one-dimensional;
7) heat transfer to the burning surface of the grains is not taken into account (heat wave movement velocity in the gunpowder equals to burning velocity);
8) material of the particles is incompressible;
9) gas parameters inside and outside of the gunpowder elements are same in this cross-section.
The system of one-dimensional hydro-mechanical equations describing the barrel system ballistics in terms of mechanics of heterogeneous medium has the form [1]:
dpmS dpmSv tt. I ^ = SG, dt dx
d6 (1 - m) S d6 (1 - m) Sw _
tt. 1 ^ = —oG,
dt dx
dpmSv dpmSv2 dp
—^--1--^-= -mS — + SGw - Stw - ncTic,
dt dx dx
d6 (1 - m) Sw d6 (1 - m) Sw2 . ^dp ^ ^
+ —-—tt-^-= - (1 - m) S^- - SGw + Stw - ncT2c, (1)
dt dx dx
(v — w)2
dpm.Se dpmSev d [mSv + w (1 — m) Sw] тт. I ^ = p ^ +
dt dx dx
Q +
I
+STw (v - w) ncTicv - ncqc,
d$ d^ S0 lit + wdxx = A0a ^ u,
p (1 - ap) = (k - 1) pe,
where t is time; x is coordinate; p is gas density; m is mixture porosity (volume of voids per unit volume); S is variable cross-sectional area of the chamber and the barrel; v,w are movement velocities for gas and solid phase in the barrel respectively; G is gas coming per unit volume; 6 is density of gunpowder material; p is pressure; tw is hydraulic resistance
of gunpowder elements per unit volume; nc is perimeter of the barrel; r\c, t2c are friction force of gas and gunpowder particles on the barrel surface per unit area respectively; £ is internal energy per unit mass of gunpowder gases; Q is heating value (potential) of gunpowder; a is covolume; qc is heat flow to the barrel surface; — is the relative proportion of burnt gunpowder; A0, S0 are the initial volume and surface of grain gunpowder; a (—) is current burning surface to original ratio; u is the linear velocity of gunpowder burning; k = — is adiabatic index; Cp, Cv are heat capacity of gas at constant volume and pressure respectively.
Received system of equations contains seven variables: gas parameters p, p,£,v, solid phase velocity w, the mixture porosity m and relative share of burnt gunpowder —.
Initial and boundary conditions:
f
0 < x < Lkm, t = 0, v = w = 0, p = pn, £ = £n = 77, — = —n;
u
x = 0, t > 0, v = 0;
x xcn, t > 0, v vcn •
Boundary movement x = xcn and projectile velocity are determined by the equations of projectile movement:
q—TT = Scn (p - Ppr) - Ft, (2)
dt
dt
= Vc
where Lkm is length of the chamber; f = RTv is gunpowder force; R is gas constant; Tv is the gunpowder burning temperature at the constant volume; U = k — 1; q is projectile mass; Ft is resistance force during the projectile movement in the barrel; Scn is area of the barrel; ppr is compressed air pressure in front of the projectile (back pressure); vcn is projectile velocity; xcn is the path of the projectile in the barrel; pn,£n,—n are the initial values of the relevant parameters.
Since the share of friction and heat transfer to the surface of the barrel in the total balance of momentum and energy of the shot is small, we can use the approximate dependencies neglecting the contact heat transfer and friction of particles on the barrel:
qc = Nu A (T - Tc) dk
Tc = Tic = 0.125£pv |v| , T2c = 0,
where Nu is Nusselt number; A is thermal conductivity coefficient of gunpowder gases; T is temperature of the gunpowder gases; Tc is surface temperature of the barrel, dk is barrel diameter; £ is resistance coefficient of the gunpowder grains in the layer. Surface temperature of the barrel is defined by the equation:
dn ~dt
2 (NuA\
CcÔcAc
V dk )
(T - Tn-^nf, n (0) = 0,
where n = (Tc — Tn), Tn is initial temperature of the barrel surface; cc,Sc,\c are heat capacity, density and thermal conductivity of barrel material.
The air initially located in the chamber and the heating particles entering the flow with the gaseous igniter burning products and creating considerable heat flows deposited
2
on the gunpowder surface are taken into account to describe the ignition process. In the gas-dynamic system of equation, which takes into account the gradual ignition, it is additionally accepted:
1) burning products of various parts of the charge, igniter and air are the homogeneous non-reactive mixture with the same velocities, pressures and temperatures;
2) the gases nature for different brands of igniter is the same;
3) grains of the igniter move with gunpowder elements;
4) inner and outer surfaces of the gunpowder elements in a given section are ignited at the same time;
5) deposition of incandescent particles in ignited part of the charge is neglected;
6) gunpowder burning velocity is determined on the basis of approximate methods of non-stationary and erosive burning [1].
2. Numerical Solution of the Model Equations
Conjunction Euler-Lagrange (CEL) method [1] is used for the numerical solution of equations (1). The scheme of this method belongs to the class of homogeneous conservative schemes that by introducing pseudoviscosity allow to keep "through" account of gas-dynamic parameters in the shock waves presence. Pseudoviscosity component included in the pressure and "smears" the shock wave front on several intervals of grid, so that the values of the functions vary continuously passing through the leap and satisfy conservation law of Rankine-Gyugonio. The upstream differences are used to approximate the convective terms. Due to shifted grids, CEL-method has the second order of accuracy, which ensures fast convergence of the solution regarding to the step along the coordinate.
To simplify the model (1), (2) we consider the processes occurring in the shot, by assumption that pressure, density and temperature of gas are taken average over the burning chamber volume. Accounting of combustion products movement will be made on the bases of some functions selected on the results of numerical simulation for solutions of equations (1), (2). For average over volume parameters we have the following system of ordinary differential equations [2]:
r^rdp kR
W-f- = --
dt 1 — ap
M (xTp — T ) +
kp
P (1 — ep)
M (1 — p ) — ScnVcnPc
Wp aT kR
T at 1 — ap
M (xTp — T ) +
dp
p (1 — ep)
de dt
= u,
dW ~dt~
S v + M
ucnucn I r- ) 0
M (1 — 0 ) — ScnVcnpc
(3)
as well as equation (2).
Here W is free volume of the chamber and the barrel; Tp is isobaric temperature; x is heat loss coefficient; pc is density behind the projectile associated with pressure pc in
back of projectile space pc = -rt ; M = ôuS is gas coming from combustion products;
a + —
Pc
S (e) is relationship between the burning surface area and thickness of the burnt body;
gunpowder burning velocity is determined by the dependence u = «ip; u is unit burning velocity.
We introduce the vector of correction functions $ = (e) , (t) , (t) , ^u (e)]T, responsible for the approximation of the model (3) to the mathematical model (1):
E = + ; (4)
x = (5)
u8(1 - fi (x))
PPc
J2 cn
pc =-^-; (6)
1 + 2f(x)
u = ui^up; (7)
82 - 1
fi (x) = 778T-X2 • (8)
(8 + x)
Here 8 = ; Wk = SkLk is combustion chamber volume; Sk, Lk are sectional area
L
and length of the combustion chamber; x = cn-k.
Lk
The controls have meaning for e < ek and therefore are accepted depending
on the size of the burnt body. The value of Sq = — in (4) is determined by the mass of
ek
charge u and body thickness ek.
During burning it is assumed that igniter has the same burning velocity as gunpowder, and the area of the burning surface depends on the thickness of the burnt body
6uw / 2e x 2
(e) = 7D 1 - _
0m Dw \ D
where uw, ôw, Dw are mass, density, diameter of igniter particulars. The typical value = 0,9. Initial values for equations (3), (8): p (0) = pn,T (0) = Tn,e (0) = 0, W (0) = Wn, Vcn (0) = 0.
Dependencies P (t) = p (t) , Pc (t) = pc (t) , VCn (t) = vcn (t) obtained by numerical solution of equations (1), (2) are used to calculate the functions Functions $ should ensure the minimum values of the deviations of variables p (t) , pc (t) , vcn (t), calculated for equations (3), (8) from values P (t) , Pc (t) , Vn (t) for the time of the shot tk:
ftk 2 Jp ($) = [p (t) - P (t)]2 dt ^ min,
Jo *
/•tk 2
Jc ($)= / [pc (t) - Pc (t)]2 dt ^ min, (9)
Jo *
/•tk 2
Jv ($) = [vcn (t) - Vcn (t)] dt ^ min .
Jo *
After substituting expressions (4) - (7) in equations (3), (8) we obtain a system of ordinary differential equations for the phase variables X = (p, T, e, W, vcn)T with the control The system of equations (3), (8) together and conditions for minimum of
functionals (9) represent the optimal control problem. The system of equations (3), (8) is solved numerically by the Runge-Kutta method with fourth-order accuracy with a step At. Controls $ are replaced by linear functions on intervals of time and the burnt body [tj, tj+i] , to = 0, = tk; [ej, ej+1] j eo = 0, eN$ = ek:
$j+i - $
$ (y) = $j + -^ (y - yj) , j = 0,N* - 1, y = (t, e) ,
Vj+i - yj
where N$ is number of intervals. Thus the optimal control problem is reduced to the finite-dimensional problem of multicriteria optimization to determine the values $j. Multicriteria optimization problem is solved by use of a hybrid multipopulation genetic algorithm [3].
We consider the following algorithm, which does not require the introduction of additional sub-populations and user interaction in the choice of optimal Pareto solutions
[3].
Suppose we are given the problem of multi-criteria optimization:
F (x) ^ min,
or
Fi (x) ^ min, i = 1, N.
Just as in the case of one criterion the population by given size R is formed. From this population the individuals (solutions x) are chosen, that are the bests for each of N criteria with criteria values
F0 (x) , i =
Individual, a leader in a given population, is defined by the rule:
x0 = arg
min I max I Fj (x) — F0 (x) I I
i=1,R \i=1,N J
(10)
Selecting for crossing is carried out by the tournament. Solution obtained in a result of a number of iterations is unique and Pareto optimal. Genetic algorithm with a real encoding is used as an optimization method.
Fig. 1 shows relationships P (t) , Pc (t), obtained from the numerical solution of equations (1), (2) and used as a standard for the selection of corrective functions Ф (black circles show pressure on the bottom of the channel, light circles show pressure on the bottom of the projectile). Fig. 2 shows dependence Vcn (t) (black circles). Relations P (t) , pc (t) , vcn (t), calculated on equations (3), (8) according to corrective functions, are shown at the same figures (bold line in fig. 1 shows pressure on the bottom of the channel, light line shows pressure on the bottom of the projectile). Velocity of the projectile vcn (t) in fig. 2 is a bold line.
As you can see, the selection of the corrective functions provides reproduction of firing process in the basic details. Change of correction functions is shown in fig. 3: ^is
dependence 1; (px{jQ is dependence 2; p^^Jr^j is dependence 3; is dependence 4.
Let us consider the following level of dimension reduction of mathematical model -approximation of point features. On the example of artillery shot it is the dependence of the maximum pressure pmax in the chamber and muzzle velocity v^f* from the initial process parameters (vector z = (pmax,vmnx)T). As the impact parameters we consider
P, Pc, MPa 500 n-
450 -
0 0.005 0.01 0.015 0.02 0.025
Fig. 1. Pressure curves at the bottom of the channel and the bottom of the projectile
Vcn, m/s
0 0.005 0.01 0.015 0.02 0.025
Fig. 2. Velocity of the projectile
vector Y = (Tp, ek, U\, S, Xo) in the range [0.95; 1.05] of nominal values. Under the random distribution of the parameters Y from the solution of equations (3), (8) for given values
Fig. 3. Correction functions
of the correction functions $ we obtained H = 1000 points (z, Y)j , j = 1,H. On these data we carry out comparative analysis of different approximators:
1) multiple regression equation;
2) support vector method;
3) multilayer neural network;
4) radial network;
5) the method of fuzzy decision trees.
Multiple regression equation is z = b + (w, Y), where the coefficients b, w are found by the method of least squares.
In support vector method z (Y) = EH=1 (Ah - Ah) K(Y, Yh) + b, summation is carried out for (Ah - Ah) = 0, where K(Y, Y') = exp (-0, 5(Y - Y')Ta-1 (Y - Y')) is a kernel function. Coefficients Ah, Ah, which are Lagrange multipliers in the quadratic programming problem, are calculated iteratively [4] on available data set, a is a covariance matrix.
In multilayer neural network the input signal is converted by the layers according to the expressions:
Nk-i
uk =
Y, wj j1, qk = g(uk), i = l, Nk, k = l, Kc, g(u) =
i
j=0
l + exp(-ви) '
Kc
q0 = Y, z = q
where Nk is the number of neurons in the kth layer; Kc is the number of layers in the neural network; w j are weighting coefficients; g (u) is activation function with optimization parameter ^.
Radial network performs conversion z (Y) = wh exp
2-h
where L is the
number of centers of radial neurons with values Ch; coefficients wh, ah are determined by training.
In the method of fuzzy decision trees the decision tree [3] is constructed on a data set. Constructed decision tree is considered as a set of fuzzy rules by the form
Rr : i/flY £ Ajr then z is Br, r =1, KR.
Condition Yj £ Ajr corresponds to the separation of the set of points on the breakdown
Yj < Wj and means entering of value Yj in the fuzzy intervals on the left and right of
parameter Wj with the given membership functions. For a given vector Y the degrees of truth for each rule ar, r = 1,KR, are defined. As a result, an aggregated output signal is determined by formula:
1 Kr / n
z (Y) = - J2ar I Wr0 + Y1 WrjYj
l^r=l ar r=l \ j=l
Coefficients wrj are determined by available training selection using the procedure of pseudoinversion.
Selection of 1000 points is divided into the training selection (750 points) and the test selection (250 points). On the test selection the standard error was calculated for each approximator. The results are given in table. The best result was shown by the radial
Table
The standard error
Approximator pmax p „.max vcn
1 10,43 5,07
2 3,23 4,32
3 3,88 4,45
4 2,28 4,05
5 6,89 5,05
network. Error for maximum pressure pmax is 2,28 MPa or 0,52%, for muzzle velocity error is 4,05 m/s or 0,43%. Number of centers of radial functions is L = 100, ah = 2. The linear dependence (1st approximator) showed the worst results. This is well illustrated by the graphs for comparing of the approximated value pmax with pmax calculated from equations (3), (8). Fig. 4 corresponds to the linear approximator, and fig. 5 corresponds to the radial approximator.
Conclusion
Application of the methods of knowledge and data extraction, and genetic algorithms allows us to construct the adequate mathematical models that can be implemented in various control systems.
Acknowledgements. The research was done under the state assignment No 1.1481.2014/E, registration number is 114072170012 from 21.07.2014.
Fig. 4. Comparison of the approximated and calculated pressures (linear approximator)
Fig. 5. Comparison of the approximated and calculated pressures (radial approximator)
References
1. Rusyak I.G., Ushakov V.M. Vnutrikamernye geterogennye protcessy v stvol'nykh sistemakh [Intrachamber Heterogeneous Processes in Barrel Systems]. Ekaterinburg, UrO RAN, 2001. 259 p.
2. Sorkin R.E. Teoriya vnutrikamernykh protcessov v raketnykh sistemakh na tverdom toplive: vnutrennyaya ballistika [Intrachamber Processes Theory in Missile Systems for Solid Fuels: Internal Ballistics]. Moscow, Nauka, 1983. 288 p.
3. Tenenev V.A., Yakimovich B.A. Practice of Genetic Algorithms. Universitas GYOR Nonprofit Kft., 2012. 279 p.
4. Kecman V. New Support Vector Machines Algorithm for Huge Data Sets. Lektcii po neyroinformatike. Po materialam Shkoly-seminara "Sovremennye problemy neiroinformatiki" [Lectures on Neuroinformatics. On Materials of School-Seminar "Modern Problems of Neuroinformatics"]. Moscow, 2007, pp. 97-176.
Received November 14, 2014
УДК 533.6.011 + 519.653 DOI: 10.14529/mmp150106
ПОСТРОЕНИЕ АППРОКСИМИРУЮЩИХ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ ПО РЕЗУЛЬТАТАМ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ
В.А. Тененев, И.Г. Русяк, В.Г. Суфиянов, М.А. Ермолаев, Д.Г. Нефедов
Математическая модель артиллерийского выстрела представлена в виде системы нестационарных одно- и двумерных дифференциальных уравнений многофазной газодинамики и теплообмена. Для численного решения газодинамических уравнений используется совместный эйлерово-лагранжев метод. Исходная математическая модель аппроксимируется системой обыкновенных дифференциальных уравнений с применением вектора корректирующих функций. Корректирующие функции находятся из решения многокритериальной задачи оптимального управления. Многокритериальная оптимизация осуществляется с применением гибридного генетического алгоритма. Полученная модель является адекватной и позволяет провести большую вычислительную серию расчетов основных параметров процесса (скорости снаряда и максимального давления) в зависимости от исходных параметров. Сравнительный анализ различных аппроксиматоров (линейная множественная регрессия, метод опорных векторов, многослойная нейронная сеть, радиальная сеть, метод нечетких деревьев решений) показал, что приемлемую точность 0.4-0.5% обеспечивают только методы нелинейной аппроксимации, такие как многослойная и радиальная нейронные сети. Построенные аппроксимирующие модели не требуют больших затрат вычислительного времени и могут быть реализованы в системах управления.
Ключевые слова: математическая модель выстрела; многофазная газодинамика; аппроксимирующие модели; многокритериальная оптимизация.
Литература
1. Русяк, И.Г. Внутрикамерные гетерогенные процессы в ствольных системах / И.Г. Русяк, В.М. Ушаков. - Екатеринбург: Изд-во УрО РАН, 2001. - 259 с.
2. Соркин, Р.Е. Теория внутрикамерных процессов в ракетных системах на твердом топливе: внутренняя баллистика / Р.Е. Соркин. - М.: Наука. Главная редакция физико-математической литературы, 1983. - 288 с.
3. Tenenev, V.A. Practice of Genetic Algorithms / V.A. Tenenev, B.A. Yakimovich. -Universitas GYOR Nonprofit Kft., 2012. - 279 p._
4. Kecman V. New Support Vector Machines Algorithm for Huge Data Sets // Лекции по
нейроинформатике. По материалам Школы-семинара «Современные проблемы нейро-
информатики>. - Москва, 2007. - С. 97-176.
Валентин Алексеевич Тененев, доктор физико-математических наук, профессор, кафедра «Высшая математика>, Ижевский государственный технический университет имени М.Т. Калашникова (г. Ижевск, Российская Федерация), tenenev@istu.ru.
Иван Григорьевич Русяк, доктор технических наук, профессор, заведующий кафедрой «Математическое обеспечение информационных систем>, Ижевский государственный технический университет имени М.Т. Калашникова (г. Ижевск, Российская Федерация), primat@istu.ru.
Вадим Гарайханович Суфиянов, кандидат физико-математических наук, доцент, кафедра «Математическое обеспечение информационных систем>, Ижевский государственный технический университет имени М.Т. Калашникова (г. Ижевск, Российская Федерация), vsufiy@mail.ru.
Михаил Александрович Ермолаев, инженер-программист I категории, кафедра «Математическое обеспечение информационных систем>, Ижевский государственный технический университет имени М.Т. Калашникова (г. Ижевск, Российская Федерация), michael6@yandex.ru.
Денис Геннадьевич Нефедов, инженер-программист, кафедра «Математическое обеспечение информационных систем>, Ижевский государственный технический университет имени М.Т. Калашникова (г. Ижевск, Российская Федерация), denisnefedov 1 @yandex. ru.
Поступила в редакцию 14 ноября 2014 г.
