The Effect of Manufacturing Quality on Rocket Precision

: The effect of manufacturing quality on rocket impact point dispersion is analyzed. The approach presented here applies to any type of rocket. Here, manufacturing quality is demon-strated for the unguided rocket, and by simulating four typical manufacturing errors: erroneously manufactured warhead, misalignment between the warhead and engine chamber, asymmetrically installed propellant, and error in nozzle manufacturing. A new methodology is proposed, which combines a 3D CAD model of the asymmetrical projectile (due to manufacturing errors) and the improved Six-degrees-of-freedom (6DOF) model of its ﬂight into a comprehensive Monte-Carlo simulation. In that way, the rocket trajectory dispersion is correlated directly to the imperfection of the manufacturing process. Three quality levels are simulated (low, standard, and high quality), and each of the analyzed manufacturing errors depends on the chosen quality. The results show how important it is to impose the highest quality on nozzle manufacturing, and if this condition is not met, reveal if strict tolerances applied to other steps of the manufacturing process can compensate for the consequential drop of precision.


Introduction
Finding the right optimum between a criterion of precision and a criterion of production cost is of the utmost importance when designing a manufacturing process. The traditional approach to defining tolerances, based on empirical methods, does not give an adequate answer. For this reason, we propose a method that combines the 3D CAD (Computer-Aided Design) model of a realistic rocket (produced with some manufacturing errors) with the adjusted Six-degrees-of-freedom (6DOF) flight model. Two models are combined into one simulation that gives analytically justified manufacturing tolerances. The statistical analysis of the resulting rocket dispersion is a step further towards a deeper understanding of the imperfect projectile behavior, comparing to previous analyses [1,2]. The authors analyze a dispersion of rocket trajectories due to selected disturbances, among others due to the variations in inertia characteristics (and also launching pitch angle, propellant burning time, air density, etc.). In [1], a 122 mm rocket is analyzed, which was also chosen here as a case study. However, the cited papers introduce assumptions that narrow the applicability of their findings. In addition, limits of disturbances are set arbitrarily (for example, the longitudinal and lateral moments of inertia change by ±2%, rocket total mass and propellant mass change by ±1%, etc.).
Manufacturing errors were analyzed previously [3,4], leading to the conclusion that a certain level of uncertainty is inevitable when it comes to mechanical systems [5]. While [3] also analyzes the effect of dynamic asymmetry within arbitrarily chosen limits, a different approach is given in [4], where the serious consequences of dynamic imbalance are shown. However, a dynamic imbalance is not generated by the non-deterministic production imperfections, but intentionally by storing artillery shells on their side at hot conditions

Materials and Methods
A parametric 3D CAD model of an imperfect projectile is introduced. An example of the CAD model for such a projectile is presented in Figure 1. The 3D model is divided into two main components: the propellant that burns and therefore has changing inertia characteristics; and the platform that includes all other parts (fuse, warhead, engine chamber, nozzle, fins, etc.) with constant geometrical and inertia characteristics. Relevant inertia characteristics for the platform: mass m L , inertia tensor I L and center of gravity (CG) position ρ L = [x L y L z L ] T are obtained from a CAD program and depend on the simulated manufacturing errors for the particular rocket.

Manufacturing Errors Simulation
For the imperfect rocket, geometrical and inertia characteristics differ from that of an ideal projectile. Four typical manufacturing errors are simulated: erroneously manufactured warhead, misalignment between the warhead and engine chamber ("warhead-tochamber error"), asymmetrically installed propellant, and the thrust force misalignment caused by the error in nozzle manufacturing ( Figure 1).
To facilitate tracking of the imperfect rocket flight, a new coordinate system called geometrical coordinate system (G-CS) is introduced. The G-CS is linked to the engine chamber axes of symmetry. The x G axis coincides with the longitudinal axis of symmetry, and transversal axes y G and z G are fixed to rocket fins. The origin is placed at the center of gravity. The attitude of the G-CS is determined relative to the vehicle-carried coordinate system (O-CS, i.e., north-east-down c.s.) by Euler angles: φ G , θ G , and ψ G [11]. The transformation matrix from O-CS to G-CS is L GO = L x (φ G ) L y (θ G ) L z (ψ G ). Manufacturing errors are expressed in regard to G-CS. To facilitate tracking of the imperfect rocket flight, a new coordinate system called geometrical coordinate system (G-CS) is introduced. The G-CS is linked to the engine chamber axes of symmetry. The xG axis coincides with the longitudinal axis of symmetry, and transversal axes yG and zG are fixed to rocket fins. The origin is placed at the center of gravity. The attitude of the G-CS is determined relative to the vehicle-carried coordinate system (O-CS, i.e., north-east-down c.s.) by Euler angles: φG, θG, and ψG [11]. The transformation matrix from O-CS to G-CS is LGO = Lx (φG) Ly (θG) Lz (ψG). Manufacturing errors are expressed in regard to G-CS.

Warhead-to-Engine Chamber Misalignment
The warhead-to-engine chamber misalignment is described by two angles, as presented in Figure 2. The angle δHC lies in the plane of disturbance xG − xH, while the angle ϕHC defines the rotation of that plane relative to the reference plane xG − yG. Both angles δHC and ϕHC occur as stochastic parameters in the rocket 3D CAD model.

