Comparison of Two Aerodynamic Models for Projectile Trajectory Simulation

Abstract

The accuracy of aerodynamically controlled guided projectile simulations is largely determined by the aerodynamic model employed in flight simulations which impacts vehicle interaction with the surrounding air. In this work, the performance of projectile path following with two distinct aerodynamic models is examined for their possible influence on trajectory following accuracy. The study incorporates the path following guidance algorithm, which enables the object to navigate along a predefined path. The simulation mathematical model is developed in the MATLAB/Simulink environment. In addition, by integrating the path-following algorithm with the two aerodynamic models, the dynamic behaviour of the aerodynamically controlled projectile can be compared. This allows for a more comprehensive analysis of the trajectory and the effects of each model on the desired flight path. Further research can explore the differences between the two models in greater detail and quantify their impact on unmanned projectile trajectory predictions, in addition to further exploring the specific characteristics and limitations of each model. This will involve analysing their assumptions, computational methods, and inputs to identify potential sources of error or uncertainty in the simulations. Moreover, these results have important implications for the design of aerodynamically controlled projectiles as well as a deeper understanding of aerodynamic mathematical modelling in flight simulation.

1. Introduction

Aerodynamically designed vehicles, such as projectiles and other unmanned airborne systems, are essential for a variety of purposes, such as exploration, transportation, and communications. The precision of their trajectory, which is a result of different aerodynamic characteristics, has a substantial impact on their performance. Aerodynamic investigation and analysis are core principles of design and development. Nearly every aspect of a projectile’s performance is influenced by aerodynamics; this ranges from range to stability and precision to manoeuvrability. Developing an effective aerial vehicle without an in-depth understanding of the effects imposed by the aerodynamic forces on moments on the projectile may prove challenging. Several popular methods were used to analyse the aerodynamic performance of flying vehicles; a popular method was the use of aerodynamic wind tunnels. The popularity of this method arose due to its ability to create a realistic flow-controlled environment with relatively precise measurements of aerodynamic forces, moments, pressure distribution and characteristics of the flow around the object. The use of this method for analysing projectile design was presented by several studies, such as [1,2,3,4]. The work conducted by G. Chen. [5] based on wind tunnel testing addressed the aerodynamic modelling problem of an axis-symmetric air-to-air projectile with small-scale canards, resulting in the development of three aerodynamic force models. However, this method presents disadvantages when it comes to the infrastructure necessary to conduct such analysis and the costs associated with it. The model of the object in the tunnel is limited by its size, and therefore scaled-down versions of the vehicle are often used for testing. This may cause a mismatch in the characterisation of flow conditions, potentially affecting the accuracy of the results.

Nowadays, the use of computational numerical methods for aerodynamic analysis is common, especially with the Computational Fluid Dynamics (CFD) approach. This approach allows the performance of detailed examinations of the flow, including complex phenomena such as shock waves, turbulent flows and flow separations. Several studies, such as in [6], employed this method to investigate the aerodynamic performance of aerial vehicles as well as perform investigations on the multi-body separation interference of kinetic energy projectile [7]. On the other hand, the CFD approach presents a key challenge concerning the required computational resources when performing highly complex scenarios. Generally speaking, the most reliable method is flight testing. This approach offers the most realistic evaluation of performance and validation, as was demonstrated in [8] for the validation of flight dynamic stability optimisation constraints with flight tests. An additional study was performed in [9], where flight tests and aerodynamic model validation were conducted for the Eurofighter Typhoon aircraft. Another notable study presented in [10] assessed the accuracy of aerodynamic wake encounter models with flight tests, and in [11], the flight tests performed for the NASA Hyper-X programme were described.

In the early stages of projectile design, researchers and engineers often use analytical [12,13] and semi-empirical methods like DATCOM (Data Compendium, a standard software tool used for estimating the aerodynamic stability and control characteristics) [14] for quick estimation of aerodynamic coefficients based on geometry and flight conditions [15]; this is useful for preliminary analysis and design comparison. The method assists in refining design elements which may significantly impact performance, leading to an efficient and effective design. In addition, it offers a simpler alternative for reasonable approximations in comparison to wind tunnel testing, Computational Fluid Dynamics and flight tests. Moreover, the method is not without drawbacks; several complex aerodynamic phenomena might not be captured as accurately as full-scale simulations and testing. This method depends on databases, which, for certain projectile configurations and flight conditions, are not well represented; it might prove to be less effective and potentially oversimplifying. The DATCOM method was utilised in the research for the design optimisation of long-range projectiles [16]. A weighted multi-objective particle swarm optimisation technique was used in research to determine the control surface sizing that maximised lift-to-drag, minimised drag, and met a static margin value for the vehicle at a specific body angle of attack. Several studies combined the outlined methods to perform a detailed analysis of aerodynamic performance and comparison between these practices, such as [17] in which wind tunnel experiments, CFD and theoretical methods were compared for projectile aerodynamics. Similarly, the work conducted in [18] presented the comparison of DATCOM, CFD and wind tunnel results for a ballistic projectile. The methods were adopted in the determination of longitudinal coefficients in supersonic flight conditions. Moreover, the research conducted in [19] compared the results of aerodynamic performance predictions using DATCOM and CFD for the SA-2 projectile while other research concentrated on predicting and validating the estimates computed by DATCOM and Aeroprecition codes with a wind tunnel [20]. Given the necessity for efficient and comparative research in aerodynamic modelling, this study adopted the DATCOM approach. The method’s ability to provide predictions of aerodynamic coefficients based on the geometry of the projectile in a fast and effective manner made it especially suitable for the purpose of this study, allowing a quick comparison of various aerodynamic models to evaluate their performance.

During the design process of such aerospace systems, it is a common and crucial practice to develop flight simulation software that allows for engineers to test and validate the performance and design of the system under various conditions and uncertainties. The practice has proven beneficial over the years in reducing costs and the probable risk arising from testing physical prototypes in real life. By employing these models, early identification of potential flaws and performance issues may be possible. Safety aspects may be taken into account as additional benefits due to the ability to perform further studies in extreme conditions, pushing the limits of the design in a controlled environment, akin to the research performed in [21] which presented computer simulation results for two types of a 122 mm calibre projectile. The MATLAB/Simulink environment was used to create various simulation mathematical models; it is well known for its resilience and practicality in dealing with complicated dynamic systems. Several studies employed this environment for modelling flight simulation during research, such as developing a flight simulation model for a modified Hypersonic X-15 aircraft [22], for the Javelin projectile system [14] and research concerning the use of a vertical cold launch system with rapid pitch manoeuvres to increase range and enhance projectile coverage [23]. Moreover, the development of the six-degrees-of-freedom model was utilised to perform a parametric guidance study for a 160 mm projectile steered by lateral thrusters [24]. Furthermore, MATLAB/Simulink might be used for the purpose of aerodynamic modelling, as presented in [25]. The work developed an algorithm capable of computing the aerodynamic derivatives of a supersonic projectile and developed a flight simulation model.

