0 Introduction

The traits such as mobility, agility and fierce firepower have made multiple launch rocket system (MLRS) a hot research topic for years. The problem of poor firing accuracy remains as one of the major hurdles for the development of MLRS for several years and has ever since a worldwide concern that researchers are responding to^[1].

MLRS is a complex multi⁃body mechanical system composed of multiple components connected by various kinematic pairs, and its shooting process is accompanied by complex mechanism movements, deformations and force transmissions, making it a typical multi⁃body dynamic system. Typically, the dynamics of MLRS are done by constructing a mechanical model of the system (keeping in view of the extension relations in the system), establishing and then solving the differential equation to get the dynamic response of the system.

Yang^[2] shed light on different features of MLRS including multiple launch tubes, spinning of projectile body, large length to diameter ratio, short⁃time, high thrust motor, curved trajectory and static stability design. He et al.^[3] worked on making the lightweight MLRS. For optimization, they have considered structural stiffness as an objective function. They succeeded in reducing the mass and increasing the stiffness of pedestal by 14.4% and 33.2%, respectively. Li et al.^[4] used the radial basis function neural network for vibration control in multiple launch rocket system. Zhu^[5] researched on calculations of collision force between a rocket and a directional tube by using explicit dynamics method. He obtained the collision force between the directional button and the front, the middle, the rear centering parts and the directional tube.

It is necessary to solve the accurate dynamic analysis of the multi⁃rigid⁃flexible body in order to realize the dynamic design of the performance of MLRS. At present, the main research methods of the complex multi‑body systems include Kane method based on analytical mechanics and vector mechanics^[6], Roberson⁃Wittenburg method based on graph theory^[7], Schiehlen method^[8,9] and so on. However, almost all of the multi⁃body system dynamics methods have two common characteristics: First, the overall dynamics equations of the system have to be established and need to be deduced again if the topological structure of the system changes; Second, the matrix order of the complex system is high and the calculation time is long because of the system matrix order of the system general dynamics equation is not less than two times of the degree of freedoms.

The need of fast calculation of complex multi⁃body system eigenvalue problem has been successfully addressed by Prof. Rui Xiaoting, who proposed a new multi⁃body system dynamics method namely the multi‑body system transfer matrix method (MSTMM)^[10] which avoids huge computational workload of large systems that are difficult to withstand with the finite element method and the computational “pathological” problems caused by large stiffness gradients. MSTMM has attracted much attention because of its special features like, realizing multi‑body system dynamics without the need of systematic overall dynamic equations, low order of system matrixes, fast calculation of kinetics, automatic derivation of the overall transfer equation, and high program and calculation accuracy.

The actual impact point of MLRS is difficult to be measured and requirs rapid trajectory prediction. Harlin et al.^[11] proposed a trajectory prediction method based on the state transition matrix for a ballistic missile. Khalil et al.^[12,13] established the modified point mass model for projectile trajectory prediction based on the discrete time transfer matrix method and further reconstructed the repose angle at any instant measuring time. However, most of the methods are based on the Kalman filter using three degree of freedom (3DOF) model^[14,15,16], while the usual ballistic testing does not install the angle sensor on the rocket and unable to use the six‑DOF (6DOF) model to extrapolate the impact point.

The traditional methods to improve firing performance of MLRS can be divided into two categories. The first type starts with the vehicle, optimizing the structure or adds additional devices. For example, Yan et al.^[17] obtained the best material distribution under frame through topology optimization in order to get a certain type rocket launcher with reasonable structure and lighter weight. He et al.^[18] designed and actualized vibration controlled system of launch guider of MLRS. The second type starts with the rocket, adding the guidance technology to the aircraft. For example, Qiao et al.^[19] modified the ability of canard wing rocket trajectory correction projectile.

Whereas, the first type needs to add a large number of instruments on the vehicle, which may lead to modifying the car body. The second type reduces the warhead space and power. Both types significantly increase operational costs.

In this paper, MSTMM is used to model the dynamics of a certain type of MLRS, which can quickly calculate the dynamic response and the initial disturbance. Then a fast trajectory extrapolation method is established to predict the impact point of the rockets quickly. Finally, the impact point’s correction method of MLRS is demonstrated through the closed⁃loop fire correction principle that can be achieved at low cost.

1 Dynamic Model of Multiple Lau⁃nch Rocket System