Warhead-to-Engine Chamber Misalignment
The warhead-to-engine chamber misalignment is described by two angles, as presented in Figure 2. The angle δ HC lies in the plane of disturbance x G − x H , while the angle ϕ HC defines the rotation of that plane relative to the reference plane x G − y G . Both angles δ HC and ϕ HC occur as stochastic parameters in the rocket 3D CAD model. To facilitate tracking of the imperfect rocket flight, a new coordinate system called geometrical coordinate system (G-CS) is introduced. The G-CS is linked to the engine chamber axes of symmetry. The xG axis coincides with the longitudinal axis of symmetry, and transversal axes yG and zG are fixed to rocket fins. The origin is placed at the center of gravity. The attitude of the G-CS is determined relative to the vehicle-carried coordinate system (O-CS, i.e., north-east-down c.s.) by Euler angles: φG, θG, and ψG [11]. The transformation matrix from O-CS to G-CS is LGO = Lx (φG) Ly (θG) Lz (ψG). Manufacturing errors are expressed in regard to G-CS.

Warhead-to-Engine Chamber Misalignment
The warhead-to-engine chamber misalignment is described by two angles, as presented in Figure 2. The angle δHC lies in the plane of disturbance xG − xH, while the angle ϕHC defines the rotation of that plane relative to the reference plane xG − yG. Both angles δHC and ϕHC occur as stochastic parameters in the rocket 3D CAD model. Angles are defined as follows: • δHC-the angle of misalignment in the plane of disturbance, dispersed according to the normal distribution (parameters of dispersion defined by the manufacturing quality level); • ϕHC-the radial angle, defining the radial position of angle δHC, 0 ≤ ≤ 360 dispersed according to the uniform distribution. Angles are defined as follows: • δ HC -the angle of misalignment in the plane of disturbance, dispersed according to the normal distribution (parameters of dispersion defined by the manufacturing quality level); • ϕ HC -the radial angle, defining the radial position of angle δ HC , 0 o ≤ ϕ HC ≤ 360 o dispersed according to the uniform distribution.

Warhead and Propellant Errors
The warhead error refers to the misalignment of the inner and outer warhead surface's axes of symmetry. Consequently, the explosive is not aligned with the warhead outer geometry (Figure 3a) [12]. The propellant installment error causes the propellant and chamber axes of symmetry not to coincide (Figure 3b). In analogy with the warhead-to-chamber error, these two errors are also defined by the corresponding pair of angles: δH and ϕH for the warhead error, and by δP and ϕP for the propellant error. Both δH and δP angles are dispersed according to the normal distribution, but not necessarily with the same standard deviations, as will be explained later in the case study analysis. The radial angles ϕH and ϕP are dispersed according to the uniform distribution, in a full circle 0-360°.
If there is a δHC angle, both the warhead and the explosive charge rotate around the warhead-chamber joint. Additionally, if there is a δH angle, the explosive charge rotates inside the warhead around the same point. As a result, there is a rotation of both parts relative to the rocket CG (CR in Figure 2) and their CGs are no longer on the xG axis. This was introduced into the 6DOF model via the modified inertia tensor (obtained from the 3D CAD model) and altered moment of aerodynamic force.

Thrust Force Misalignment
It is already known that, for unguided rockets, the thrust force misalignment results in high impact point dispersion and therefore must be kept within very strict boundaries [13]. Misalignment of the thrust force may occur due to the asymmetrical combustion, but also if the nozzle is not perfectly aligned with the xG axis of symmetry. Here the combustion is assumed to be perfect, so the thrust force misalignment occurs only due to the nozzle manufacturing error.
The nozzle is a rotationally symmetrical body. A new auxiliary coordinate system is introduced, called the nozzle coordinate system (N-CS); it is aligned with the nozzle longitudinal axis of symmetry xN. The plane xG − xN is called the nozzle disturbance plane. The origin of N-CS is located at the nozzle center with a radial position = [ 0 0] . In analogy with the warhead-to-chamber error, the nozzle error is defined over two angles: In analogy with the warhead-to-chamber error, these two errors are also defined by the corresponding pair of angles: δ H and ϕ H for the warhead error, and by δ P and ϕ P for the propellant error. Both δ H and δ P angles are dispersed according to the normal distribution, but not necessarily with the same standard deviations, as will be explained later in the case study analysis. The radial angles ϕ H and ϕ P are dispersed according to the uniform distribution, in a full circle 0-360 • .
If there is a δ HC angle, both the warhead and the explosive charge rotate around the warhead-chamber joint. Additionally, if there is a δ H angle, the explosive charge rotates inside the warhead around the same point. As a result, there is a rotation of both parts relative to the rocket CG (CR in Figure 2) and their CGs are no longer on the x G axis. This was introduced into the 6DOF model via the modified inertia tensor (obtained from the 3D CAD model) and altered moment of aerodynamic force.