Following this, MATLAB/Simulink is chosen for the development of the flight simulation model in this study, which examines the dynamic behaviour of a guided projectile under varied flight conditions. This is achieved by combining a path-following algorithm with two independent aerodynamic models. The comparison of these models is critical for understanding their impact on projectile performance and flight simulations, as it enables a thorough assessment of trajectory and performance implications. Such comparative studies are essential for advancing the understanding of aerodynamic models in flight simulations and for identifying potential inaccuracies or uncertainties. Building upon established methodologies, this study introduces innovative aspects to the field of aerodynamically controlled unmanned projectiles. By uniquely integrating the path-following algorithm with two distinct aerodynamic models within the MATLAB/Simulink environment, in this study, we facilitate a direct and quantitative comparison of their effects on projectile trajectory and performance. The modified aerodynamic model in particular represents a significant advancement in accurately simulating canard aerodynamics and understanding the effects of aerodynamic coefficients and damping under various flight conditions, including the presence of wind. This enhanced model provides a more realistic and precise representation of projectile dynamics, especially in controlled flight scenarios. The application of this model in comparative trajectory analysis under different test scenarios is a notable contribution, enhancing the accuracy of simulations and providing valuable insights for the design and optimisation of aerodynamically controlled projectiles, thereby pushing the boundaries of current aerodynamic modelling techniques.

The organisation of the remainder of this work is described as follows: Section 2 presents a modified guided version of a projectile that is used as the test vehicle for the study. Next, the flight simulation and the investigated aerodynamic models are outlined in Section 3. This is followed by Section 4, in which the guidance and control system is presented. Section 5 outlines the investigated research scenarios, which is then followed by Section 6 for result presentation and discussion. Finally, the study ends with conclusions and future suggestions.

2. Test Bed Vehicle

The object used in the study (depicted in Figure 1) was based on the parameters of an existing vehicle presented in [26]. The primary modifications involved replacing the original curved fins with standard swept fins of a similar span to ensure static stability. In addition, canard control surfaces were added forward of the centre of mass to facilitate the desired aerodynamic control necessary for an aerodynamically guided projectile. These modifications were introduced to tailor the object to the specific requirements of the research. Furthermore, the modified object retained the same inertial parameters, characteristic lengths (see

Table 1. Projectile’s main parameters [26].

Parameter	Value	Units
$m_f$	60.51	kg
$m_0$	35.71	kg
L	2.7	m
$d_{ref}$	0.122	m
$I_{x{x}_f}$	0.144	$\mathrm{kg}·{\mathrm{m}}^2$
$I_{y{y}_f}=I_{z{z}_f}$	34.037	$\mathrm{kg}·{\mathrm{m}}^2$
$I_{x{x}_0}$	0.098	$\mathrm{kg}·{\mathrm{m}}^2$
$I_{y{y}_0}=I_{z{z}_0}$	23.355	$\mathrm{kg}·{\mathrm{m}}^2$

), and propulsion thrust curve (see Figure 2).

3. Flight Simulation Model

The equations of motion describing the object’s non-linear dynamics were resolved with respect to the centre of mass using the MATLAB/Simulink six-degrees-of-freedom (6 DOF) mathematical model in the flat earth ($NED$) frame of reference [27]. The quaternion-based variant of the model was chosen for attitude representation in space relative to the body-fixed frame coordinate system $O_b{x}_b{y}_b{z}_b$ in terms of mathematical computations. The choice of quaternion representation enhances accuracy in three-dimensional simulations and gimbal lock avoidance. However, for the computation of gravity force components and for the clarity of result presentation, Euler angles $\varphi$, $\theta$ and $\psi$ were employed due to their intuitive interpretability. The relationship between the Euler angles, the body coordinate system and the $NED$ reference coordinates is depicted in Figure 3. Furthermore, Equations (1) and (2) which were developed on the basis of change theorems of linear momentum and angular momentum define the projectile’s equations of motion.

$${\overrightarrow{\dot{V}}}_b={\displaystyle \frac{1}{m}}({\overrightarrow{F}}_A+{\overrightarrow{F}}_G+{\overrightarrow{F}}_T)-{\overrightarrow{\omega }}_b\times {\overrightarrow{V}}_b$$

(1)

$${\overrightarrow{\dot{\omega }}}_b=I^{-1}{\overrightarrow{M}}_A-{\overrightarrow{\omega }}_b\times \left(I{\overrightarrow{\omega }}_b\right)$$

(2)

Meanwhile, Equation (3) computes the velocity of the object relative to the ground coordinate system by transforming the object’s velocity from the body frame to the reference frame. This required a transformation matrix to convert the orientation of the vehicle from the body-fixed frame to the North–East–Down ($NED$) frame of reference. The directional cosine matrix $C_B^{NED}$, which captures the rotational transformations governed by the Euler angles, is defined by Equation (4).

$${\overrightarrow{\dot{r}}}_B=C_B^{NED}·{\overrightarrow{V}}_b$$

(3)

$$C_B^{NED}=\left[\begin{array}{ccc}cos\theta cos\psi & cos\theta sin\psi & -sin\theta \\ sin\varphi sin\theta cos\psi -cos\varphi sin\psi & sin\varphi sin\theta sin\psi +cos\varphi cos\psi & sin\varphi cos\theta \\ cos\varphi sin\theta cos\psi +sin\varphi sin\psi & cos\varphi sin\theta sin\psi -sin\varphi cos\psi & cos\varphi cos\theta \end{array}\right]$$

(4)

Finally, the time derivative of the quaternion which describes the mathematical orientation of the object in a three-dimensional space, is described by Equation (5). Moreover, Equation (6) is utilised to convert the orientation representation from quaternion to Euler angles.

$$\dot{\mathbf{q}}=\frac{1}{2}\left(\left(\begin{array}{c}0\\ {\overrightarrow{\omega }}_b\end{array}\right)\times \mathbf{q}\right)$$

(5)

$$\left[\begin{array}{c}\varphi \\ \theta \\ \psi \end{array}\right]=\left[\begin{array}{c}arctan\left({\displaystyle \frac{2(q_0{q}_1+q_2{q}_3)}{1-2(q_1^2+q_2^2)}}\right)\\ arcsin\left(2(q_0{q}_2-q_3{q}_1)\right)\\ arctan\left({\displaystyle \frac{2(q_0{q}_3+q_1{q}_2)}{1-2(q_2^2+q_3^2)}}\right)\end{array}\right]$$

(6)