As shown in Fig.1, the dynamic model of the 40 MLRS is a multi⁃rigid⁃flexible system, which is composed of 28 rigid bodies and 40 flexible ones. The body and hinge elements are numbered. The six wheels are numbereded as 2, 5, 8, 11, 14, and 17 sequentially, and the vehicle chassis, the rigid gyration part and the rigid pitching part are numbered 19, 21, and 23 in turn. The tail of the l th ( l=1,2,⋯,40 is the directional tube number where the rocket is being fired) directional tube is numbered as 21+5l . The part number between the l th directional tube’s front and back support frame is numbered as 22+5l . The number of the front part of the support frame is numbered as 23+5l . The elements 1, 3, 4, 6, 7, 9, 10, 12, 13, 15, 16, 18, 20, 22, 19+5l and 20+5l are the elastic hinges of the coupling body elements’ the ground or body elements, and the body elements. The elements 2, 5, 8, 11, 14, 17, 19, 21, 23 and 21+5l are rigid bodies. The elements 22+5l and 23+5l are flexible beams.

Fig.1 Launch dynamics model of a new MLRS

The system has a total of eight boundary points, including the contact point between the six wheels and the ground and the free boundary at both ends of the directional tube. The boundary point number is 0.

According to the automatic deduction theorem of overall transfer equation and MLRS model topology, as shown in Fig.2, the multiple rockets system transfer equation can be automatically written as^[20]

Fig.2 Topological diagram of a multiple launch rocket dynamics model

where

where O is a zero element; C=I6O6×6O6×6-I6 ; T _{a - b} is the sequential continuous product of all element transfer matrices on the path from the root point a to the system tip point b; G _{c - b} the sequential continuous product of all element transfer matrices on the path from the element point c to the system tip point b.

According to Eq.(1), the unknown state vector Z20+5l in Zall* can be further eliminated and get

Applying the boundary conditions for multiple launch rockets to Eq.(2), removing the zero element in Zall and results in Z¯all , removing the column in Uall corresponding to the zero element in Zall , and resulting in a square matrix U¯all , we get Eq.(4) from Eq.(2)

Then the characteristic equation of the multiple rockets system is

The natural vibration frequency ωk(k=1,2,3,⋯,n) of MLRS can be obtained by solving Eq.(5). Under the given normalization conditions (for example, the element with the maximum absolute value of all elements in Z¯all is equal to 1), the solution of Eq.(4) can obtain Zall corresponding to the natural vibration frequency ωk . By using the transfer equation of the elements, the state vectors of ωk corresponding to the point of the directional tube are obtained, and the mode shapes of each element in the k th order mode of MLRS are obtained.

According to MSTMM^[21], the first 14 order natural frequencies of MLRS with elevation and azimuth angle of 50°and 0° are calculated under the full load condition, as shown in Table1. The matrix order and computation time of the two methods are listed in Table2which shows that the calculation speed of MSTMM is 2 806.7 times faster than the finite element (FEM) method under the same conditions.

MSTMM is used for fast dynamic modeling of MLRS, and the dynamic response and the initial disturbance of MLRS can be quickly obtained with the launch dynamics equations. Then the initial disturbance is substituted into the 6DOF rigid rocket external ballistic equations^[1]. The impact points and the firing accuracy of the MLRS are calculated by using Monte Carlo simulation. We establish a numerical simulation system for launch and flight dynamics of MLRS. The simulation flow chart is shown in Fig.3.

Fig.3 Flowchart of launch and flight dynamics simulation system for MLRS

2 Ballistic Trajectory Extrapolation for MLRS

The large rocket extrapolation method uses GPS/IMU to locate positions, velocities, and angles of the rocket after launch, the Kalman filter and other filtering methods to obtain ballistic parameters, and 6DOF external ballistic equations to calculate the impact points. However, due to multi⁃factors like high temperature, high pressure, high overload and cost during firing multiple rockets, it is hardly to install angle sensors on all rockets in MLRS and to test positions by historical radar data of rockets. The angle information cannot be measured and therefore the 6DOF equations is no longer applicable.

In order to improve the accuracy of extrapolation of multiple rockets for using 6DOF equations, a method of ballistic extrapolation for MLRS is established by using the Davidon‑Fletcher‑Powell (DFP) method. Through historical radar data, the angular information of the passive section of the rocket is obtained quickly. Then the impact points are extrapolated by the 6DOF ballistic equations. The results of the extrapolation points agree well with actual ones.

2.1 Equations of motion of the ordinary rocket

Considering the dynamic unbalance, mass eccentric and gust of wind, the equations of motion of the ordinary rocket can be written as