Thrust Force Misalignment
It is already known that, for unguided rockets, the thrust force misalignment results in high impact point dispersion and therefore must be kept within very strict boundaries [13]. Misalignment of the thrust force may occur due to the asymmetrical combustion, but also if the nozzle is not perfectly aligned with the x G axis of symmetry. Here the combustion is assumed to be perfect, so the thrust force misalignment occurs only due to the nozzle manufacturing error.
The nozzle is a rotationally symmetrical body. A new auxiliary coordinate system is introduced, called the nozzle coordinate system (N-CS); it is aligned with the nozzle longitudinal axis of symmetry x N . The plane x G − x N is called the nozzle disturbance plane. The origin of N-CS is located at the nozzle center with a radial position ρ G N = [l N 0 0] T . In analogy with the warhead-to-chamber error, the nozzle error is defined over two angles: δ N and ϕ N , as in Figure 4. The vector of the thrust force → F T coincides with the x N axis and is deflected from the x G by the angle δ N .  Components of the thrust force in the G-CS are: and its moment in the geometrical coordinate system is where FT is the thrust force intensity, variable during the propellant combustion.
In the earlier industry-related experience, the nozzle error has been unsuccessfully corrected by imposing the highest level of manufacturing quality to other steps of production. Even if the danger of erroneous nozzle manufacturing is known, it is not known to what extent can its detrimental effect be compensated. As a method, the Monte Carlo simulation is chosen since it allows examining not only the impact of individual manufacturing errors but also their correlation that might otherwise be missed. Furthermore, it allows the analysis of the effect from a very wide spectrum of parameters (wind, mass, launching angle, etc.) that are not covered during the initial analysis, since it is focused on the production quality.

Monte-Carlo Simulation
The analysis of the effect that four manufacturing errors have on the impact points dispersion is performed by the Monte Carlo simulation. The value of Monte Carlo simulation for determining the manufacturing parameters is extensively described, for example in [14] where it is presented as a beneficial tool for analysis of tolerances used in mechanical assemblies. Its main advantage is flexibility and the ability to use various nonnormal input or output distributions. In that way, it can describe a wider spectrum of problems and input variables compared to some other analytical methods, for example, the Failure Modes and Effects Analysis (FMEA) method, although this one also results in preventively identifying problems and gives their gradation based on severity, probability of occurrence, and probability of detection [15]. The scheme of simulation that is used for this case study (previously presented in [16]) is shown in Figure 5: Components of the thrust force in the G-CS are: and its moment in the geometrical coordinate system is where F T is the thrust force intensity, variable during the propellant combustion.
In the earlier industry-related experience, the nozzle error has been unsuccessfully corrected by imposing the highest level of manufacturing quality to other steps of production. Even if the danger of erroneous nozzle manufacturing is known, it is not known to what extent can its detrimental effect be compensated. As a method, the Monte Carlo simulation is chosen since it allows examining not only the impact of individual manufacturing errors but also their correlation that might otherwise be missed. Furthermore, it allows the analysis of the effect from a very wide spectrum of parameters (wind, mass, launching angle, etc.) that are not covered during the initial analysis, since it is focused on the production quality.

Monte-Carlo Simulation
The analysis of the effect that four manufacturing errors have on the impact points dispersion is performed by the Monte Carlo simulation. The value of Monte Carlo simulation for determining the manufacturing parameters is extensively described, for example in [14] where it is presented as a beneficial tool for analysis of tolerances used in mechanical assemblies. Its main advantage is flexibility and the ability to use various non-normal input or output distributions. In that way, it can describe a wider spectrum of problems and input variables compared to some other analytical methods, for example, the Failure Modes and Effects Analysis (FMEA) method, although this one also results in preventively identifying problems and gives their gradation based on severity, probability of occurrence, and probability of detection [15]. The scheme of simulation that is used for this case study (previously presented in [16]) is shown in Figure 5:

Statistical Process Control, Quality, and Process Capability Indices
The basic concept of statistical process control relates to the comparison of obtained data and calculated control within specification limits. Statistical process control is recognized as a modern method for the analysis of process capabilities. The process capability assessment measures efficiency and effectiveness of the process when there are no special causes of variation, in other words in the state of statistical control. When the process is under control, it will be less likely for parameters of the observed process to go beyond the control limits.
If there is a need to determine whether the process is within the control limits, three statistical instruments are suggested [17,18]: control charts, histograms, and mathematical analysis of distribution. Before assessing the capability of the process, it is necessary to choose a critical parameter or variable that will be controlled for. More about the control charts parameters, processes that are either in control or out of control, type of quality control policy (including the frequency and duration of sampling), and the cost of a quality program can be found in [19]. Due to the inability to obtain a larger sample of data from the actual manufacturing system, some key assumptions are introduced. Using literature sources [20] and available data from manufacturers, assumptions of the quality levels were made, and the quality of manufacturing is divided into three groups: low-