By performing numerical integration of Equations (1)–(3) and (5), with the additional utilisation of Equations (4) and (6), one can determine the object’s states [28,29,30,31].

3.1. Environment

Environmental factors are crucial in determining the performance, effectiveness, and operational planning of airborne vehicles. The effects are most significant when it comes to the flight path trajectory and behaviour of the projectile; this, in turn, requires the simulation model to take these factors into account. Accurately modelling all environmental influences such as gravity, electromagnetic interferences, atmospheric and weather conditions presents a challenge due to their complexity in real-life scenarios. A common approach was to concentrate on the most crucial factors of the environment, which satisfy the main objectives of the study, during the development of the flight simulation model. For the purposes of this study, only gravity, atmospheric parameters and wind were included. To that end, atmospheric conditions were modelled using the International Standard atmospheric model [32]. The model computes atmospheric parameters such as speed of sound, air density, state pressure and temperature at a given flight altitude of the vehicle (see Equations (7)–(10)).

$$\rho \left(H\right)={\displaystyle \frac{p\left(H\right)}{RT\left(H\right)}}$$

(7)

$$p\left(H\right)=p_0{e}^{-{{\displaystyle \frac{a_g}{RT}}}_0(H-H_0)}$$

(8)

$$T\left(H\right)=T_0+L_0(H-H_0)$$

(9)

$$a\left(H\right)=\sqrt{kRT\left(H\right)}$$

(10)

The Taylor series WGS84 gravity model was chosen to compute the gravitational acceleration $a_g$ [33]. This model requires the position in geodetic latitude, longitude, and altitude as inputs. The gravitational load vector in the body-fixed frame was then computed by Equation (11), taking into account the orientation of the projectile and the computed gravitational acceleration.

$${\overrightarrow{F}}_G=\left[\begin{array}{c}-m{a}_{g}sin\left(\theta \right)\\ m{a}_{g}cos\left(\theta \right)sin\left(\varphi \right)\\ m{a}_{g}cos\left(\theta \right)cos\left(\varphi \right)\end{array}\right]$$

(11)

Finally, the Horizontal Wind Model 14 was used for wind effects [34]. The model provides zonal and meridional winds for a specified geodetic position. The wind speed was converted to the body-fixed frame, and the relative velocity of the vehicle to the wind was computed to determine the aerodynamic forces and moments.

3.2. Aerodynamics

The aerodynamic models employed in the flight simulation model are crucial as they govern the interaction of the vehicle with its surroundings, altering its flight path, performance and precision. It is essential for identifying the optimal shape of the object and control surface configurations. Figure 4 presents the notation used for control surfaces. It additionally ensures accurate prediction of stability and control, which are essential for safety and reliability aspects. As previously mentioned, two aerodynamic models were employed in the flight simulation model. This comparison was aimed at assessing their strengths and weaknesses under a variety of flight conditions. The total aerodynamic coefficients were obtained via the DATCOM semi-empirical method [35,36] for the geometry of the object used in the study. This method provided estimates across a range of flight conditions, in particular for Mach numbers ranging from 0 to 4 and angles of attack ranging from 0 to 15 degrees. The static and dynamic derivatives for each model were fitted using the least squares method [37], a statistical method to reduce the discrepancies, ensuring both models align well with the outputs from DATCOM.

The first investigated aerodynamic model, denoted as the “Classic” model, adopts a traditional approach to modelling aerodynamic forces, moments and the influence of control surfaces. This model incorporates fundamental aerodynamic principles [38,39], enhanced by Taylor series expansions for aerodynamic coefficients similarly to [28,40,41], which is a common approach in modelling the aerodynamic coefficients of aerospace vehicles. The equations defining the aerodynamic forces and moments for the classic model are presented by Equations (12) and (13). Moreover, the values of the coefficients were computed in the aerodynamic coordinate system; this frame rotates with respect to the body-fixed frame by incidence angles $\alpha$ and $\beta$ as presented in Figure 5.

$${\overrightarrow{F}}_A={\displaystyle \frac{1}{2}}\rho v_A^2{S}_{ref}{C}_A^B\left[\begin{array}{c}-C_D\\ C_S\\ -C_L\end{array}\right]$$

(12)

$${\overrightarrow{M}}_A={\displaystyle \frac{1}{2}}\rho v_A^2{S}_{ref}{d}_{ref}\left[\begin{array}{c}{C}_l\\ C_m\\ C_n\end{array}\right]$$

(13)

The detailed modelling of the aerodynamic coefficients for the classic model is defined by Equations (14)–(19). These equations take into account various aspects of aerodynamic forces and moments, such as drag, side force, lift, rolling moment, pitching moment and yawing moment. The model offers a strong foundation for comprehending the projectile’s aerodynamic behaviour, especially in relation to the static and dynamic derivatives as a function of flight conditions, the influence of control surfaces and aerodynamic damping on the flight dynamics of the object.

$$C_D(M,\alpha ,\delta )=C_{D_0}\left(M\right)+{\displaystyle \frac{C_L^2(M,\alpha ,\delta )}{e·\Lambda ·\pi }}+\sum _{i=1}^{N_{\delta }}{\displaystyle \frac{\partial C_D}{\partial {\delta }_i}}\left|{\delta }_i\right|$$

(14)

$$C_Y(M,\beta ,\delta )=C_{Y_{\beta }}\left(M\right)·\beta -{\displaystyle \frac{\partial C_Y}{\partial {\delta }_2}}\left(M\right){\delta }_2+{\displaystyle \frac{\partial C_Y}{\partial {\delta }_4}}\left(M\right){\delta }_4$$

(15)

$$C_L(M,\alpha ,\delta )=C_{L_{\alpha }}\left(M\right)·\alpha -{\displaystyle \frac{\partial C_L}{\partial {\delta }_1}}\left(M\right){\delta }_1+{\displaystyle \frac{\partial C_L}{\partial {\delta }_3}}\left(M\right){\delta }_3$$

(16)

$$C_l(M,p)=\sum _{i=1}^{N_{\delta }}{\displaystyle \frac{\partial C_l}{\partial {\delta }_i}}{\delta }_i$$

(17)

$$\begin{array}{c}\hfill C_m(M,\alpha ,\dot{\alpha },q,\delta )=C_{m_{\alpha }}\left(M\right)·\alpha +C_{m_q}\left(M\right){\displaystyle \frac{d_{ref}·q}{2|v_A|}}+C_{m_{\dot{\alpha }}}\left(M\right){\displaystyle \frac{d_{ref}·\dot{\alpha }}{2|v_A|}}-{\displaystyle \frac{\partial C_m}{\partial {\delta }_1}}\left(M\right){\delta }_1+{\displaystyle \frac{\partial C_m}{\partial {\delta }_3}}\left(M\right){\delta }_3\end{array}$$

(18)