where x,y,z are the coordinates of the rocket geometry center in the ground coordinate system； vx,vy,vz the components of the velocity of the rocket geometry center in the ground coordinate system； m is the mass of the rocket； A the equatorial rotational inertia of the rocket； C the polar moment of inertia of the rocket； g the gravitational acceleration； Fp the thrust force of rocket, Ir the specific impulse of engine； vr the relative atmospheric velocity； wx the velocity of the longitudinal wind； wz the velocity of the cross wind； v˙r the relative acceleration； Rx the drag force； Ry the lift force； Rz the Magnus force； ψr the relative deflection angle； θr the relative ballistic inclination angle； δr the relative angle of attack； δr1 the projection of the relative angle of attack on the vertical plane and δr2 the projection of the relative angle of attack on the lateral plane； ψ˙r the relative declination angle velocity； θ˙ the trajectory inclination angle velocity； αr the relative aiming angle； kxD the polar damping moment coefficient； kzD the equator torque damping coefficient； Mz the Magnus moment； Mz the stable moment； βDη and βDζ are the projections of dynamic unbalance angle βD in axis η and ζ of the o1ξηζ , respectively； Lmη and Lmζ are the projections of mass eccentricity Lm in axis η and ζ of the o1ξηζ , respectively^[1].

2.2 Algorithm for estimating rocket attitude angles

When computing the 6D trajectory model of the rocket, we need to provide 12 independent initial values including (x,y,z,vx,vy,vz,γ,ϕa,ϕ2,γ˙,ϕ˙a,ϕ˙2) . However, only the first six (x,y,z,vx,vy,vz) can be obtained from the radar test data, and the initial value of the spin angle γ has no influence on the motion of the rocket. In this paper, we reconstruct the angular velocity of spinning, two angles of oscillation, and two angular velocities of oscillation (γ˙,ϕa,ϕ2,ϕ˙a,ϕ˙2) will be the matter that it requires to be researched.

If the radar data at time instant ti and ti+1 are known, the positions of a rocket is expanded using a second⁃order Taylor series at ti as

where Xi and Xi+1 are the position coordinates of the rocket at ti and ti+1 ， respectively; Vi is the velocities of the rocket at ti ； Δt=ti+1-ti the difference between ti and ti+1 ； Ai the accelerations of the rocket at ti .

Further, it can be deduced that Ai=2Δt2(Xi+1-Xi-ViΔt) . Ai can also be expressed as dvxdt,dvydt,dvzdti , where (x,y,z,vx,vy,vz)i at ti and (x,y,z,vx,vy,vz)i+1 at ti+1 can be obtained from radar measurement. On the right side of 6DOF equations, only the angular velocity of spinning and the two angles of oscillation (γ˙,ϕa,ϕ2)i are unknown and independent, therefore, dvxdt,dvydt,dvzdt can be considered as a nonlinear function with respect to (γ˙,ϕa,ϕ2)i at ti , namely, Fi(γ˙,ϕa,ϕ2)=Ai-2Δt2(Xi+1-Xi-ViΔt) . In order to solve (γ˙,ϕa,ϕ2)i , the problem is transformed into solving a ternary nonlinear equations as follows

There are many solutions to the system of nonlinear equations. Typically, the nonlinear equations are numerically solved by iteration methods. Common methods include the fixed point iteration, Newton iterative and Quasi Newton iterative.

In order to improve the precision of ballistic trajectory prediction, this paper proposes an algorithm to estimate the rocket attitude angles based on the 6DOF trajectory model. Hereby, the algorithm utilizes DFP method in solving nonlinear equations and radar trajectory test information containing only position coordinates of the rocket to reconstruct the angular displacements.

DFP formula provide the solution to the secant equation that is the closest to current estimate and satisfies the curvature condition. It is the first quasi⁃Newton method to generalize the secant method to a multi‑dimensional problem. This update maintains the symmetry and positive definiteness of the Hessian matrix^[22].

DPF avoids obtaining inverse matrix in its iteration formula, thus increases the computation speed^[23].

In a similar way, the radar data at ti+1 and ti+2 , (γ˙,ϕa,ϕ2)i+1 can be calculated. Hereby, (ϕ˙a,ϕ˙2)i at ti can be approximately evaluated by

2.3 Simulation and test for ballistic trajectory extrapolation method

A ballistic trajectory is calculated by the 6DOF model in order to validate the attitude angle algorithm， because the radar test data does not include angle information. The trajectory in 3—4 s of the calculating position information is applied to the algorithm to estimate the attitude angle information. As shown in Figs.4—6, the difference between the estimated and the 6DOF⁃calculated is very small.

Fig.4 Pitch angle of the estimated and the 6DOF during flight

Fig.5 Yaw angle of the estimated and the 6DOF during flight