Statistical Process Control, Quality, and Process Capability Indices
The basic concept of statistical process control relates to the comparison of obtained data and calculated control within specification limits. Statistical process control is recognized as a modern method for the analysis of process capabilities. The process capability assessment measures efficiency and effectiveness of the process when there are no special causes of variation, in other words in the state of statistical control. When the process is under control, it will be less likely for parameters of the observed process to go beyond the control limits.
If there is a need to determine whether the process is within the control limits, three statistical instruments are suggested [17,18]: control charts, histograms, and mathematical analysis of distribution. Before assessing the capability of the process, it is necessary to choose a critical parameter or variable that will be controlled for. More about the control charts parameters, processes that are either in control or out of control, type of quality control policy (including the frequency and duration of sampling), and the cost of a quality program can be found in [19]. Due to the inability to obtain a larger sample of data from the actual manufacturing system, some key assumptions are introduced. Using literature sources [20] and available data from manufacturers, assumptions of the quality levels were made, and the quality of manufacturing is divided into three groups: low-quality, standard-quality, and a high-quality level, according to the amount of variability in the manufacturing process.
The calculation is made on the basis of the potential process capability index Cp, which estimates the ability of the process to produce the output (good parts) in the case when the mean is centered between the specification limits. The general assumption is that the manufacturing process is perfectly centered: to achieve more realistic conditions, the probability of the mean shift will be included in the simulation model. Table 1 shows the potential process capability and related manufacturing quality. DPMO is defect per million opportunities and Z-value is the value of how many standard deviations fit inside the specification limits. More about the quality analysis and the need to quantify the influential parameters for their later assessment can be found in [21]. Since the index assumes the normal or near-normal data, the link between Cp and process variation σ can be made as follows: C p = (USL − LSL)/6σ where USL is the upper specification limit, and LSL is the lower specification limit. The USL and LSL values are assumed according to the literature (for example [3]) and are unique for all types of manufacturing errors.
Using chosen USL, LSL, and the quality level, the estimation of the process deviation is calculated as: More on the assessment of six-sigma quality indices can be found in [22].

Simulating Manufacturing Errors
Manufacturing errors are simulated in the rocket 3D CAD model. For each error, the appropriate pair of angles (nondeterministic variables) is introduced: δ which defines the angle of error or misalignment, and radial angles ϕ as explained in Section 2 and presented in  Radial angles ϕ are dispersed according to the uniform probability distribution: while δ angles are dispersed according to the normal probability density function: The parameter σ for each error and each quality level is determined as in Equation (3), according to the appropriate USL and LSL. Values are given in Table 2.
As can be seen from Table 2, the upper and lower specification limits (USL and LSL) for warhead error are higher than USL and LSL for other error limits, as preliminary analysis has proven that the effect of this error is the weakest. As opposed to that, USL and LSL for the nozzle error are lower compared to the limits for other errors. The probability distribution function (PDF) is defined using the values listed in Table 2. For example, for the warhead-to-chamber error (and the low manufacturing quality) specification limits are set as follows: The standard deviation for δ HC dispersion is then: σ HC = 0.66/6·0.8335 = 0.132, since Cp = 0.8335 is chosen as the average for the low-quality manufacturing level. Finally, if µ = 0 and σ HC = 0.132, the PDF for δ HC becomes: The same procedure is done for the other two quality levels, resulting in the same form of PDFs but using the parameter σ HC = 0.094 for the standard quality and σ HC = 0.060 for the high quality level. PDFs for δ HC are presented in Figure 6.
As can be seen from Table 2, the upper and lower specification limits (USL and LSL) for warhead error are higher than USL and LSL for other error limits, as preliminary analysis has proven that the effect of this error is the weakest. As opposed to that, USL and LSL for the nozzle error are lower compared to the limits for other errors. The probability distribution function (PDF) is defined using the values listed in Table 2. For example, for the warhead-to-chamber error (and the low manufacturing quality) specification limits are set as follows: USL = 0.33°, LSL = −0.33°.
The standard deviation for δHC dispersion is then: The same procedure is done for the other two quality levels, resulting in the same form of PDFs but using the parameter σHC = 0.094 for the standard quality and σHC = 0.060 for the high quality level. PDFs for δHC are presented in Figure 6. The same calculation is performed (but with different specification limits, as presented in Table 2) for δH, δP, and δN angles, giving the probability density function for each error and each of three quality level.

Flight Model Adapted for Imperfect Projectile
In classical flight mechanics [11], equations of motion are written in a frame coordinate system, the axes of which coincide with principal axes of inertia, or in its nonspinning variant, the ballistic coordinate system (B-CS). However, neither of these coordinate systems is suitable to describe the flight of the rocket with asymmetrical geometry or asymmetrical mass distribution. Therefore, the flight model is developed in the newly introduced geometric coordinate system (G-CS) and will be called G6DOF for distinction. The same calculation is performed (but with different specification limits, as presented in Table 2) for δ H , δ P , and δ N angles, giving the probability density function for each error and each of three quality level.

Flight Model Adapted for Imperfect Projectile
In classical flight mechanics [11], equations of motion are written in a frame coordinate system, the axes of which coincide with principal axes of inertia, or in its non-spinning variant, the ballistic coordinate system (B-CS). However, neither of these coordinate systems is suitable to describe the flight of the rocket with asymmetrical geometry or asymmetrical mass distribution. Therefore, the flight model is developed in the newly introduced geometric coordinate system (G-CS) and will be called G6DOF for distinction.