$$\begin{array}{c}\hfill C_n(M,\beta ,\dot{\beta },r,\delta )=C_{n_{\beta }}\left(M\right)·\beta +C_{n_r}\left(M\right){\displaystyle \frac{d_{ref}·r}{2|v_A|}}+C_{n_{\dot{\beta }}}\left(M\right){\displaystyle \frac{d_{ref}·\dot{\beta }}{2|v_A|}}-{\displaystyle \frac{\partial C_n}{\partial {\delta }_2}}\left(M\right){\delta }_2+{\displaystyle \frac{\partial C_n}{\partial {\delta }_4}}\left(M\right){\delta }_4\end{array}$$

(19)

Following the description of the “Classic” aerodynamic model, the second model investigated in the study was an enhanced aerodynamic model more suitable for projectile simulations by taking into account the periodicity of the aerodynamic incidence angles for abnormal flight conditions. This model was denoted as the “Modified” aerodynamic model in the study. Furthermore, the model offers a more detailed and accurate representation of control surface dynamics by taking into account the interference effects of control surface vortex shedding. This is crucial for enhanced predictions of a projectile’s manoeuvrability and response to control inputs, with a full description of the model found in [6]. The equations defining the aerodynamic forces and moments for the modified model are Equations (20) and (21).

$${\overrightarrow{F}}_A={\displaystyle \frac{1}{2}}\rho v_A^2{S}_{ref}\left(C_{A^B}\left[\begin{array}{c}-C_D\\ C_S\\ -C_L\end{array}\right]+\sum _{i=1}^{N_{\delta }}{C}_{C_i}^B{C}_{A{C}_i}^{C_i}\left[\begin{array}{c}-\Delta C_D\\ \Delta C_S\\ -\Delta C_L\end{array}\right]\right)$$

(20)

$${\overrightarrow{M}}_A={\displaystyle \frac{1}{2}}\rho v_A^2{S}_{ref}{d}_{ref}\left(\left[\begin{array}{c}{C}_l\\ C_m\\ C_n\end{array}\right]+\sum _{i=1}^{N_{\delta }}{C}_{C_i}^B\left[\begin{array}{c}\Delta C_l\\ \Delta C_m\\ \Delta C_n\end{array}\right]\right)$$

(21)

The relationship between the aerodynamic coordinate system of this model and the object’s reference frame is presented in Figure 6. The modelling of the aerodynamic coefficients was performed using Equations (22)–(27) while the control surface aerodynamic coefficients were modelled using Equations (28) and (33).

$$C_D(M,{\alpha }^{\prime })=C_{D_0}\left(M\right)+C_{D{\alpha }^{\prime 2}}\left(M\right)si{n}^2{\alpha }^{\prime }$$

(22)

$$C_Y(M,{\alpha }^{\prime },p)=C_{Y_0}\left(M\right)+C_{Y_{p{\alpha }^{\prime }}}\left(M\right)sin{\alpha }^{\prime }·\frac{p{d}_{ref}}{2v_A}$$

(23)

$$C_L(M,{\alpha }^{\prime },\dot{{\alpha }^{\prime }},q^{*})=C_{L_0}\left(M\right)+\frac{1}{2}{C}_{L_{{\alpha }^{\prime }}}\left(M\right)sin2{\alpha }^{\prime }+C_{L_{q^{*}}}\left(M\right)\frac{q^{*}{d}_{ref}}{2v_A}+C_{L_{{\alpha }^{\prime }}}\left(M\right)\frac{{\alpha }^{\prime }{d}_{ref}}{2v_A}$$

(24)

$$C_l(M,p)=C_{l0}\left(M\right)+C_{l_p}\left(M\right)\frac{p{d}_{ref}}{2v_A}$$

(25)

$$\begin{array}{c}\hfill C_m(M,{\alpha }_{xz},\dot{{\alpha }_{xz}},{\phi }^{\prime },p,q)=C_{m0}\left(M\right)+C_{m_{{\alpha }_{xz}}}\left(M\right)sin{\alpha }_{xz}+C_{m_{{\alpha }_{xz}^3}}\left(M\right)si{n}^3{\alpha }_{xz}+\\ \hfill C_{m_q}\left(M\right)\frac{q{d}_{ref}}{2v_A}+C_{m_{\dot{{\alpha }_{xz}}}}\left(M\right)\frac{\dot{{\alpha }_{xz}}{d}_{ref}}{2v_A}-C_{n{p}_{\alpha }}\left(M\right)sin{\alpha }^{\prime }sin{\phi }^{\prime }\frac{p{d}_{ref}}{2v_A}\end{array}$$

(26)

$$\begin{array}{c}\hfill C_n(M,{\alpha }_{xy},\dot{{\alpha }_{xy}},{\phi }^{\prime },p,r)=C_{n0}\left(M\right)+C_{n_{{\alpha }_{xy}}}\left(M\right)sin{\alpha }_{xy}+C_{n_{{\alpha }_{xy}^3}}\left(M\right)si{n}^3{\alpha }_{xy}+\\ \hfill C_{n_r}\left(M\right)\frac{r{d}_{ref}}{2v_A}+C_{m_{\dot{{\alpha }_{xy}}}}\left(M\right)\frac{\dot{{\alpha }_{xy}}{d}_{ref}}{2v_A}-C_{n{p}_{\alpha }}\left(M\right)sin{\alpha }^{\prime }cos{\phi }^{\prime }\frac{p{d}_{ref}}{2v_A}\end{array}$$

(27)

$$\Delta C_D=C_{D_{\delta 2}}{sin}^2\left({\delta }_{eq,i}\right)$$

(28)

$$\Delta C_S=C_{S_{\delta 2}}{sin}^2\left({\delta }_{eq,i}\right)$$

(29)

$$\Delta C_L=\frac{1}{2}{C}_{L_{\delta }}sin\left(2{\delta }_{eq,i}\right)$$

(30)

$$\Delta C_l=\frac{1}{2}{C}_{l_{\delta }}sin\left(2{\delta }_{eq,i}\right)$$

(31)

$$\Delta C_l=\frac{1}{2}{C}_{m_{\delta }}sin\left(2{\delta }_{eq,i}\right)$$

(32)

$$\Delta C_n=C_{n_{\delta 2}}{sin}^2\left({\delta }_{eq,i}\right)$$

(33)

A key dynamic change during the projectile’s boost phase is the shift in its centre of mass ($\Delta c_g$) resulting from the consumption of the propellant. This shift presents a critical aspect of projectile dynamics and affects its control and stability. As the propellant burns and the mass at the aft part of the projectile decreases, the centre of mass shifts forward. This results in an increase in static stability and a more pronounced nose-down pitching moment, which leads to a reduction in manoeuvrability as the object becomes less responsive to control inputs. Moreover, it is worth noting for the study that a forward shift of 0.172 m occurred after 3.6 s of motor operation; this was interpolated during the simulation. Modelling this behaviour in the flight simulation increases the accuracy in predicting and managing the vehicle’s flight dynamics; this modelling was achieved using Equation (34), which quantifies the change in aerodynamic moments relative to propellant consumption at each simulation time step (denoted by i).