Fig.6 Spin rate of the estimated and the 6DOF during flight

The impact points can be obtained by substituting the flight parameters of the last points of the rockets into the 6DOF ballistic equations.

Define the impact point error (%) as

The errors between the extrapolation and the real trajectory impact point are listed in Table3.

The real flight data validate that the proposed algorithm shows good prediction effect using the 6DOF trajectory model.

3 Correction of Firing Precision for MLRS

3.1 Firing correction algorithm

Firing accuracy is an important technical indicator of MLRS, including two aspects: Firing precision and firing dispersion. Firing precision is characterized by the deviation of the mean impact point of the rocket relative to the target; while the firing dispersion indicates the density of the impact point for the mean impact point.

The point at which a rocket hits a target, land or air blast (whether or not it hits a target) is called the impact point. The average of impact points of MLRS is called the mean impact point. A virtual trajectory corresponding to the mean impact point is called the average trajectory. It reflects the average of multiple trajectories, as shown in Fig.7.

Fig.7 Average trajectory of multiple launch rockets

Due to various environmental factors in actual firing of rockets, there is an error between the mean impact point and the target point. Many scholars use guidance methods on rockets. However, due to financial constraints, low⁃cost firing corrections are widely used. According to the closed⁃loop fire correction principle, this paper considers correcting the range of the longitudinal X⁃direction and the drift of the lateral Z⁃direction of the ballistic error on the geodetic coordinate system XOY on the same ground plane.

（1） Range

Elevation angle can be adjusted by using Newton⁃Raphson method for the iterative calculation.

If the average rate of change f(θk)-f(θk-1)θk-θk-1 of function f(θk) is used to approximate the instantaneous rate of change f'(θk) , the iterative format for the elevation angle that converges to the mean point of impact is

where θk-1 ， θk are the elevation values at the begin of the iteration； θk=θk-1+Δθ , Δθ is the iteration step at the beginning of iteration； f(θk-1) ， f(θk) ， f(θM) are the corresponding ranges of elevation angles θk-1 , θk and the mean point of impact M， respectively.

（2） Drift

Azimuth angle can be adjusted by the angle between two vectors of the target point and the impact point of the final iteration

where A is the impact point of changed the elevation angle； O the origin point of the coordinate； T the target point, and ∠AOT is the correction to the azimuth.

The correction steps of elevation and azimuth angles are as follows:

(1) Input the rocket data, aerodynamic data, wind data, and initialize the firing parameters.

(2) Given elevation and azimuth, use Monte Carlo launch dynamics calculation and 6DOF ballistic equations to get the mean point of impact.

(3) Determine the value of the elevation step and calculate the corresponding ranges f(θk-1) , f(θk) of iteration at starting points θk-1 , θk ； calculate the elevation iterative; end the iteration of the elevation angle when the error of the iteration of the range is less than the given error.

(4) Calculate the azimuth angle. When the iteration error of the drift is less than the given error, end the iteration of the azimuth angle.

(5) Record the correction of elevation and azimuth, and analyze the error situation.

The flow chart of elevation and azimuth angle correction algorithm is shown in Fig.8.

Fig.8 Flowchart of elevation and azimuth angle correction algorithm

3.2 Results of firing correction

In the actual firing test, the elevation and azimuth angles are (19.8°, 0°). The initial mean impact point of Monte Carlo simulation before iteration is (13 811.202 3 m, -147.312 0 m). The coordinates of the target point are (14 000 m, 0 m). The corrections of elevation and azimuth are obtained by putting the error between the mean point of impact and the target point into the angles iterative formulas. The results of the iteration are shown in Table4.

The mean impact point after correction is (14 000.898 5 m, -0.540 6 m), and the corresponding elevation and azimuth angles are (20.356 6°, 0.588 7°). The corrections of the elevation angle and azimuth angle are (0.356 6°, 0.588 7°). The error between the mean impact point after corrected and target point meets the engineering error, which effectively improves the firing precision of MLRS.

4 Conclusions

The firing angle of MLRS is improved by keeping the cost as low as possible. First, the dynamic model of MLRS has been built using MSTMM, and the corresponding dynamic characteristics, the initial disturbance and the impact points are obtained. Then a method for predicting the angle information of rockets is established, which can quickly and effectively extrapolate the impact points by using radar test data. Finally, the closed⁃loop fire correction principle is employed for adjusting the elevation and the azimuth angles to improve the firing precision of the MLRS at low cost. This research can provide a theoretical basis and technical means in the field of MLRS to significantly enhance the firing accuracy.

References