Governing Equations
Differential equations for G6DOF model are obtained from: • Extension of Newton's laws of motion for a body with variable mass. Includes the aerodynamic force F G A , thrust force F G T , gravity mg O and Coriolis force ma G c .
. V G K is a derivation of flight velocity, in G-CS: • Derivation of angular momentum where R is given as • Coordinates in the geocentric coordinate system (E-CS, Earth-fixed spherical coordinate system): .
with R being the radius of Earth and Ω E being its angular velocity.
The state vector X G for a rocket with mass asymmetry is defined according to Equations (7)-(11): A difference between the classical 6DOF and the new G6DOF flight model is visible, because in X G , components of H G are introduced instead of Ω G G . The aerodynamic velocity V G = [u v w] T is also defined in the G-CS. Total angle of attack is σ = arctan √ v 2 +w 2 u , and angle between the airflow plane and vertical plane around x G is φ = arctan(v/w). Derivations of σ and φ are: Wind is given only in horizontal plane, with no vertical wind: Angular velocity of the rocket is determined from the angular momentum: Ω G G = inv I G ·H G . The angular velocity of O-CS relative to the Earth is: with the following transformation matrix from E-CS to O-CS: L OE = L y 3π 2 + φ E · L z (λ). Projectile's angular velocity about axes of the airflow coordinate system (S-CS) can be deter- Here φ S is the attitude of the airflow around the longitudinal axis.
The presented model can be used for mid-range rockets flying at high altitudes. A change in gravity with altitude [13] could also be included into the system of equation.

Inertia Model
The inertia characteristics for propellant and platform can be obtained from the CAD model. A large number of assumptions that are common in the previously cited papers have been omitted, and the calculation is performed using the actual inertia characteristics. For example, the assumption is dismissed that the projectile inertia characteristics remain axis-symmetrical during the flight (in our method: I y = I z ), which corresponds better to the real case. In addition, products of inertia are no longer assumed to be equal to zero (in our method: I xy = I yz = I zx = 0). Even when the limits of production errors are chosen conservatively, a small change of inertia characteristics causes the rocket to deviate significantly from the target. It is no longer assumed that the main axes of inertia coincide with the outer geometry axes of symmetry, nor that the center of gravity (CR in Figure 2) remains on the x G -axis.
The analyzed rocket has the propellant shaped as a hollow cylinder, and the combustion occurring simultaneously over the external and internal surface. The propellant mass decreases linearly: m P (t) = m P0 (1 − t/t a ), where m P0 is the initial propellant mass and t a is the total combustion time. Platform mass m L and its C.G. position → ρ L are known. Because of dynamic imbalance, centers of gravity for the platform and the complete rocket are not on the axis of symmetry (Figure 2). Moments of inertia for the propellant change proportionally with its mass: I Px (t) = (m P (t)/m P0 )I Px0 . Inertia tensor of propellant is given relative to its principal axes and center of gravity (denoted CP in Figure 2): This inertia tensor must be recalculated relative to the entire rocket center of gravity CR, and relative to G-CS axes. Parallel axes theorem applies, and the transformation matrix L GP must be used: In analogy, the platform's inertia tensor also must be recalculated for the G-CS and translated to the projectile's center of gravity (giving I G L ). If platform axes of symmetry are not parallel to G-CS, then the appropriate transformation matrix should be used. Finally, tensor of inertia for the entire rocket is: I G R = I G P + I G L .

Aerodynamics of an Asymmetrical Rocket
To complete the G6DOF model, components of aerodynamic force and aerodynamic moment are required. If harmonics can be neglected, then aerodynamic coefficients are independent of the roll angle φ (the angle between the airflow plane and the x G -z G plane): with air density ρ, aerodynamic velocity V, and projectile reference area S. Components of aerodynamic moment for reference point RP in G-CS are: If external surface of the rocket is rotationally symmetrical, the following relations among aerodynamic gradients exist: α , C nr = C mq . In Equation (14), p * , q * , r * are non-dimensional components of the rocket angular velocity. In the absence of data for the Magnus effect, and considering that these coefficients are not relevant for this analysis, the following is assumed: C yα = 0, C zβ = 0, C mβ = 0, and C nα = 0. For an ideal rocket, components of aerodynamic force along G-CS in Equation (18) are: and for the reference point RP (defined with the position vector ρ RP ), the aerodynamic moment in Equation (19) is given as The aerodynamic moment in projectile's center of gravity (defined with position vector The warhead-to-chamber misalignment creates the additional aerodynamic normal force N H in the plane of disturbance, and the moment M H , as shown in Figure 7.
If external surface of the rocket is rotationally symmetrical, the following relations among aerodynamic gradients exist: In Equation (14), * , * , * are non-dimensional components of the rocket angular velocity. In the absence of data for the Magnus effect, and considering that these coefficients are not relevant for this analysis, the following is assumed: = 0, = 0, = 0, and = 0. For an ideal rocket, components of aerodynamic force along G-CS in Equation (18) are: and for the reference point RP (defined with the position vector ), the aerodynamic moment in Equation (19) is given as The warhead-to-chamber misalignment creates the additional aerodynamic normal force NH in the plane of disturbance, and the moment MH, as shown in Figure 7. The intensity of additional normal force is = ( /2) . The distance from the rocket apex to the point of application of the force NH is xHn, and thus its moment in the projectile center of gravity is = ( − ) where xR is the distance from the apex of the rocket to its center of gravity. The aerodynamic gradient and the position of the force application xHn are determined according to [23]. Since both NH and MH are tied to the warhead geometry, they should be transformed into the G-CS: The intensity of additional normal force is N H = ρV 2 S/2 C Nδ δ HC . The distance from the rocket apex to the point of application of the force N H is x Hn , and thus its moment in the projectile center of gravity is M H = N H (x R − x Hn ) where x R is the distance from the apex of the rocket to its center of gravity. The aerodynamic gradient C Nδ and the position of the force application x Hn are determined according to [23]. Since both N H and M H are tied to the warhead geometry, they should be transformed into the G-CS: Besides additional normal force N H , the warhead-to-chamber manufacturing error causes the additional induced drag. The aerodynamic derivative is C xδ and it is dependent on the total angle of attack σ, as well as on the angle of misalignment δ HC , as presented in Figure 8. Besides additional normal force NH, the warhead-to-chamber manufacturing error causes the additional induced drag. The aerodynamic derivative is Cxδ and it is dependent on the total angle of attack σ, as well as on the angle of misalignment δHC, as presented in Figure 8.
The solution is to replace the effect of these two angles with the resultant angle: where ∆α and ∆β are components of δHC angle: Δ = The aerodynamic coefficient of the total drag, incorporating the influence of the misalignment angle δHC, would be usually defined as in Equation (18)