$$\left[\begin{array}{c}{C}_{l_{i+1}}\\ C_{m_{i+1}}\\ C_{n_{i+1}}\end{array}\right]=\left[\begin{array}{c}{C}_{l_i}\\ C_{m_i}\\ C_{n_i}\end{array}\right]+\left({\displaystyle \frac{{\left[\begin{array}{ccc}-\Delta c_{g_x}& 0& 0\end{array}\right]}^T}{d_{ref}}}\times C_A^B\left[\begin{array}{c}-C_{D_i}\\ C_{Y_i}\\ -C_{L_i}\end{array}\right]\right)$$

(34)

In the determination of aerodynamic characteristics for this study, the reference area for the calculation of aerodynamic forces was defined by the nominal cross-sectional area, derived from the reference diameter of the projectile. This approach affirmed that the forces were normalised in a manner related to the object’s physical dimensions. Moreover, the reference diameter also served as the reference length for determining aerodynamic moments. Utilising the diameter as the reference length, moment coefficients may be consistently and meaningfully compared, which aided in the comprehension of how the projectile rotates under different flight conditions. The aerodynamic derivatives with common values across both models as function of Mach number were presented in

Table 2. Static and dynamic derivatives common for both models as function of the Mach number.

M	${\mathit{C}}_{\mathit{D}0}$	${\mathit{C}}_{{\mathit{m}}_{\mathit{q}}}={\mathit{C}}_{{\mathit{m}}_{\dot{\mathit{\alpha }}}}={\mathit{C}}_{{\mathit{n}}_{{\dot{\mathit{\alpha }}}_{\mathbf{xz}}}}={\mathit{C}}_{{\mathit{n}}_{\mathit{r}}}={\mathit{C}}_{{\mathit{n}}_{\dot{\mathit{\beta }}}}={\mathit{C}}_{{\mathit{n}}_{{\dot{\mathit{\alpha }}}_{\mathbf{xy}}}}$
0.2	0.388	−789.640
0.4	0.334	−775.414
0.6	0.320	−760.700
0.8	0.443	−744.123
1	0.696	−767.924
1.2	0.789	−692.697
1.4	0.641	−578.439
1.6	0.558	−498.128
1.8	0.504	−445.794
2	0.465	−416.681
2.2	0.442	−383.414
2.4	0.419	−365.026
2.6	0.370	−335.349
2.8	0.351	−317.579
3	0.333	−300.399
3.2	0.321	−285.587
3.4	0.307	−279.702
3.6	0.294	−272.634
3.8	0.286	−262.003
4	0.276	−254.941

. Moreover, the Classic model’s coefficients in the form of static, dynamic, and control surface derivatives as a function of the Mach number are shown in

Table 3. Classic aerodynamic model static and dynamic derivatives as function of the Mach number.

M	${\mathit{C}}_{{\mathit{Y}}_{\mathit{\beta }}}={\mathit{C}}_{{\mathit{L}}_{\mathit{\alpha }}}$	${\mathit{C}}_{{\mathit{m}}_{\mathit{\alpha }}}={\mathit{C}}_{{\mathit{n}}_{\mathit{\beta }}}$	${\displaystyle \frac{\mathit{\partial }{\mathit{C}}_{\mathit{D}}}{\mathit{\partial }{\mathit{\delta }}_{\mathit{i}}}}$	${\displaystyle \frac{\mathit{\partial }{\mathit{C}}_{\mathit{Y}}}{\mathit{\partial }{\mathit{\delta }}_{\mathit{i}}}}={\displaystyle \frac{\mathit{\partial }{\mathit{C}}_{\mathit{L}}}{\mathit{\partial }{\mathit{\delta }}_{\mathit{i}}}}$	${\displaystyle \frac{\mathit{\partial }{\mathit{C}}_{\mathit{l}}}{\mathit{\partial }{\mathit{\delta }}_{\mathit{i}}}}$	${\displaystyle \frac{\mathit{\partial }{\mathit{C}}_{\mathit{m}}}{\mathit{\partial }{\mathit{\delta }}_{\mathit{i}}}}={\displaystyle \frac{\mathit{\partial }{\mathit{C}}_{\mathit{n}}}{\mathit{\partial }{\mathit{\delta }}_{\mathit{i}}}}$
0.2	21.181	−124.737	0.321	1.814	0.932	14.365
0.4	21.475	−125.694	0.330	1.839	0.945	14.565
0.6	22.129	−128.947	0.347	1.880	0.966	14.888
0.8	23.061	−137.943	0.371	1.936	0.995	15.333
1	25.588	−155.583	0.621	2.162	1.112	17.128
1.2	25.798	−164.733	0.808	2.189	1.126	17.343
1.4	22.036	−129.364	0.915	1.942	0.998	15.382
1.6	19.286	−100.611	0.919	1.740	0.894	13.773
1.8	17.451	−84.123	0.880	1.555	0.798	12.302
2	16.421	−72.245	0.856	1.417	0.728	11.214
2.2	15.213	−63.688	0.789	1.272	0.653	10.066
2.4	14.547	−56.996	0.772	1.187	0.609	9.392
2.6	13.470	−49.429	0.733	1.106	0.568	8.749
2.8	12.811	−44.469	0.702	1.035	0.531	8.182
3	12.168	−40.372	0.655	0.955	0.490	7.554
3.2	11.607	−37.287	0.645	0.917	0.470	7.248
3.4	11.391	−32.316	0.612	0.853	0.438	6.743
3.6	11.127	−30.966	0.630	0.841	0.432	6.652
3.8	10.720	−27.443	0.585	0.781	0.401	6.177
4	10.451	−24.934	0.588	0.762	0.391	6.025

, while the coefficients for the Modified model are similarly outlined in

Table 4. Modified aerodynamic model static and dynamic derivatives as function of the Mach number.

