Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 51, 1, pp. 235-246, Warsaw 2013 EQUATIONS OF MOTION OF A SPIN-STABILIZED PROJECTILE FOR FLIGHT STABILITY TESTING Leszek Baranowski Military University of Technology, Faculty of Mechatronics and Aerospace, Warsaw, Poland e-mail: leszek.baranowski@wat.edu.pl The paper presents amathematical model of motion of a balanced spin-stabilized projectile considered as a rigid bodywith 6 degrees of freedom.Themodeling uses coordinate systems conforming to Polish and International Standard ISO 1151. The design of kinematic equ- ations describingmotion around the center of mass uses the system of Tait-Bryan angles or Euler parameters.The total angle of attack and aerodynamic roll angle express aerodynamic forces andmoments. Key words: spin-stabilized projectile, flight stability, exterior ballistics, equations of motion of projectile 1. Introduction Now, ballistic computations usemathematical models of projectile motion of varying degrees of simplification dependingon thepurposeof themodel.Oneof the following twomodels is used for developing firing tables: the point-mass model (with 2 degrees of freedom) describingmotion of the center of projectile mass with an underlying assumption that the projectile becomes ideally stabilized on its trajectory and the effect of aerodynamic forces can be substituted with the effect of drag force or the modified point-mass model (with 4 degrees of freedom) with one of its implementation contained in STANAG 4355. Amodel representing the projectile as a rigid body is one of themost often used for testing dynamic properties of the projectile. In this model, aviation angles (Ψ,Θ,Φ) are used to de- termine angular position of the projectile relative to the ground-fixed system, and the angle of attack α and angle of sideslip β to determine angular position of the projectile relative to air flow. Flight stability testing requires a projectilemotionmodel inwhich the projectile is represen- ted as a rigid body with 6 degrees of freedom, addressing the effect of full aerodynamic force, including specifically Magnus force and moment, to enable stimulation of actual atmospheric flight, particularly for large quadrant elevation QE. This relates to the fact that the inclination angle of the projectile Θ often comes up to 90◦ in the final flight phase and the total angle of attack αt (contained between the projectile axis and the relative velocity vector) can become large (40 or more degrees) near the vertex. To develop such a mathematical model, the work uses standard coordinate systems confor- ming to Polish and International Standard ISO 1151, provided that the transformation matrix between the ground-fixed system Oxgygzg and the body-fixed system Oxyz uses the new sys- tem of Tait-Bryan angles (Θn,Ψn,Φn) instead of the conventional aviation angles (Ψ,Θ,Φ) for the avoidance of singularities in kinematic equations for projectile motion around the center of mass. In addition, the paper proposes kinematic equations of motion of the projectile as a rigid body based on Euler parameters. Toeliminate theerror fromcomputationof components of theaerodynamic forceandmoment in the case when the projectile axis deviates significantly from the relative velocity vector, the 236 L. Baranowski paperproposes expressing the aerodynamic forces andmomentswith the total angle of attack αt and the aerodynamic roll angle ϕ rather than with the conventional angles: angle of attack α and angle of sideslip β. 2. Using the new system of Tait-Bryan angles in designing kinematic equations for motion of the projectile as a rigid body To avoid singularities in kinematic equations of motion around the center of projectile mass in a simulation of firing at maximum quadrant elevations (where the inclination angle of the projectile Θ often comes up to 90◦ in the final flight phase), the transformationmatrix between the ground-fixed system Oxgygzg and the body-fixed system Oxyz was derived using the new rotation system (system of Tait-Bryan angles) shown in Fig. 1 (Roberson and Shwertassek, 1988; Wittenburg, 2008) instead of the conventional aviation angles (ISO 1151, 1988): azimuth angle Ψ, inclination angle Θ and bank angle Φ. The first rotation is around the horizontal axis of the ground-fixed system Oyg by the new angle of inclination Θn, the second rotation is around instantaneous axis Oz′g by the new angle of azimuth Ψn and the third rotation is around the axis Oz′′g by the new angle of bank Φn. The transformationmatrix between the ground-fixed system Oxgygzg and thebody-fixedsystem Oxyz using the new system ofTait-Bryan angles can be obtained from the following dependence LΦnΨnΘn =LΦnLΨnLΘn (2.1) Using the formulas for elementary matrices (Fig. 1) will provide the following LΦnΨnΘn = (2.2)    cosΘncosΨn sinΨn −sinΘncosΨn sinΘn sinΦn− cosΘn sinΨncosΦn cosΨncosΦn cosΘn sinΦn+sinΘn sinΨncosΦn sinΘncosΦn+cosΘn sinΨn sinΦn −cosΨn sinΦn cosΘncosΦn− sinΘn sinΨn sinΦn    The angular velocity of the body-fixed system Oxyz relative to the ground-fixed system (see Fig. 1) can be expressed with vectors of angular velocities of the new system of Tait-Bryan angles as Ω= Ψ̇n+Θ̇n+Φ̇n, and its components along the axis of the body-fixed system Oxyz can be expressed with the following dependence    p q r    =LΦnΨnΘn    0 Θ̇n 0    +LΦnΨn    0 0 Ψ̇n    +LΦn    Φ̇n 0 0    (2.3) Resolving appropriate matrix multiplications in equation (2.3) will provide the following    p q r    =    0 sinΨn 1 sinΦn cosΨncosΦn 0 cosΦn −cosΨn sinΦn 0       Ψ̇n Θ̇n Φ̇n    (2.4) Using the concept of inverse matrix, the equation for derivatives of transformation angles can be expressed as follows    Ψ̇n Θ̇n Φ̇n    =    0 sinΨn 1 sinΦn cosΨncosΦn 0 cosΦn −cosΨn sinΦn 0    −1   p q r    (2.5) Equations of motion of a spin-stabilized projectile... 237 Fig. 1. (a) Rotation of ground-fixed system Oxgygzg around axis Oyg by angle Θn, (b) rotation of instantaneous system Ox′gy ′ gz ′ g around axis Oz ′ g by angle Ψn, (c) rotation of instantaneous system Ox′′gy ′′ gz ′′ g around axis Ox ′′ g by angle Φn 3. Alternative method of determining projectile position relative to air flow for computation of aerodynamic forces and moments Determination of aerodynamic forces and moments affecting the projectile in flight requires computation of the angular position of the projectile relative to air flow (or projectile velocity vector with respect to the air V). The most popular method of computing aerodynamic force components consists of deter- mining the angle of attack α and angle of sideslip β. For axial-symmetric artillery projectiles however, it is not themost convenient one because large spin produces continuous change of the angles even if the angles that the forces andmoments really depend on do not change so fast. 238 L. Baranowski Aerodynamic forces affecting spinning projectiles operate in the plane of drag and perpen- dicularly to the plane of drag (Magnus force) independently of the projectile bank angle. Accordingly, for axial-symmetric flying objects, it is better (for the determination of the angular position of the projectile relative to the vector of velocity V) to use angles that are independent of the projectile spin, such as the total angle of attack αt and aerodynamic roll angle ϕ (Baranowski, 2006). The values of the angles, shown in Fig. 2, can be computed from the following equations αt =arctan √ (wK −wW)2+(vK −vW)2 uK −uW ϕ=arctan wK −wW vK −vW (3.1) Fig. 2. Illustration of spatial position of angles αt and ϕ The total aerodynamic force RA and total aerodynamic moment MAO acting on axial- symmetric spinning projectiles can be presented as follows (Fig. 3) R A =RAα +R A Ω M A O =M A Oα+M A OΩ (3.2) where RAα – aerodynamic force operating in the plane of drag, resulting from the effect of air on non-spinning projectile, the longitudinal axis of which is inclined from the air flow direction by the angle αt RAΩ – aerodynamic force acting perpendicularly to the plane of drag, resulting from the projectile spin and angle αt (Magnus force) MAOα – aerodynamic moment acting on a non-spinning projectile MAOΩ – aerodynamic moment resulting from the projectile spin and angle αt. To facilitate thedeterminationof components of theaerodynamic forceandmomentactingon the spinning projectile, the body-fixed system Oxyz uses the splitting of aerodynamic force RAα acting in the plane of drag into a component following the longitudinal axis of the projectile XA = CAXSρV 2/2 and a component perpendicular to the longitudinal axis of the projectile PA =CAN(M,αt)SρV 2/2 (M –Mach number), Fig. 3. The aerodynamic moment MAOα produced by the force P A is referred to as: • overturning moment, for spin-stabilized projectiles; • or stabilizing moment, for fin-stabilized projectiles. For artillery projectiles, it can be expressed with the coefficient of overturning moment Cm(M,αt) as follows MAOα =Cm(M,αt) ρV 2 2 Sl (3.3) Equations of motion of a spin-stabilized projectile... 239 Fig. 3. Components of aerodynamic force acting on the projectile in flight Using the angle ϕ, projections of the force PA and moment MAOα on the axes Oy and Oz of the body-fixed system Oxyz can be expressed as follows YAα =C A N(M,αt) ρV 2 2 S cosϕ ZAα =C A N(M,αt) ρV 2 2 S sinϕ MAα =Cm(M,αt) ρV 2 2 Slsinϕ NAα =−Cm(M,αt) ρV 2 2 Slcosϕ (3.4) Using notations conforming to ISO 1151 (1988), components of the total aerodynamic for- ce RA in the body-fixed system Oxyz take the following form: — axial force XA =−[CAX0(M)+C A Xα2(M)α 2 t ] ρV 2 2 S (3.5) — transverse force YA = [−CAZα(M)αt cosϕ+C A Ypα(M)p ∗αt sinϕ] ρV 2 2 S (3.6) — normal force ZA = [−CAZα(M)αt sinϕ−C A Ypα(M)p ∗αtcosϕ] ρV 2 2 S (3.7) In turn, components of the total aerodynamicmoment MAO in the body-fixed system Oxyz can be expressed as follows LA =CAlp(M) ρV 2 2 p∗Sl MA = [CAmα(M)αt sinϕ+C A mq(M)q ∗+CAnpα(M)p ∗αtcosϕ] ρV 2 2 Sl NA = [−CAmα(M)αtcosϕ+C A mq(M)r ∗+CAnpα(M)p ∗αt sinϕ] ρV 2 2 Sl (3.8) In accordance with equations (3.5)-(3.8), determination of themain aerodynamic properties of ground artillery projectiles consists of computing the following quantities: • axial force coefficient for αt =0 – CAX0(M) • derivative of the axial force coefficient – CA Xα2 (M) 240 L. Baranowski • derivative of the normal force coefficient – CAZα(M) • derivative of theMagnus force coefficient – CAYpα(M) • derivative of the spin dampingmoment coefficient – CAlp(M) • derivative of the overturning moment coefficient – CAmα(M) • derivative of the pitch dampingmoment coefficient – CAmq(M) • derivative of theMagnus moment coefficient – CAnpα(M) where CAXα2 = ∂2CAX ∂α2 CAZα = ∂CAZ ∂α · · · CAmα = ∂CAm ∂α CAlp = ∂CAl ∂p∗ CAmq = ∂CAm ∂q∗ · · · CAYpα = ∂2CAY ∂p∗α CAnpα = ∂2CAn ∂p∗α CAnp = ∂CAn ∂p∗ p∗ = pd 2V q∗ = qd 2V r∗ = rd 2V (3.9) In the case when the Magnus moment coefficient shows strong non-linear reliance on the total angle of attack, the following equation can be used CAnpα(M)p ∗αt =C A np(M,αt)p ∗ (3.10) 4. Mathematical model of motion of the projectile as a rigid body There are two groups ofmethods for the development ofmathematicalmodels ofmotion of flying objects based on the principles of classical and analytical mechanics. In the group ofmethods of analytical mechanics, one can distinguishmethods based on inertial generalized coordinates and referring directly to the Hamilton principle and Lagrange equations (Koruba et al., 2010) and the methods consisting in applying the equations of analytical mechanics in quasi-coordinates, e.g. Boltzman-Hamel equations (Ładyżyńska-Kozdraś and Koruba, 2012). Classical mechanics uses the law of change of themomentum and angular momentum of a rigid body (Gacek, 1997; Kowaleczko and Żyluk, 2009). Based on the principles of classicalmechanics, spatialmotion of the projectile as a rigid body in the frame moving together with the projectile, with the origin of coordinates located in the center of mass of the projectile, can be described with the following vector equations: — vector dynamic equations of motion m (δVK dt +Ωr×VK ) =RA+G+Fc δKO dt +Ωr×KO =MAO (4.1) — vector kinematic equations of motion drK dt =VK Ω= Ψ̇n+Θ̇n+ Φ̇n (4.2) where VK – vector of the projectile velocity with respect to the ground KO – vector of the projectile angular momentum relative to its center of mass Ω – vector of the projectile angular velocity Ωr – vector of the angular velocity of the frame moving together with the projectile respect to the ground-fixed system Oxgygzg ΩZ – vector of the angular velocity of the Earth Equations of motion of a spin-stabilized projectile... 241 rK = [xg,yg,zg] – initial position vector of the center of mass of the projectile and its components in the ground-fixed system Oxgygzg RA = [XA,YA,ZA] – vector of the aerodynamic force and its components in the body-fixed system Oxyz G= [Gxg,Gyg,Gzg] – vector of the gravity force and its components in the ground- fixed system Oxgygzg Fc =−2m(ΩZ ×VK) – vector of the Coriolis force due to rotation of the Earth MAO = [L A,MA,NA] – vector of the aerodynamicmoment and its components in the body-fixed system Oxyz. The scalar formof the foregoing vector – dynamic and kinematic – equations (in appropriate coordinate systems), together with the complementary equations, represents a mathematical model of motion of the projectile as a rigid body. 4.1. Scalar form of equations of motion in the body-fixed system Oxyz In its final vector-matrix form, the mathematical model of motion of ground artillery pro- jectiles as rigid bodies contains the following groups of equations: — dynamic differential equations of motion of the projectile center of mass in the body-fixed system Oxyz    u̇K v̇K ẇK    =    XA/m YA/m ZA/m    +LΦnΨnΘn    gxg +FCxg/m gyg +FCyg/m gzg +FCzg/m    +    0 r −q −r 0 p q −p 0       uK vK wK    (4.3) where components of the Coriolis force in the ground-fixed system Oxgygzg have the following form    FCxg FCyg FCzg    =    2Ω(cos(lat)sin(AZ)wKg− sin(lat)vKg) 2Ω(cos(lat)cos(AZ)wKg+sin(lat)uKg) −2Ω(cos(lat)cos(AZ)vKg+cos(lat)sin(AZ)uKg)    (4.4) for a spherical model of the Earth, components of the gravitational acceleration in the ground- fixed system Oxgygzg can be expressed as follows (STANAG 4355 Ed.3, 2009)    gxg gyg gzg    = gn    −xg/Rz −yg/Rz 1+2zg/Rz    (4.5) and Ω = 7.292115 · 10−5 rad/s – angular speed of the Earth, gn = 9.80665(1 − 0.0026cos(2 lat))m/s2 – magnitude of acceleration due to gravity at the mean sea level, Rz = 6356766m – radius of the sphere, locally approximating the geoid, lat – latitude of the launch point, for the southern hemisphere lat is negative [deg], AZ – azimuth (bearing) of the xg axis measured clockwise from true North [mil]. — kinematic differential equations of motion of the projectile center of mass    ẋg ẏg żg    =    uKg vKg wKg    =LTΦnΨnΘn    uK vK wK    (4.6) —dynamic differential equations of rotational motion about the projectile center ofmass in the body-fixed system Oxyz overlapping with the principle central axes of inertia    Ix 0 0 0 Iy 0 0 0 Iz       ṗ q̇ ṙ    =    LA MA NA    +    0 r −q −r 0 p q −p 0       Ix 0 0 0 Iy 0 0 0 Iz       p q r    (4.7) 242 L. Baranowski —kinematic differential equations of rotational motion about the projectile center of mass    Ψ̇n Θ̇n Φ̇n    =    0 sinΦn cosΦn 0 cosΦn/cosΨn −sinΦn/cosΨn 1 −cosΦn tanΨn sinΦn tanΨn       p q r    (4.8) — equation for the total angle of attack αt αt =        π 2 if uK −uW =0 arctan √ (wK −wW)2+(vK −vW)2 uK −uW otherwise (4.9) — equation for the aerodynamic roll angle ϕ sinϕ=        0 if √ (wK −wW)2+(vK −vW)2 wK −wW √ (wK −wW)2+(vK −vW)2 otherwise (4.10) and cosϕ=        1 if √ (wK −wW)2+(vK −vW)2 vK −vW √ (wK −wW)2+(vK −vW)2 otherwise (4.11) — complementary equations γ =arcsin wKg VK χ=arctan vKg uKg u=uK −uW v= vK −vW w=wK −wW V = √ u2+v2+w2 VK = √ u2Kg+v 2 Kg+w 2 Kg (4.12) where: u,v,w – components of the vector of projectile velocity with respect to the air V in the body-fixed system Oxyz uK,vK,wK – components of the vector of projectile velocity with respect to the ground VK in the body-fixed system Oxyz uKg,vKg,wKg – components of the vector of projectile velocity with respect to the ground VK in the ground-fixed system Oxgygzg uW ,vW ,wW – components of the vector of wind velocity with respect to the gro- und VK in the body-fixed system Oxyz γ,χ – path inclination angle and path azimuth angle, respectively. A comparison of components of the matrix LΦnΨnΘn, Eq. (2.2), with the matrix LΦΘΨ (Baranowski, 1998; Gacek, 1997) reveals relations between the aviation angles and the new Tait-Bryan angles sinΘ=sinΘncosΨn sinΨ = sinΨn cosΘ sinΦ= cosΘn sinΦn cosΘ (4.13) 4.2. Using Euler parameters in designing kinematic equations of motion for the projectile as a rigid body Also Euler parameters (Gajda, 1990; Roberson and Shwertassek, 1988) in the form of qu- aternions can be used for the avoidance of singularities in kinematic equations ofmotion around Equations of motion of a spin-stabilized projectile... 243 the center of mass. According to Eulers theorem, an object placed in 3D space can be moved from any starting position to any target position with a single rotation around a single axis. So, to define spatial orientation of amovable axis system (e.g. the body-fixed system Oxyz) relative to a fixed system (e.g. the ground-fixed system), it is enough to specify three direction cosines of the axis of rotation (for instance, using Tait-Bryan angles: αE, βE,γE as parameters describing the instantaneous position of the axis of rotation) and the value of the angle of rotation around the axis δE. These four numbers (3 direction cosines and the angle of rotation) are known as Euler parameters and can bewritten in the formof quaternions λ (Gajda, 1990; Gosiewski andOrtyl, 1995; Wittenburg, 2008). Quaternions are defined as vector quantities with 4 degrees of freedom λ=λ0+λ1i+λ2j+λ3k (4.14) where i, j, k as imaginary numbers meet the following conditions i2 = j2 = k2 =−1 ij=−ji= k jk=−kj= i ki=−ik= j and λ0 =cos δE 2 λ1 =cosαE sin δE 2 λ2 =cosβE sin δE 2 λ3 =cosγE sin δE 2 Quaternion components have to meet an additional combining equation (requirement for orthogonality) λ20+λ 2 1+λ 2 2+λ 2 3 =1 (4.15) The transformation matrix T between the body-fixed system Oxyz and the ground-fixed system Oxgygzg can be presented in two ways: — using aviation angles: Ψ,Θ, Φ (Baranowski, 1998; Gacek, 1997) T=L−1 ΦΘΨ =    cosΘcosΨ −cosΦsinΨ+sinΦsinΘcosΨ sinΦsinΨ+cosΦsinΘcosΨ cosΘsinΨ cosΦcosΨ+sinΦsinΘsinΨ −sinΦcosΨ+cosΦsinΘsinΨ −sinΘ sinΦcosΘ cosΦcosΘ    (4.16) — using quaternions (Gosiewski and Ortyl, 1995) T=2    (λ20+λ 2 1−λ 2 2−λ 2 3)/2 λ1λ2−λ0λ3 λ1λ3+λ0λ2 λ1λ2+λ0λ3 (λ20−λ 2 1+λ 2 2−λ 2 3)/2 λ2λ3−λ0λ1 λ1λ3−λ0λ2 λ2λ3+λ0λ1 (λ20−λ 2 1−λ 2 2+λ 2 3)/2    (4.17) Using quaternions, in the case of deriving equations of motion for the projectile as a rigid body in the body-fixed system Oxyz: — the kinematic differential equations of motion of the projectile center of mass (4.6) are as follows    ẋg ẏg żg    =    uKg vKg wKg    =T    uK vK wK    (4.18) 244 L. Baranowski — the kinematic differential equations of rotational motion about the projectile center of mass (4.8) have the following form (Baranowski et al., 2005; Baranowski, 2006) dλ0 dt = 1 2 (−λ1p−λ2q−λ3r) dλ1 dt = 1 2 (λ0p−λ3q+λ2r) dλ2 dt = 1 2 (λ3p+λ0q−λ1r) dλ3 dt = 1 2 (−λ2p+λ1q+λ0r) (4.19) Unlike the description using Euler and Tait-Bryan angles, it is a system of four differential equations in which the solution remains within the [−1,1] range, which facilitates numerical computations. The main computational problem in the quaternion model is the meeting of combining equation (4.15). Thequaternions are “improved” in order to satisfy this equation.The improving algorithm has the following form (Ortyl, 2000)      λ̇0 λ̇1 λ̇2 λ̇3      = 1 2      0 −p −q −r p 0 r −q q −r 0 p r q −p 0           λ0 λ1 λ2 λ3      −εw      λ0 λ1 λ2 λ3      (4.20) where εw – rate of violation of the combining equation (ideally εw =0) εw =λ 2 0+λ 2 1+λ 2 2+λ 2 3−1 (4.21) Theuse of quaternions relates todifficultywithphysical interpretation of quaternions as they relate to the orientation of the axis of rotation rather than the orientation of the object itself. Therefore, to interpret the computation results correctly, we need to transform these parameters onto aviation angles, which are natural coordinates defining the position of the flying object in space. Taking advantage of the fact that individual components of the matrix T are equal one to another, based on equations (4.16) and (4.17), the following relations can be established between the aviation angles and quaternions sinΘ=−T31 =2(λ0λ2−λ1λ3) − π 2 ¬Θ¬ π 2 tanΨ = T21 T11 = 2(λ1λ2+λ0λ3) λ20+λ 2 1−λ 2 2−λ 2 3 −π<Ψ ¬π tanΦ= T32 T33 = 2(λ2λ3+λ0λ1) λ20−λ 2 1−λ 2 2+λ 2 3 0<Φ¬ 2π (4.22) and between quaternions and aviation angles (Gajda, 1990) λ0 =cos Ψ 2 cos Θ 2 cos Φ 2 +sin Ψ 2 sin Θ 2 sin Φ 2 λ1 =cos Ψ 2 cos Θ 2 sin Φ 2 − sin Ψ 2 sin Θ 2 cos Φ 2 λ2 =cos Ψ 2 sin Θ 2 cos Φ 2 +sin Ψ 2 cos Θ 2 sin Φ 2 λ3 =sin Ψ 2 cos Θ 2 cos Φ 2 − cos Ψ 2 sin Θ 2 sin Φ 2 (4.23) Using quaternions instead ofTait-Bryan angles in kinematic equations ofmotion for artillery projectiles can provide the following benefits: • complete elimination of singular points in computing projectile position, Equations of motion of a spin-stabilized projectile... 245 • shorter computation time, • less numerical errors during simulation thanks to algebraic computations instead of nu- merical computations of trigonometric functions in the form of expansion in Taylor series with omission of terms of higher orders (which is the casewith solving kinematic equations of motion for projectile incorporating Tait-Bryan angles). 5. Summary and conclusions The paper presents a complete mathematical model of motion of a balanced spin-stabilized projectile considered as a rigid bodywith 6 degrees of freedom in coordinate systems conforming to ISO 1151. The resulting scalar equations of motion for the projectile, free from singularities, enable simulation of the flight of projectiles fired at the whole range of gun quadrant elevation (0