Results and Discussion
To demonstrate the model, an unguided rocket similar to the 9M22U 122 mm projectile for the BM-21 GRAD multiple rocket launcher is analyzed as a case study (Figure 1). The rocket has the following characteristics: Given data are valid for an ideal rocket, while the non-ideal rocket (i.e., being erroneously manufactured) will have altered inertia characteristics. For such a rocket, the assumption Iy0 = Iz0 is no longer valid.
The inertia and motor characteristics are taken from technical data sheets, while aerodynamic gradients are determined by the method described in [24] and confirmed in [1]. The coefficient Cx0(Ma) from Equation (20) is determined according to the etalon CD58(Ma) and adjusted so that the ideal rocket has the same maximum range as stated in the Firing Tables (in the NATO standard atmosphere). For the coefficient of induced drag, it is assumed Cxσ = Czα.
The target is at 13700 m (which presents 2/3 of the maximum range). The rocket flight is simulated by the 6DOF model as defined by Equations (8)- (11), while geometric and The solution is to replace the effect of these two angles with the resultant angle: σ Hn = (α + ∆α) 2 + (β + ∆β) 2 (24) where ∆α and ∆β are components of δ HC angle: ∆α = δ HC ·sinϕ HC , ∆β = δ HC ·cosϕ HC . The aerodynamic coefficient of the total drag, incorporating the influence of the misalignment angle δ HC , would be usually defined as in Equation (18) but now with the application of σ Hn : C x (Ma) = C x0 (Ma) + C xσ 2 (Ma)σ Hn 2 .