M	${\mathit{C}}_{{\mathit{D}}_{{\mathit{\alpha }}^{\prime 2}}}$	${\mathit{C}}_{{\mathit{L}}_{{\mathit{\alpha }}^{\prime }}}$	${\mathit{C}}_{{\mathit{m}}_{{\mathit{\alpha }}_{\mathbf{xz}}}}={\mathit{C}}_{{\mathit{n}}_{{\mathit{\alpha }}_{\mathbf{xy}}}}$	${\mathit{C}}_{{\mathit{m}}_{{\mathit{\alpha }}_{\mathbf{xz}}^3}}={\mathit{C}}_{{\mathit{n}}_{{\mathit{\alpha }}_{\mathbf{xy}}^3}}$	${\mathit{C}}_{{\mathit{D}}_{{\mathit{\delta }}^2}}$	${\mathit{C}}_{{\mathit{S}}_{{\mathit{\delta }}^2}}$	${\mathit{C}}_{{\mathit{L}}_{\mathit{\delta }}}$	${\mathit{C}}_{{\mathit{l}}_{\mathit{\delta }}}$	${\mathit{C}}_{{\mathit{m}}_{\mathit{\delta }}}$	${\mathit{C}}_{{\mathit{n}}_{{\mathit{\delta }}^2}}$
0.2	41.455	22.289	−124.697	−22.241	1.175	6.438	−1.914	0.983	−15.152	50.830
0.4	42.553	22.599	−125.654	−22.411	1.209	6.528	−1.940	0.997	−15.363	51.537
0.6	45.041	23.288	−128.906	−22.991	1.269	6.673	−1.983	1.019	−15.703	52.679
0.8	48.696	24.270	−137.899	−24.595	1.357	6.872	−2.042	1.050	−16.173	54.255
1	48.696	26.934	−155.533	−27.741	1.357	7.677	−2.280	1.172	−18.066	60.606
1.2	102.382	27.155	−164.680	−29.372	2.957	7.773	−2.309	1.187	−18.293	61.368
1.4	112.302	23.190	−129.322	−23.066	3.350	6.894	−2.049	1.053	−16.224	54.428
1.6	111.166	20.293	−100.579	−17.939	3.369	6.173	−1.835	0.943	−14.527	48.734
1.8	110.101	18.360	−84.095	−14.999	3.225	5.514	−1.640	0.842	−12.976	43.530
2	113.433	17.275	−72.222	−12.881	3.136	5.026	−1.495	0.768	−11.828	39.681
2.2	110.860	16.003	−63.667	−11.356	2.889	4.511	−1.342	0.689	−10.617	39.681
2.4	113.266	15.301	−56.978	−10.162	2.826	4.209	−1.252	0.643	−9.906	33.231
2.6	107.479	14.168	−49.413	−8.813	2.683	3.921	−1.167	0.599	−9.228	30.956
2.8	106.327	13.473	−44.455	−7.929	2.572	3.667	−1.091	0.560	−8.630	28.951
3	104.125	12.797	−40.359	−7.198	2.397	3.386	−1.008	0.517	−7.968	26.730
3.2	102.148	12.206	−37.275	−6.648	2.362	3.249	−0.967	0.496	−7.645	25.647
3.4	105.309	11.979	−32.306	−5.762	2.237	3.022	−0.900	0.462	−7.112	23.860
3.6	107.108	11.701	−30.956	−5.521	2.305	2.981	−0.888	0.455	−7.016	23.536
3.8	105.628	11.272	−27.434	−4.893	2.140	2.769	−0.824	0.423	−6.516	21.858
4	106.246	10.990	−24.926	−4.446	2.149	2.700	−0.804	0.412	−6.355	21.318

. It is important to note that the methodology used for calculating the aerodynamic coefficients computes the sum of $C_{m_q}+C_{m_{\dot{\alpha }}}$. Due to the combined nature of this calculation, it was assumed that this sum was evenly divided between the two components. Furthermore, for the purpose of simplifying the study, the projectile was considered to be perfectly axially symmetrical. This assumption, along with other reasons outlined in [6], led to several variables being omitted from the table because their values were assumed to be zero.

4. Guidance and Control

The Simulink way point follower block was used for directing the projectile path through a set of way points defining the commanded trajectory. The block implements the guidance algorithm which computes the look-ahead coordinates and the desired heading [42]. Key parameters of the algorithm, like the minimum look-ahead distance and the look-ahead transition radius, were set to 300 and 50 m, respectively. These parameters are essential for the algorithm and block’s operation. An autopilot system, incorporating PID controllers for roll, pitch, and yaw, was developed to stabilise and navigate the aerodynamically controlled projectile, as shown in Figure 7. The roll channel’s primary function was to maintain lateral stability by keeping the desired roll angle at zero. This was achieved by computing the commanded roll rate using the first PID controller; then, the error between the commanded roll rate and the actual roll rate was inputted to the second PID controller to compute the normalised commanded roll signal. The pitch channel comprised three PID controllers; the first computed the desired pitch angle based on altitude deviation from the guidance algorithm’s command. This was followed by the second PID controller, which outputted the required pitch rate, and the third generated the normalised command for pitch control. For yaw control, the structure followed the roll channel, with the first PID controller calculating the yaw rate and the second producing the yaw control signal based on the guidance algorithm’s commanded yaw angle. The output from these controllers, a normalised roll, pitch and yaw control signal, was then transformed into commanded control surface deflection angles by a mixer matrix (see Equation (35)) to achieve the desired flight path. The values of PID coefficients for both aerodynamic model controllers are presented in

Table 5. Classic aerodynamic model PID autopilot coefficients for roll, pitch and yaw channels.

PID Coefficients	$\mathit{\varphi }$ PID	p PID	Altitude PID	$\mathit{\theta }$ PID	q PID	$\mathit{\psi }$ PID	r PID
P	5	2.5	0.6	4	4.2	5.5	2.8
I	0.1	0.05	0.001	0	0	0.3	0.02
D	0.3	0.1	5	0.05	4	0.8	0.04

and

Table 6. Modified aerodynamic model PID autopilot coefficients for roll, pitch and yaw channels.

PID Coefficients	$\mathit{\varphi }$ PID	p PID	Altitude PID	$\mathit{\theta }$ PID	q PID	$\mathit{\psi }$ PID	r PID
P	5	2.5	0.5	2	0.5	5.3	3
I	0.1	0.05	0	0	0	0.3	0.01
D	0.3	0.1	0.5	0.05	0.3	0.5	0.05

.

$$\left[\begin{array}{c}{\delta }_1\\ {\delta }_2\\ {\delta }_3\\ {\delta }_4\end{array}\right]=\left[\begin{array}{ccc}1& -1& 0\\ 1& 0& -1\\ 1& 1& 0\\ 1& 0& 1\end{array}\right]·\left[\begin{array}{c}{\delta }_{{\varphi }_{cmd}}\\ {\delta }_{{\theta }_{cmd}}\\ {\delta }_{{\psi }_{cmd}}\end{array}\right]$$

(35)

