Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 51, 3, pp. 511-522, Warsaw 2013 FEASIBILITY ANALYSIS OF THE MODIFIED POINT MASS TRAJECTORY MODEL FOR THE NEED OF GROUND ARTILLERY FIRE CONTROL SYSTEMS Leszek Baranowski Military University of Technology, Faculty of Mechatronics and Aerospace, Warsaw, Poland e-mail: leszek.baranowski@wat.edu.pl The work shows computational results of the effect of a simplified mathematical model of the flight dynamics of a spin-stabilized projectile on the accuracy and speed of calculation parameters for the projectile trajectory (in the entire range of operating quadrant elevation of the barrel andpossible perturbations of the conditions firing) to determine an appropriate mathematical model for the automated fire control systems of ground artillery. Key words: exterior ballistics, equations of motion, spin-stabilized projectile, fire control systems,MPMTM 1. Introduction One of the most challenging issues during the process of designing the ground (field) artillery fire control system is a model of the trajectory. The model must be simple in the context of computation (time consuming process), however very precise from the point of view of firing in real atmospheric conditions. Presented in the earlier Author’s work (2011) the projectile trajec- tory models for computation, which were based on a rigid bodymodel of the projectile motion (in short 6DOF model), are time consuming and cannot be accepted in condition of live-firing on the battle field, where time is the critical issue. Principles of classical mechanics are used in developing a 6DOF model. They are often used when constructing the equations of motion of flying objects (Gacek, 1997; Kowaleczko and Żyluk, 2009) which are an alternative to me- thods based on the principles of analytical mechanics. Among thesemethods,we can distinguish ones based on inertial generalized coordinates and referring directly to Hamilton’s principle or Lagrange equations (Koruba et al., 2010), andmethods involving the use of equations of analy- tical mechanics in quasi-coordinates e.g. Boltzman-Hamel equations (Ładyżyńska-Kozdraś and Koruba, 2012). To solve the problem of the rightmodel, NATOalliances recommend a simplified one, called the Modified Point Mass Trajectory Model (in short the MPMTM), described in STANAG 4355, a NATO official standardization agreement. The MPMTM can be used in fire control systems and to create artillery firing tables as well. However, application of the MPMTM in the fire control system requires elaboration or estimation of the specific projectile data such as mass-inertia, aerodynamic characteristics and fitting factors of the projectile: fL, i, QD, QM. The mentioned fitting factors are used to adjust the projectile impact point (elaborated from computer simulations) to live firing results (from the fire range). STANAG 4144 (Edition 2, 2005) presents the elaboration method of the fitting factors in function of the quadrant elevation angle (in short QE) taken from the firing test design. One of tests is the fall of shotmethodwhich requires aminimumof 50 rounds per velocity zone fired in 5 roundgroups at 5 quadrant elevations on a total of 10 occasions; 5with oneweapon/propellant lot combination and5with another.Then,measured coordinates of the average projectile impact points are used for calibration of themathematical model, whereas the fitting factors in function 512 L. Baranowski of the QE are approximated by a polynomial of the proper degree. STANAG 4355 (Edition 2, 1997) recommendsapolynomial of degree2 (quadratic) however, apolynomial of degree 3 (cubic) is recommended in the updated 3rd edition (2009 issue). For the need of evaluation and possible application of theMPMTM in ground artillery fire control systems, computer simulation of the projectile flight time and accuracy of its flight trajectory in real atmospheric conditions were carried out. Sets of aerodynamic characteristics and fitting factors (needed for the MPMTM) were computed and elaborated in accordance with STANAG4144 (Edition 2, 2005) procedures. Due to high cost of live firing tests, the Author used results of computer simulation of a three dimensional trajectory for the test projectile to verify the MPMTM model. The investigations were done with the following procedure: 1) Coordinates of the projectile impact point for five quadrant elevations angles (QE =350, 550, 750, 950, 1150mil) were taken from computer simulation of firing with the Denel 155mm projectile (Assegai M2000 Series) – considered as a rigid body (6DOF model). The precise description of the 6DOF model can be found in the Author’s work (2011), Carlucci and Jacobson (2008), Gacek (1997), McCoy (1999), while the values of mass- -inertia and aerodynamic characteristics for test projectile are presented in the Author’s work (2011); 2) Fitting factors, required for the simplified model, were taken in such a way as to get the compatibility results for the 6DOF with MPMTM model, obtained from computer simulation, in the context of the fire range Xend and drift yend (mass and aerodynamic characteristics for theMPMTMmodel were granted as for the 6DOFmodel); 3) Polynomials: quadratic and cubic for approximation of the fitting factors (i and fL of theMPMTMmodel) in function of QE angle, were elaborated. The simplifiedmodel was called the calibrated model (for the need of this paper); 4) Using both models: the 6DOF and MPMTM calibrated, a considerable number of firing simulations were done, including full range quadrant elevations angles (0◦ ¬ QE ¬ 80◦) at standard firing conditions as well with some perturbations. 2. The Modified Point Mass Trajectory Model TheModified PointMass Trajectory Model (MPMTM), also known as the 4 degree-of-freedom model and the Lieske model (after the American ballistician R. Lieske, 1966, who initiated its widespread usage) is described in STANAG 4355. This model represents the flight of spin- -stabilized projectiles that are dynamically stable and possess at least trigonal symmetry. Its basis is a conventional point-mass model but, in addition, the instantaneous equilibrium yaw is calculated at each time step along the trajectory so as to provide estimates of yaw, drag, drift andMagnus force effects resulting from the yaw of repose. Themathematical modeling is accomplished mainly by: (a) including only the most essential forces andmoments, (b) approximating the actual yaw by the yaw of repose neglecting transient yawing motion, (c) applying fitting factors to some of the above forces to compensate for their negligence or approximations for other forces andmoments. The vector of the yaw of repose αe (Fig. 1) has themagnitude (sinαt), where αt is the total angle of attack; αe is perpendicular to the vector velocity v and lies in the plane of drag. It is defined as follows (Lieske and Reiter, 1966) αe = iV × (x× iV ) (2.1) Feasibility analysis of the modified point mass trajectory... 513 where: iV is the unit vector in the direction of the velocity vector with respect to the air v,x – unit vector along the projectile axis of symmetry (see Fig. 1), |αe| = sinαt, for small angles, one can assume that |αe| ∼=αt. Fig. 1. The yaw of repose αe, and the total angle of attack αt 2.1. The vector form of the MPMTM In its final vector form, themathematical model of motion of the ground artillery projectiles contains the following equations (STANAG 4355 Edition 3, 2009): — dynamic differential equation of motion of the projectile center of mass mu̇=DF+LF+MF+mg+mΛ (2.2) where DF=− πρid2 8 [ CD0 +CDα2(QDαe) 2 ] v ·v − vector of the drag force LF= πρd2fl 8 ( CLα +CL α3 α2e ) v2αe − vector of the lift force MF=− πρd3QMpCmag−f 8 (αe×v) − vector of theMagnus force mg=−mg0    X1/R 1−2X2/R X3/R    − vector of the gravity force mΛ=2m(ω×u) − vector of the Coriolis force u=u0+ t ∫ 0 u̇ dt − vector of projectile velocity with respect to the ground v=u−w − vector of projectile velocity with respect to the air — dynamic equation of motion around the projectile axis of symmetry dp dt = πρd4vCspin 8Ix p (2.3) — equation for the vector of the yaw of repose αe =− 8Ixp(v× u̇) πρd3 ( CMα +CM α3 α2e)v 4 (2.4) — equation for the position of the projectile relative to the ground-fixed coordinate system X=X0+ t ∫ 0 u dt (2.5) 514 L. Baranowski 2.2. The scalar form of the MPMTM To determine the location of the projectile relative to the Earth in STANAG 4355 (2009) is used a right-handed, orthonormal, ground-fixed, Cartesian coordinate systemwith the origin of coordinates located at the gun muzzle (Fig. 2). In this coordinate system, the 1-3 plane is tangent to the earth’s surface at the launch point, the 1-axis points downrange, the 2-axis points vertically upward through the launch point, and the 3-axis points to the right, when looking downrange. Fig. 2. The ground-fixed coordinate systemwith the unit vectors (1, 2 and 3) The scalar equations of motion, which were used in the simulation model, were obtained by projecting vector equations (2.2)-(2.5) on the axes of the ground-fixed coordinate system O0123. The thus construedmathematicalmodel is a systemof differential-algebraic equations including: — scalar dynamic equations of motion for the centre of mass derived from equation (2.2) du1 dt =− πρid2 8m [ CD0 +CD α2 (QDαe) 2 ] vv1+ πρd2fL 8m ( CLα +CL α3 α2e ) v2αe1−g0 X1 R − πρd3QMpCmag−f 8m (αe2v3−αe3v2)−2Ω[sin(lat)u3+cos(lat)sin(AZ)u2] du2 dt =− πρid2 8m [ CD0 +CDα2(QDαe) 2 ] vv2+ πρd2fL 8m ( CLα +CL α3 α2e ) v2αe2−g0 ( 1− 2X2 R ) (2.6) − πρd3QMpCmag−f 8m (αe3v1−αe1v3)+2Ω[cos(lat)sin(AZ)u1+cos(lat)cos(AZ)u3] du3 dt =− πρid2 8m [ CD0 +CDα2(QDαe) 2 ] vv3+ πρd2fL 8m ( CLα +CL α3 α2e ) v2αe3−g0 (X3 R ) − πρd3QMpCmag−f 8m (αe1v2−αe2v1)−2Ω[cos(lat)cos(AZ)u2− sin(AZ)u1] — scalar kinematic equations of motion for the centre of mass derived from equation (2.6)1 dXi dt =ui i=1,2,3 (2.7) — scalar dynamic equation ofmotion around the projectile axis of symmetry derived from (2.3) dp dt = πρd4pvCspin 8Ix (2.8) Feasibility analysis of the modified point mass trajectory... 515 —equations for components of the yaw of repose vector derived from equation (2.4) αe1 =− 8Ixp(v2u̇3−v3u̇2) πρd3 ( CMα +CM α3 α2e ) v4 αe2 =− 8Ixp(v3u̇1−v1u̇3) πρd3 ( CMα +CM α3 α2e ) v4 αe3 =− 8Ixp(v1u̇2−v2u̇1) πρd3 ( CMα +CM α3 α2e ) v4 (2.9) — equations for components and magnitude of the vector of projectile velocity with respect to the air vi =ui−wi i=1,2,3 v= √ v21 +v 2 2 +v 2 3 (2.10) where: g0 =9.80665(1−0.0026cos(2lat))m/s2 – acceleration due to gravity at mean sea level, Ω=7.292115 ·10−5 rad/s – angular speed of the earth, Rz =6356766m – radius of the sphere, locally approximating the geoid, lat – latitude of the launch point; for the southern hemisphe- re lat is negative [deg], AZ – azimuth (bearing) of 1 axis measured clockwise from the true north [mil], m=43.7kg – fuzzed projectile mass, Ix =0.1444kgm2 – axial moment of inertia. In the work, the ISO Standard Atmosphere (1975) (air temperature, density, Mach number and their distribution in function of altitude) was taken as the projectile trajectory model simulation test. 3. Study of the effect of the polynomial on the precision of trajectory calculation Results of the range fitting factors (as form factor i) as well the deflection fitting factor (as lift factor fL) by mean-square approximation with polynomials of the 2nd and 3rd degree, are depicted in the figures (see below), in form of the regress line in shade of nodal points within the approximation section (for QE=350, 550, 750, 950, 1150mil) as well as beyond the approximation section (for QE = 150 and 1244.4mil). In addition, in the figures, the form of the approximation polynomial and determination correlation coefficient R2 were depicted. The closer to 1, the coefficient of the regression line is better adjusted (see STANAG 4144 Ed. 2, 2005). Two variants of firings were considered: for minimal – VK0 =319m/s (Fig. 3) andmaximal – VK0 =935m/s (Fig. 4) initial (muzzle) velocity of the projectile. Fig. 3. Results of approximations of the fitting factors (i and fL) versus quadrant elevation for the variant of firing with the minimum initial (muzzle) velocity VK0 =319m/s 516 L. Baranowski Fig. 4. Results of approximations of the fitting factors (i and fL) versus quadrant elevation for the variant of firing with the maximum initial (muzzle) velocity VK0 =935m/s For the need of visualization of the used approximated polynomial for the fitting factors approximation (i and fL), and their effect on the accuracy of the projectile trajectory (its coordinates) during a live fire test (in standard condition), simulation was conducted. The simulation was done for the simplifiedMPMTMmodel, and three following variants were taken into account: • without fitting factors, • with fitting factors approximated by the polynomial of the 2nd degree (quadratic), • with fitting factors approximated by the polynomial of the 3rd degree (cubic). Then, using the simulation of the projectile trajectory (by the 6DOF model), the following percentage errors were calculated (Figs. 5-7): — relative range error ∆rXend = Xend(6DOF)−Xend(MPMTM) Xend(6DOF) ·100% (3.1) — relative drift error ∆ryend = yend(6DOF)−yend(MPMTM) yend(6DOF) ·100% (3.2) — relative vertex error ∆rHmax = Hmax(6DOF)−Hmax(MPMTM) Hmax(6DOF) ·100% (3.3) The below presented simulations and verifications confirmed that the fitting factors are necessary in theMPMTMmodel. The factors reduce the trajectory error significantly, especially in the context of the range and drift of the impact point (Fig. 5) as well. Moreover, it can be pointed out that the approximation fitting factors by the cubic poly- nomial reduce the error more than an adequate approximation by the quadratic polynomial. Relatively significant errors in the case of the range simulation (Fig. 5b) and side deflection – drift (Fig. 6) for the barrel quadrant elevations outside the approximation section (QE = 155 and 1244.4mil) confirmed the need for the extension of live fire tests beyond STANAG 4144 Ed. 2 recommendation, it means for additional barrel angles: maximum angle for a specific cannon and for angles below 350mil. Feasibility analysis of the modified point mass trajectory... 517 Fig. 5. Computational errors of the range for the minimum initial velocity VK0 =319m/s (a) and for the maximum initial velocity VK0 =935m/s (b) Fig. 6. Computational errors of the drift for theminimum initial velocity VK0 =319m/s (a) and for the maximum initial velocity VK0 =935m/s (b) Fig. 7. Computational errors of the vertex the minimum initial velocity VK0 =319m/s (a) and for the maximum initial velocity VK0 =935m/s (b) 4. Comparison of the calculation results of the nutation angle with the yaw of repose angle For more comprehensive verification of both models (6DOF and MPMTM), in Figs. 8 and 9 the following parameters are depicted: nutation angle δ(t) for the 6DOF model and the yaw of repose angle αe(t) for theMPMTMmodel (the approximation of the fitting factors by the cubic degree polynomial was used). In the sameway as itwas done in previous investigations, two variants of firewere considered (in standardcondition):minimum VK0 =319m/s (Fig. 8)andmaximum VK0 =935m/s (Fig. 9) projectile initial velocities, both with six barrel quadrant elevations (QE =350, 550, 750, 950, 1150, 1244.4mil) of the cannon. 518 L. Baranowski Fig. 8. Comparison of the nutation angle δ(t) and the yaw of repose angle αe(t) during firings with the minimum initial velocity VK0 =319m/s Generally, except for the initial phase of the fire, where the projectile due to its angular velocity (existing in the cross section of the muzzle-barrel end) reaches a local maximum angle of nutation about 4deg, the yaw of repose angle αe is similar to the nutation angle δ (Fig. 9a-d) or goes very close (Fig. 9e,f), or can be its approximation to some extent (Fig. 8c-f). 5. Investigation of the effect of perturbations of the firing conditions on the trajectory precision calculation Finally, analysis of the computational accuracy (by the MPMTM model) of the vertex coordi- nates and impact point in the case of so-called perturbed (non-standard) conditions was done. Typical perturbations of the live firing were investigated, including: Feasibility analysis of the modified point mass trajectory... 519 Fig. 9. Comparison of the nutation angle δ(t) and the yaw of repose angle αe(t) during firings with the maximum initial velocity VK0 =935m/s • projectile mass error: ∆m=2% • initial velocity deviation: up to ∆VK0 =5% • air density error: ∆ρ=2% • difference in altitude of the target and the gun: 200m • longitudinal (range) wind, its velocity: uWg =−10m/s • side (cross) wind, its velocity: vWg =10m/s. Below, selected results of computations are presented (in form of absolute errors) for such parameters as: range, deflection, and vertex for the following variants of firing: minimum initial velocity VK0 =319m/s (Fig. 10),maximum initial velocity VK0 =935m/s (Fig. 11), five barrel quadrant elevations within the approximation sections (QE =350, 550, 750, 950, 1150mil) and one outside the approximation section (QE =150mil). 520 L. Baranowski The absolute errors such as range, deflection and vertex (marked on the graphs by the symbol ∆) were calculated on the basis of the following expressions: • absolute range error ∆Xend = |Xend(6DOF)−Xend(MPMTM)| • absolute deflection error ∆yend = |yend(6DOF)−yend(MPMTM)| • absolute vertex error ∆Hmax = |Hmax(6DOF)−Hmax(MPMTM)| Fig. 10. The effect of perturbations of firing conditions on computation errors of the range, deflection and vertex during firings with the minimum initial velocity VK0 =319m/s; (a) difference in the projectile mass, (b) difference in the initial velocity, (c) difference in the air density, (d) difference in the altitude of the target and gun, (e) the longitudinal (range) wind, (f) the side (cross) wind 6. Summary and final recommendations The following conclusions can be formulated, from the analysis of the final results: • approximation of the fitting factors by the cubic polynomial generally gives better results (reduces error) than approximation by the quadratic polynomial (Fig. 3 and Fig. 4), • relatively considerable errors for the range and deflection simulations (Fig. 6) for the barrel quadrant elevations outside the approximation section (QE =155 and 1244.4mil) pointed out the need for the extension of live fire tests beyond the STANAG 4144 Ed. 2 (2005) recommendation. Itmeans for additional barrel angles:maximumangle for a specific cannon and for angles below 350mil, Feasibility analysis of the modified point mass trajectory... 521 Fig. 11. The effect of perturbations of firing conditions on computation errors of the range, deflection and vertex during firings with the maximum initial velocity VK0 =935m/s; (a) difference in the projectile mass, (b) difference in the initial velocity, (c) difference in the air density, (d) difference in the altitude of the target and gun, (e) the longitudinal (range) wind, (f) the side (cross) wind • diagrams of angles αe and δ (Fig. 8 and Fig. 9) confirmed a satisfactory computatio- nal compatibility for both numerical models (6DOF and MPMTM) of projectile flight simulation in standard and non-standard conditions, • in most cases of live firing perturbations, the range and deflection errors (Fig. 10 and Fig. 11) do not exceed a fewmeters, whereas the highest errors in the range (up to 25m) are associated with the steep firing trajectory in strong longitudinal wind condition, • it can be concluded, from the results of investigations, that a properly calibratedMPMTM is satisfactorily accurate and gives a fast-time effective model (approx. ten times faster than the 6DOF fastest model), whatmakes it applicable to the algorithms used in the fire control computers of the ground artillery. Acknowledgement The research work was supported by the Polish Ministry of Science and Higher Education in the years 2010-2013within the grant 423/B0/A. References 1. Baranowski L., 2011, Modeling, Identification and Numerical Study of the Flight Dynamics of Ballistic Objects for the Need of Field Artillery Fire Control Systems, Military University of Tech- nology,Warsaw, p. 258 [in Polish] 522 L. Baranowski 2. Carlucci D.E., Jacobson S.S., 2008,Ballistics: Theory and Design of Guns and Ammunition, CRCPress, Taylor & Francis Group 3. Gacek J., 1997, Exterior Ballistics. Part I. Modeling Exterior Ballistics and Flight Dynamics, Military University of Technology,Warsaw, p. 352 [in Polish] 4. Koruba Z., Dziopa Z., Krzysztofik I., 2010, Dynamics of a controlled anti-aircraft missi- le launcher mounted on a moveable base, Journal of Theoretical and Applied Mechanics, 48, 2, 279-295 5. Kowaleczko G., ŻylukA., 2009, Influence of atmospheric turbulence on bomb release, Journal of Theoretical and Applied Mechanics, 47, 1, 69-90 6. LieskeR.F., ReiterM.L., 1966,Equations ofMotion for aModified PointMass Trajectory,Me- morandumReport No. 1314, U.S. ArmyBallistic Research Laboratory,Aberdeen ProvingGround 7. Ładyżyńska-Kozdraś E., Koruba Z., 2012, Model of the final section of navigation of a self- guidedmissile steeredbya gyroscope,Journal of Theoretical andAppliedMechanics,50, 2, 473-485 8. McCoy R.L., 1999,Modern Exterior Ballistics. The Launch and Flight Dynamics of Symmetric Projectiles, Schiffer Publishing 9. Procedures toDetermine theFireControl Inputs foruse in IndirectFireControlSystems, STANAG 4144 (Ed. 2), 2005 10. The ISO Standard Atmosphere, ISO 2533, 1975 11. TheModifiedPointMass andFiveDegrees of FreedomTrajectoryModels, STANAG4355 (Ed. 3), 2009 Manuscript received June 13, 2012; accepted for print July 26, 2012