Results and Discussion
To demonstrate the model, an unguided rocket similar to the 9M22U 122 mm projectile for the BM-21 GRAD multiple rocket launcher is analyzed as a case study (Figure 1). The rocket has the following characteristics: Given data are valid for an ideal rocket, while the non-ideal rocket (i.e., being erroneously manufactured) will have altered inertia characteristics. For such a rocket, the assumption I y0 = I z0 is no longer valid.
The inertia and motor characteristics are taken from technical data sheets, while aerodynamic gradients are determined by the method described in [24] and confirmed in [1]. The coefficient C x0 (Ma) from Equation (20) is determined according to the etalon C D58 (Ma) and adjusted so that the ideal rocket has the same maximum range as stated in the Firing Tables (in the NATO standard atmosphere). For the coefficient of induced drag, it is assumed C xσ = C zα .
The target is at 13,700 m (which presents 2/3 of the maximum range). The rocket flight is simulated by the 6DOF model as defined by Equations (8)- (11), while geometric and inertia characteristics for each (non-ideal) rocket are obtained from the 3D CAD model. Figure 9 shows some parameters of the motion for an ideal projectile. inertia characteristics for each (non-ideal) rocket are obtained from the 3D CAD model. Figure 9 shows some parameters of the motion for an ideal projectile. The flight time is 34 seconds, and the trajectory apex is at 1367 m. The maximum velocity is 702 m/s, and the final velocity is 307 m/s. These data comply well with data gathered from the Firing Tables and the available technical documentation.
Because of simulated manufacturing errors, impact points are dispersed around the target. Simulations are used here as an approximation for expensive polygon tests, which is a common practice today due to financial limitations [25].
Trajectories have been simulated for 1000 rockets (for each quality level) manufactured with non-deterministic warhead-to-chamber errors, and non-deterministic nozzle misalignment errors. Parameters for errors dispersion are explained in previous sections, with 99.7% probability: • error angles δH, δHC, δP and δN dispersed according to the normal distribution, with standard deviations depending on the production quality as given in Table 2; • 0 ≤ ≤ 360 for all radial angles, dispersed according to the uniform distribution. Figure 10 shows results of Monte Carlo simulations for three different quality levels. The flight time is 34 s, and the trajectory apex is at 1367 m. The maximum velocity is 702 m/s, and the final velocity is 307 m/s. These data comply well with data gathered from the Firing Tables and the available technical documentation.
Because of simulated manufacturing errors, impact points are dispersed around the target. Simulations are used here as an approximation for expensive polygon tests, which is a common practice today due to financial limitations [25].
Trajectories have been simulated for 1000 rockets (for each quality level) manufactured with non-deterministic warhead-to-chamber errors, and non-deterministic nozzle misalignment errors. Parameters for errors dispersion are explained in previous sections, with 99.7% probability: • error angles δ H , δ HC , δ P and δ N dispersed according to the normal distribution, with standard deviations depending on the production quality as given in Table 2; • 0 o ≤ ϕ i ≤ 360 o for all radial angles, dispersed according to the uniform distribution. Figure 10 shows results of Monte Carlo simulations for three different quality levels.  It is interesting to compare parameters of motion around the CoG for an ideal rocket, with a rocket produced with considerable errors (Case No. 1, a whole rocket made in low quality, production errors: δN = 0.151°, δH = 0.056°, δHC = 0.173°, δP = 0.091°). As presented in Figure 11, the angle of attack α and the angle of sideslip β for the non-ideal rocket ( Figure 11b) are larger in comparison to the results for the ideal rocket (Figure 11a). The stability of the non-ideal rocket is not endangered, not even if produced with considerable production errors as the one analyzed here. In Figure 11c, the angle of attack for both projectiles is directly compared for the first five seconds, since later the angle of attack diminishes, and differences are difficult to be graphically presented. It is interesting to compare parameters of motion around the CoG for an ideal rocket, with a rocket produced with considerable errors (Case No. 1, a whole rocket made in low quality, production errors: δ N = 0.151 • , δ H = 0.056 • , δ HC = 0.173 • , δ P = 0.091 • ). As presented in Figure 11, the angle of attack α and the angle of sideslip β for the non-ideal rocket (Figure 11b) are larger in comparison to the results for the ideal rocket (Figure 11a). The stability of the non-ideal rocket is not endangered, not even if produced with considerable production errors as the one analyzed here. In Figure 11c, the angle of attack for both projectiles is directly compared for the first five seconds, since later the angle of attack diminishes, and differences are difficult to be graphically presented. The drift and range results are independent of each other, which is proved by calculation of variance, covariance, and correlation coefficient consequently. The sample correlation coefficients are negligibly small, ranging from 0.01 to 0.06. Furthermore, using the test of significance, with a significance level of 5%, there is insufficient evidence to conclude that there is a correlation of the given data. As for normality of the response data (range and drift), the test was performed before the prediction interval calculation. The following is concluded: • Descriptive location parameters (mean, mode and median) coincide. • The absolute values of skewness and kurtosis for the analyzed variables are in the range of 0-0.1 and 0-0.5, respectively. • In addition, the assumption of normality was tested using the normal-probability plots and Shapiro Wilk test.
Since the above-mentioned facts point out that the variables are normally distributed, the prediction interval is appropriate for describing the single-point estimation.
As Figure 10 and Table 3 show, there is a significant reduction of impact point dispersion area when higher quality levels are applied. The area of impact points dispersion corresponds well with the Firing Tables for GRAD rocket, used here as a corrective element since it obtains data from previously conducted extensive experiments and test firings. Results in Table 3 show 95% prediction intervals for impact points by the range and drift, and the impact point dispersion area for each of three simulated quality The drift and range results are independent of each other, which is proved by calculation of variance, covariance, and correlation coefficient consequently. The sample correlation coefficients are negligibly small, ranging from 0.01 to 0.06. Furthermore, using the test of significance, with a significance level of 5%, there is insufficient evidence to conclude that there is a correlation of the given data. As for normality of the response data (range and drift), the test was performed before the prediction interval calculation. The following is concluded:

•
Descriptive location parameters (mean, mode and median) coincide.