Following this process, the servos for each canard were modelled as a closed-loop control system to manage the control surface position, taking into account its maximum deflection angle, deflection speed and time delay; this is depicted in Figure 8. Initially, the system computed the error between the actual and the commanded deflection angles. The error signal was moderated by a first-order transfer function, incorporating a time delay of $T=0.05$ s representing the first-order lag element of the system. Subsequently, the actuation rate was saturated to remain within the control surface’s physical capabilities. This rate of change was integrated to determine control surface deflection, which was again subjected to saturation to prevent exceeding the maximum deflection angle. For this study, the canards’ deflection speed was set to ±65 degrees per second and the canard deflections to ±25 degrees. Finally, the system’s actual controlled position was continuously fed back into the control loop, dynamically adjusting the response to the commanded inputs. This feedback mechanism ensured adherence to the physical constraints of the servomechanism model, maintaining the balance between the system’s response and its operational limits.

5. Investigated Research Scenarios

To obtain a cohesive understanding of the influence of the aerodynamic models on the results and performance of the guided projectile, several trajectory analysis scenarios were defined and investigated. Each scenario presents unique challenges and objectives to assess the versatility and effectiveness of the projectile in varying conditions while allowing for the advantages and limitations of each of the pre-described aerodynamic models to be highlighted. Four research analysis cases were selected for the study; three of them required predefined flight paths comprising smoothly interpolated way points as the main input to the study scenario, while one of the cases analysed an uncontrolled flight. The first investigated scenario was the research flight with a disabled control system and no wind presence (see Figure 9). The objective of this case was to examine the basic aerodynamic performance and range of the projectile.

The second case depicted in Figure 10 presents the “Altitude hold and descent” scenario with the presence of wind as well as an active guidance and control system. The investigation aimed at assessing the object’s ability to maintain a specific altitude for a designated duration before initiating a controlled descent towards the ground. The analysis has the potential to highlight altitude control accuracy, energy efficiency during the hold and precision when initiating the descent phase. In addition, the impact of the modelled aerodynamic control surface effects comes into play and the accuracy while capturing the dynamic effects of the canards on the trajectory is visible when processing the results.

The ability to navigate around obstacles such as terrain while following a defined trajectory provides more testing flexibility. This type of research case examines the vehicle’s ability to adjust its course by adjusting its motion in the longitudinal and directional planes, which involves the proper combinations of control surface deflections to create the desired aerodynamic forces and moments to carry out the desired study case. Figure 11 presents the desired flight path trajectory for the “Obstacle avoidance” test case.

The last scenario that was investigated required the object to perform more complicated manoeuvres (see Figure 12). This was performed in order to investigate the aerodynamic effects as a result of the mathematical models of the aerodynamic forces and moments, especially the effects of canard deflections. Meanwhile, this provides additional insight into the evaluation of the control algorithm’s ability in terms of response and agility of the projectile.

6. Results

The 6 DOF projectile mathematical model was developed in a MATLAB/Simulink 2023a environment. The solution to the equations of motion was obtained utilising the Runge–Kutta (RK4) solver with a frequency of 500 Hz. To optimise simulation time, no navigation sensors such as the Inertial Navigation System (INS) or the Inertial Measurement Unit (IMU) were modelled for the purpose of simulations. This was performed in order to increase the computational speed of the calculations while focusing on the main objective of the study. The output states were directly used as inputs to the non-linear model for the following time steps. Furthermore, the simulation initial conditions for all investigates scenarios are presented in

Table 7. Simulation model initial conditions.

Parameter	Value	Unit
Latitude	52.178	$°$
Longitude	20.956	$°$
Altitude	0	m
$V_0$	20	m/s
${\varphi }_0$	0	$°$
${\theta }_0$	50	$°$
${\psi }_0$	90	$°$
$p_0=q_0=r_0$	0	$°/\mathrm{s}$

.

6.1. Ballistic Flight Results

The results from the uncontrolled flight case for both of the investigated aerodynamic models were compared with existing results from [26]. To perform such comparison, identical initial velocity, pitch, inertial parameters and thrust curve were used in this study. However, the aerodynamic model and coefficients were not shown or explained in the cited literature.

During the uncontrolled flight scenario, the object reached an apogee of approximately 17 km with a maximum range of 41 km, as seen in Figure 13a. The classical and modified aerodynamic models provide identical results in terms of altitude versus range, while a slight deviation from the results presented in [26] can be observed. The maximum velocity was achieved in the first 5 s of flight with a magnitude of 1200 m/s, which was equivalent to 3.5 times the speed of sound (see Figure 13b). The magnitude of the velocity reduced to about 350 m/s as the vehicle neared its maximum apogee, followed by a slight increase to approximately 400 m/s due to gravitational effects, and then a steady decrease until ground contact.

Moreover, Figure 14a–d showcases a comparison between the states of the object for the classic and modified aerodynamic models in the uncontrolled flight scenario. It can be seen that the results are identical. This can be explained by the fact that in an uncontrolled flight with no wind, the behaviour of the object is largely governed by its initial speed, direction, propulsion forces and gravity. While the similar values of the aerodynamic coefficients throughout the phase of the flight presented in Figure 15a,b are likely due to the total values of the coefficients being the same as function of the Mach number and aerodynamic incidence angles, the difference in both models comes from the interaction between the control surfaces and the flow being modelled differently in the two models. In addition, dynamic damping is limited in a ballistic flight scenario. Such effects are significant and have been demonstrated in the other investigated scenarios.

6.2. Altitude Hold and Descent Results

The second scenario required the vehicle to gradually climb to an altitude of 5 km above ground level for a distance of 8 km. This was then followed by a descent phase, completing 30 km of the total desired distance to be travelled. It is worth noting that in the entire phase of the flight, a steady course and no deviation in roll angle were desired. Meanwhile, the control system was activated in the presence of the indicated wind profile (see Figure 16).

The forward velocity of the object reached a maximum value of 1100 m/s (≈3.2 Mach) in a similar time duration to that of the first investigated case (see Figure 17a). The effect of the wind is evident when one looks at the side velocity component; a minimum value of −4 m/s was reached after 15 s of flight for the classical model, with significant oscillations for a few seconds. The modified model achieved half the magnitude of oscillation in approximately the same time duration, while the vertical component for both models reached a minimum and maximum value of −60 m/s and slightly higher than 20 m/s, respectively, with similar tendencies, except for variation in oscillations. The orientation of the object was presented in Figure 17b. The roll angle oscillated in the range of ±0.05 degrees while the projectile attempted to maintain roll stability with similar oscillations in the two models. The pitch of the projectile was rather identical, with a slight dip in the classical model at about 17 s of flight duration. The desired heading of 90 degrees for the vehicle to fly due east was relatively maintained, with the object’s control system making an effort to encounter the presence of the wind. The classical model showed similar behaviour to that of the side velocity oscillations.

Maximum control surface deflections of ±20 degrees were achieved in the longitudinal plane (see Figure 17c). A slight oscillation was observed between the first 10 and 20 s of flight, with an increase in oscillation frequency seen from 40 s throughout the remaining flight duration. The canards in the xy plane (${\delta }_2$ & ${\delta }_4$) oscillated during the entire flight with a higher frequency and lower amplitude to maintain roll stability and counteract wind effects. A clear distinction in the precision of the flight trajectory is observed in Figure 17d–f. The classic model resulted in a maximum error of 2.5 m in the north coordinates compared with an error of 0.7 of a metre for the modified model. In addition, an overshoot in the commanded altitude can be seen in Figure 17e for both models. However, the modified model resulted in significantly less overshoot than the classical model.

6.3. Obstacle Avoidance Results

The obstacle avoidance case required the guided vehicle to perform a series of pitch and yaw control manoeuvres to meet research analysis requirements.

The velocity profile presented in Figure 18a resulted in similar values in the forward and vertical components compared to the altitude hold and descent case, while the side velocity experienced a greater change with peaks of 10 m/s and −10 m/s at approximately 7 and 12 s, respectively; it was followed by a smaller peak at 30 s for both models with a value of about −5 m/s, which indicates the direction of the side velocity. The orientation was additionally affected by a significant increase in the initial roll angle, nearly reaching a value of 1 degree counterclockwise (see Figure 18b). The stabilisation algorithm succeeded in stabilising the roll after 10 s of flight time. The changes in pitch and yaw are expected for this kind of commanded trajectory. A minor delay is observed in Figure 18c between the canard deflections in the classic model in comparison to the modified model with slightly higher amplitudes. Moreover, this time around, the deflections in the pitch channel depicted by ${\delta }_1$ and ${\delta }_3$ were bound to maximum values of −10 degrees and 20 degrees in the first 18 s of flight duration. While the canards providing directional control oscillated between ±5 degrees due to the nature of the manoeuvres performed, no noticeable difference in the motion of the test vehicle was observed in Figure 18d. On the other hand, Figure 18e presents the side view of the trajectory, where a slight error between the commanded and actual flight path was seen for both models; this error, however, is greater for the classical aerodynamic model. The three-dimensional trajectory for the obstacle avoidance analysis can be seen in Figure 18f.

6.4. Evasive Manoeuvre Results

The final investigated research case study provided a greater insight into the behaviour of the research projectile when performing a series of manoeuvres; such manoeuvres are unconventional and uncommon for an aerodynamically controlled object of such configuration. However, the main objective of this research case was to highlight critical aspects of control surface deflection models and damping derivatives in addition to their influence on the motion of such object.

The U component of the velocity remained identical to all the investigated control scenarios, as seen in Figure 19a. The changes in the side and vertical components of the velocity were remarkable, with an increase in the magnitude of V reaching an upper value of 40 m/s and a lower value of −20 m/s up to 20 s of flight duration. This was followed by large frequency oscillations between both models in the time frame of 40–55 s, with a noteworthy delay in the classical model of 5 s. Moreover, the modified model achieved much smoother overall characteristics. For the W component, the test vehicle reached maximum and minimum vertical velocity values of ±50 m/s. In this regard, the classical model exhibited greater peaks at approximately 2 s, while the modified model showed a significant peak at 10 s. Additionally, a noticeable delay in the change in vertical velocity was observed in the classical model, which also tended to have larger oscillations overall. Roll stabilisation was largely effective throughout the investigated research scenario, with the exception of oscillations with high-frequency oscillations and an amplitude of 5 degrees clockwise and counterclockwise rotations for both models. Although a delay of approximately 5 s was noticed in the response of the classic model in Figure 19b, additional delays in the same model were observed in the pitch and yaw responses with slightly greater oscillations compared to the modified model. Maximum control surface deflections in the range of ±20 degrees were achieved by all canards due to the demanding corrections required by the guidance and control algorithms. However, Figure 19c depicts time delay in the motion of the canards. When it came to the classic model, special attention was paid to the rather sharp deflections of the canards in this model. Generally speaking, stronger deflections were observed for both models between 20 and 70 s of flight when performing the rather tight turns of the scenario. A slightly larger deviation from the commanded trajectory in the xy plane of the coordinate system is seen in Figure 19d. However, this is more noticeable in Figure 19e,f. The deviation from the mission profile was evidently minor for the modified aerodynamic model.

7. Conclusions

The mathematical flight simulation model for the guided aerodynamically controlled unmanned projectile with four control surfaces in the form of canards was shown. The six-degrees-of-freedom simulation with predefined trajectory analysis cases was developed successfully in MATLAB/Simulink, which allowed for comparative analysis. In addition, the study implemented the waypoint following guidance algorithm and developed a control algorithm capable of guiding the projectile with success. Two aerodynamic models were demonstrated and implemented in the flight simulation model for each of the outlined research cases. The outcomes of the uncontrolled flight test scenario showed no noticeable difference between the classic and modified aerodynamic models; the fact that both models provide similar results in terms of altitude versus range is significant. It suggests that for the uncontrolled case without the presence of wind, the models are consistent with each other and potentially with established physical principles, which is a valuable finding. The comparison between the two aerodynamic models would primarily revolve around their prediction of drag and stability in an uncontrolled flight. The performance of the two models started to deviate when the control system was activated and in the presence of wind, starting from the altitude hold and descent test case. A slight delay in response was observed by the classic model compared to the modified model, with larger oscillations by the classic model in the side velocity and Euler angles. The modified model’s performance advantages were additionally noticed in the smaller errors produced when following the intended flight path, with the classic model producing higher values of deviations from the commanded trajectory. The response oscillations of both models were largely similar in the obstacle avoidance scenario, with the exception of a slight delay in the classic model. The manoeuvring trajectory analysis case was challenging for both models as intended, and the findings were interesting for this situation when the test object was required to perform more complicated manoeuvres. Significantly higher oscillations by the classic model were observed in the velocity components in addition to a delay in the response. Similarly, this was seen in the results of the Euler angles and control surface deflections. Finally, a higher error was observed in the trajectory plots by the classic model; these errors were quite significant in comparison to the performance of the modified model.

In summary, the modified aerodynamic model resulted in higher precision due to enhanced modelling of canard aerodynamic behaviour and to the overall enhanced modelling of aerodynamic coefficients and damping effects due to the object’s motion and wind presence. The oscillations and delays caused by the classic model resulted in larger trajectory tracking errors; these errors may provide the designer of the projectile with false conclusions regarding its performance. Therefore, the study’s findings have important implications for the design of aerodynamically controlled projectiles. This study hopes to contribute to the development of more accurate and efficient unmanned projectile designs by improving the knowledge of these aerodynamic models. Further research with the implementation of sensor models, turbulent wind models and non-flat earth frames should be considered in order to expand the understanding and implications of aerodynamic modelling in aerodynamically controlled projectile design and simulations.