•
The absolute values of skewness and kurtosis for the analyzed variables are in the range of 0-0.1 and 0-0.5, respectively. • In addition, the assumption of normality was tested using the normal-probability plots and Shapiro Wilk test.
Since the above-mentioned facts point out that the variables are normally distributed, the prediction interval is appropriate for describing the single-point estimation.
As Figure 10 and Table 3 show, there is a significant reduction of impact point dispersion area when higher quality levels are applied. The area of impact points dispersion corresponds well with the Firing Tables for GRAD rocket, used here as a corrective element since it obtains data from previously conducted extensive experiments and test firings. Results in Table 3 show 95% prediction intervals for impact points by the range and drift, and the impact point dispersion area for each of three simulated quality levels. Calculation of 95% prediction intervals, as well as areas of the impact point dispersion, shows that there is a significant reduction of 51% and 79% when standard or high quality level are applied, compared to that with the low quality of manufacturing. Knowing that the nozzle error, (i.e., the consequential thrust force misalignment) is causing a significant impact point dispersion, the analysis was expanded with three more cases, again using a Monte Carlo simulation with 1000 iterations for each case. Three new simulations examine the overall projectile dispersion when the nozzle error is manufactured at the low level, as in the already presented Case No. 1, but now with other rocket components manufactured at the standard level (Case No. 4), at the high level (Case No. 5) and even perfectly (Case No. 6). Table 4 gives dispersion characteristics for all six cases (No. 1-No. 6), given by range and drift.
Finally, three more cases (numbered 7-9) are analyzed to inspect what happens if the nozzle error is removed as a contributing factor (δ N = 0). An appropriate Monte Carlo simulation is prepared with 1000 iterations for each of three cases. Results show that in such a case, the dispersion area multifold decreases, regardless of whether other rocket parts are manufactured at a standard (Case No. 8) or even low quality (Case No. 7). If all parts are manufactured at high quality and δ N = 0 (Case No. 9), a dispersion area matches the one of a guided rocket [26].
As expected, a dispersion of impact points decreases for better manufacturing quality. However, cases 4-6 indicate that the erroneously manufactured nozzle prevents a significant decrease of dispersion area, even if all other steps of production are raised to the highest quality level. Raising other steps of production from low to the standard quality level results in only 7% smaller dispersion area, and even raising it to the highest quality level gives a decrease of only 11%. In other words, if the nozzle is manufactured erroneously, all efforts and additional costs connected to the higher quality level are useless. This is valid even if a maximum nozzle error is only 2 mrad.
In Figure 12, two charts are presented, showing how critical nozzle manufacturing is. Both charts give the impact point dispersion if the entire rocket, except for the nozzle, is manufactured at the standard quality level. On the first chart, the nozzle is manufactured with low quality, and on the second, the nozzle error = 0. A difference between the two charts is a direct effect of nozzle imperfection, and it clearly causes the largest portion of the overall rocket dispersion. A preliminary DOE was conducted to determine the effects of different errors (factors). The significance of the factors was obtained using the ANOVA method, and standardized effects (presented via coefficient estimation) were calculated as presented in Tables 5 and 6: The results confirm the conclusion that nozzle error has the highest effect on range and drift. As expected, a dispersion of impact points decreases for better manufacturing quality. However, cases 4-6 indicate that the erroneously manufactured nozzle prevents a significant decrease of dispersion area, even if all other steps of production are raised to the highest quality level. Raising other steps of production from low to the standard quality level results in only 7% smaller dispersion area, and even raising it to the highest quality level gives a decrease of only 11%. In other words, if the nozzle is manufactured erroneously, all efforts and additional costs connected to the higher quality level are useless. This is valid even if a maximum nozzle error is only 2 mrad.
In Figure 12, two charts are presented, showing how critical nozzle manufacturing is. Both charts give the impact point dispersion if the entire rocket, except for the nozzle, is manufactured at the standard quality level. On the first chart, the nozzle is manufactured with low quality, and on the second, the nozzle error = 0. A difference between the two charts is a direct effect of nozzle imperfection, and it clearly causes the largest portion of the overall rocket dispersion.
(a) (b) Figure 12. Impact points dispersion with (a) and without nozzle error (b). Figure 12. Impact points dispersion with (a) and without nozzle error (b).

Conclusions
A method for estimating the effect of manufacturing errors on rocket precision is presented. The severity of manufacturing errors is examined by a Monte Carlo simulation, with 1000 iterations for each of the nine cases analyzed. These cases mutually differ according to the quality level imposed on particular parts of the rocket.
The Monte Carlo simulation consists of the non-deterministic manufacturing errors generator, a 3D CAD model of the rocket which provides changes of the inertia characteristics, and an adjusted 6DOF model of rockets flight which gives the position of impact points, dispersed due to simulated manufacturing errors. A new coordinate system (geometric CS) is introduced, helping to develop an improved 6DOF flight model (G6DOF) that facilitates tracking of imperfect or asymmetric rocket flight.
Monte Carlo simulations (denoted as Case No.1-9) analyze different combinations of four simulated manufacturing errors and three manufacturing quality levels. Extensive analysis singled out the nozzle error as particularly critical, which is in line with the conclusions of previous papers. However, it is shown not only that this error is the dominant cause of dispersion, but also that its impact is so strong, that even raising the manufacturing quality of all other parts of the rocket to the highest quality level does not compensate for the detrimental effect of the nozzle error.
That means that it is possible to achieve a reduction in the rocket cost by using a standard or even low quality level for less critical stages of production. As for the nozzle manufacturing, the strictest tolerances must be imposed, because by increasing the potential process capability index Cp, the projectile dispersion area decreases by up to 79%.
Analytically selected manufacturing tolerances can therefore achieve a major improvement in rocket accuracy, with a minimal increase in production costs. Alternative solutions regularly fail to comply with at least one of the desired criteria: either they result in a too expensive projectile or give impact points dispersion that is too large. Future research of the price of quality could further improve the algorithm for setting analytically validated manufacturing tolerances.