Development of the Flight Dynamic Model (FDM) Using Computational Fluid Dynamic (CFD) Simulations for an Unknown Aircraft

The usage of computational fluid dynamics (CFD) has enhanced 10-fold since the last decade, especially in the area of aerospace science. In this chapter, we will focus on determining the feasibility and validity of CFD results that are plugged in flight dynamic model (FDM) to that of actual flight of an aircraft. Flight data of an actual aircraft is used to determine the aerodynamic performance of the designed FDM. In addition to this, FDM consist of various systems integration of an aircraft; however, this study will focus on aerodynamic parameter optimization. Relative analysis is carried out to validate the FDM. This will enable readers to know how CFD can be a great tool for designing FDM of an unknown aircraft.


Introduction
Stable flight dynamic modeling and designing of an aircraft is a crucial phase faced by the aviation industry. From these perspective, when it comes to on-ground training of pilots, simulative training is required before they can face actual dynamics of real flight. Smooth and stable flight is a necessary fact as the concerned pilot is not always alone, and also high risks are involved if the pilot is untrained with regard to aircraft dynamics. Therefore, ground training on flight simulators is given to pilots for a particular aircraft. However, for the particular aircraft to run on flight simulators, initially its flight dynamic model (FDM) is designed.

Why one should design FDM?
It's a question of interest, why one should design FDM?, because the FDM is a heart for flight simulator [1]. The current dire need of flight simulators is captured from the high rate of hiring of new pilots and indulgence of the airlines whether it is in Pakistan or international forum. As the hiring and training of new pilots is expensive, for this reason internationally, ground training with respect to particular aircraft is being catered by flight simulators. For this reason, high-fidelity FDM is a concerning aspect. In line with the above, it is significantly necessary to train new pilots from a realistic approach, keeping in mind the existing piloting reviews from old fellows of that aircraft to adopt for a dynamically changing environment whether that be in terms of standard operating procedures (SOPs) of flight, systems, or navigation [2]. Moreover, the need of designing a realistic FDM for flight simulators can also reduce the amount of actual flight time pilots put on aircrafts by which fuel and CO 2 emission can be saved. This effort will inline flight organizations to act according to the ICAO Programme of Action on International Aviation and Climate Change, which enforces the ways to save material cost, i.e., fuel economy, and protect the environment which are the key concerns [1] that can be a driving factor for designing flight simulators.

What actually is FDM?
The era of FDM has changed the pathway of flight simulations. FDM consist description of flight model of a certain aircraft, which are the propulsion, navigation, controls, avionics, and aerodynamic data. The major part of this FDM is the aerodynamic data, as it handles the attitude behavior of the aircraft during flight. This aerodynamic data is acquired from either wind tunnel testing or computational fluid dynamics (CFD). Acquiring from wind tunnel is a difficult task as it is timeconsuming and requires scaling down of the design, as it is hard to test actual size model since they are bigger than the wind tunnel test sections. Therefore, CFD is preferred [2] (Figure 1).

CFD as a visual and graphical quantification tool
CFD is a tool for testing and quantifying fluid dynamics over an internal or external body where fluid flow is involved. CFD saves material cost and manufacturing time for analyzing aerodynamics of an initial design concept that has been created [3]. Aerodynamics is a study of airflow in either the internal or external side of the body. The aerodynamics concerned in this chapter is related to the external flow over an unknown aircraft body. CFD helps to simulate actual size of computer-aided design (CAD) model, in an enclosed control volume. Any property in the control volume is controlled using Reynolds transport theorem as shown in Eq. (1), and further this approach is applied on velocities over three-dimensional space and time using Eulerian technique as shown in Eq. (2): Partial differential equations are used for describing system of fluids (i.e., gasses and liquids) that are represented by the general laws of conservation of mass, momentum, and energy [4].
The principle of mass balance is used in light of law of conservation of mass for fluid element, and it is written in Eq. (3) [5].
where ∂ρ dt with derivation of density with time change and ∇ Â ρ Â v ð Þis time rate of change of volume of moving fluid.
The momentum equations in the x-, y-, and z-axes, respectively, are expressed in Eq. (4).
The energy equations used were derived from Navier-Stokes which depends upon the first law of thermodynamics [6]. The derivation for conservation of energy on a finite fluid element consists of a single equation which is expressed in Eq. (5).
CFD is a cost-effective and easy to use method which empowers engineers to virtually simulate and visualize the experiments carried out using wind tunnels. As far as the visualization of flow is concerned, CFD helps in depicting pattern of the fluid flow, which is difficult with regular wind tunnel experiments. However, wind tunnel experiments are expensive to conduct, and their real flow characteristics are hard to analyze due to the limitation of size of the test section for which scaling down of the geometry is required. Moreover, to determine the forces and moments in a wind tunnel, several pressure orifices are required and mounted over the model of interest to determine the pressure distribution on the surface of the model [7], due to which it is hard to set up the experiment as compared to CFD. Now keeping in mind our application, i.e., related to aerospace industry, Menter's shear stress transport (SST) model initially developed by F.R. Menter in 1994 is suitable. Moreover, according to the F.R. Menter study [8], it is noted that the SST Kω model outperforms and predicts the reduction of kinematic eddy viscosity quantity due to the adverse pressure gradient profiles in very good agreement for all x-station of a flat plate with that of the experiments. Moreover, in his study, SST Kω model is capable of predicting the accurate velocity profile charts as acquired from experimental study. In addition to this, SST Kω solves two equations, viz., turbulent kinetic energy "k" and the eddy dissipation rate "ω" which are given in Eqs. (6) and (7), where the variables are in italic and constants are in non-italic format: Variables in Eqs. (6) and (7) are as follows: i.e., ρ is the density; t is the time; u j and x j are the velocity and position vectors, respectively; μ t is the eddy viscosity; ν t is the kinematic eddy viscosity; γ is the intermittency factor; F 1 is a blending function; and β * , σ t , σ ω are constants. It is basically a combination of the KÀε model in the freestream and KÀω model near the walls of the geometry and is well suited for external aerodynamic flows around complex geometries and highly separated flows like airfoils at high angles of attack. In this study Siemens STAR-CCM+ software is used to carry out CFD analysis. This software includes numerous fluid dynamic models that are widely used in industry-level simulation requirements. Moreover, in STAR-CCM+ different wall treatment methods like all y+, low y+, and high y+ for treating boundary layer formation can be incorporated with the SST Kω model for true shear stress depiction.
For conducting CFD analysis, numerous test settings are permutated for gathering the aerodynamic data. However, out of different test scenarios, significance is given to cruise profile for which assumptions that are adopted for constraining our simulation are as follows: Altitude: 1000 m Air density: 0.9075 kg/m 3 Air viscosity: 1.581 m 2 /s Velocity: 115 m/s $ 220 knots Attitude cases: 372 Table 1 shows that the total number of cases conducted for CFD simulations were 372. In Table 1, "All" in first column first row means all control surfaces at zero deflection level. "Elevator" in first column second row shows the number of elevator deflection cases conducted. "Rudder" in first column third row shows the number of rudder deflection cases conducted, and finally "Aileron" in first column fourth row shows the number of aileron cases conducted. These cases were carried out with properties already mentioned before in Table 1.

CFD setup
Generic settings required for setting up any CFD requires the three basic processes: • Preprocessing

Preprocessing
The object of the study was C-130, and its CAD model was acquired from FlightGear database, which is an open-source platform. Preprocessing involves CAD import, generating and optimizing mesh using various techniques, physics, and environmental settings. The CAD model is imported using inbuilt feature which only supports listed file formats. In this study ".stl" format file is used from CAD software.
Note: It is worth noting here that STAR-CCM+ requires considerable computer hardware resources to work in a faster pace. Loading times, mesh generation, and simulation times are significantly reduced with improvement in hardware. It has been tested by running same simulations at different desktop configuration machines and noticed significant reduction with respect to elapsed time.
After successful importing, the CAD model is visible in current scene. The next step is to generate mesh. "STAR-CCM+ has all-around mesh generating feature that creates unstructured form fitted finite volume meshes of fluid and solid domains. Software is designed such that mesh generation is automatically informed by the surface tessellation and CAD elements defining the geometry, such as local curvature, surface proximity, and retained feature elements, and is further controlled by user-specified meshing parameters. The latter are organized into a hierarchy of global specifications and local refinements that enables precise control to achieve cell quality metrics, such as skewness, connectivity, conformity, near-wall cell properties, and growth rate with the smallest practical mesh even for exceptionally intricate geometries. STAR-CCM+ also employs a face-based solver technology uniquely designed to recognize arbitrary polyhedral cell topology [9]." For meshing user-specified parameters that were assigned are shown in Table 2, these were values adopted for a four-engine turboprop transporter aircraft. Tables 2 and 3 values were evaluated after repetitive attempts to achieve a welldefined and fine mesh at zero angle of attack steady flight. After running meshing method, software generated approximately 12 million cells in the mesh. The process of meshing was repeated with mild tweaks in values for every case of angle of attack and control surface deflections. For each case of flight profile, the number of cells in a mesh increases with changes in angle of attack, sideslip, and various configurations of control surfaces.

Processing
At this stage, the preprocessed settings are evaluated or computed using the solver which is tailored to our specific requirements with the options selected earlier. The processing depends on the stopping criteria for resolution of different number of iterations. Once the correct models and settings are chosen for physics conditions, then the simulations are processed for examination of the applied conditions. Here, the air properties, other physical conditions, result extraction for aerodynamic coefficients, and graphical depiction of iterated data were selected. Physical condition and fluid dynamics models implemented are listed below: • All y+ wall treatment The aerodynamic coefficients generated in the final report were as follows: These parameters were set up, and each case of simulation had certain control surface deflection and angle of attack. For a four-engine turboprop transporter aircraft, 372 cases were chosen. These cases proved to be sufficient for attaining high-fidelity flight dynamic model for aircraft simulator. The air speed was kept constant at 220 knots and at an altitude of 5000 feet. Air density and viscosity were set up accordingly. After all these steps, the CFD software is ready for analysis. By simply clicking on run button on the top pane, STAR-CCM+ starts running iteration and plots the results simultaneously.
The CFD software as mentioned earlier requires sufficient computer hardware to function properly. The elapsed time for analysis was almost 6 hours for a highend desktop configuration machine in year 2017. Compared to this, same analysis was done within 3 h on a high-end machine in year 2018 with new specifications. For better and faster results, cluster computers and supercomputers are used to run CFD simulations for acquiring tremendous amount of data sets.

Post-processing
The results attained for aircraft aerodynamics from the CFD simulations are then used by the FDM's aerodynamic module. Results obtained are of forces and moment coefficients for six different axes, i.e., drag, lift, side force, roll, pitch, and yaw axes, with deflected surfaces at different angles.

Mesh independency study
To acquire accuracy in attaining the CFD results, and keeping in mind the computational power, it is necessary to analyze the geometry for mesh independency study. For this the model's physical properties were defined with similar inlet velocities, and physical boundary conditions were set similar for the different cases. The model chosen was SST Kω model. Grid convergence analysis was conducted on coarse, medium, and fine mesh specifications at which C D , C L , and C M were analyzed. This is conducted to determine the effect of mesh quality on CFD results. The number of cells and simulation time for three different cases was simulated. First set was conducted with 0°angle of attack with all control surfaces non-deflected. The second set was with 0°angle of attack with elevator deflected by +5°. The results collected shown in Table 4 depict that the number of cells has a huge impact on the mesh independence and time period required for conducting simulations. The six meshes shown in Table 4 demonstrate that mesh independency was acquired as per deviation from coarse level mesh to fine settings has achieved a level of stagnation for estimated parameters of C L , C D , and C M . For the fine mesh, base size reduction for different geometrical parts was set around 6% of actual geometry.

Mesh resolution Coarse mesh Medium mesh Fine mesh
At +5°elevator on 0°AoA Conducted CFD simulations for above mesh settings are demonstrated in Figure 2.
It can be noted that the CFD results demonstrated significant velocity profiles and depiction of wake generation was considered but dominance was given to shearing stress, i.e., related to the near-wall stresses. Special focus was given to the surface shearing stress, as to capture the precise effect caused due to control surface deflections. Moreover, wake dominance can be optimized further by deploying more number of cells at the aft side of aircraft geometry with more computational power.

CFD to FDM integration
The FDM file is then processed using the FlightGear flight simulating software; this file is in .xml format. Moreover similar procedure is followed for checking under the JSBSim stand-alone module designed by Jon Berndt in 2004 [10]; however, it is embedded to the external image generation tool for visual effects. Input devices for aircraft control loading system used with JSBSim simulations were similar to the actual flight control loading system; however, with FlightGear simulation, Logitech extreme 3d edition flight joystick was used, this is further aggregated in Table 5.

Properties
Actual flight FlightGear JSBSim After the hardware is set up and FDM .xml scripted file is tested, we are able to collect different flight data of interest for quantification of flight maneuvers. In the next section, results of responses generated for angular rates that are dependent on the aerodynamic coefficients plugged in FDM file are discussed.

Results and discussion
The pitch rate variation attained from flight data recorder (FDR) data, FlightGear, and JSBSim data shows similarities from the maxima and minima values. The pitch rate data of actual aircraft, i.e., PITR, from Figure 3 depicts that during the cruise phase, aircraft pitching rate was within 1.5-1°/s. However, from the FDM perspective, i.e., for JSBSim and FG, it shows pitching motion rate maxima and minima between À2.25 to 3.5°/s and À2.25 to 1.5°/s, respectively. The depiction gives a clear understanding that the aircraft modeled using FDM is replicating the motion dynamics of the actual aircraft but with some deviation. Simulation results demonstrate greater rates than the actual aircraft because the CFD results plugged in the aircraft dynamic file are overpredicting the pitching rate of actual aircraft. Nevertheless, it does give insight to the reader that CFD can be helpful in the initial designing of FDM aerodynamic table [11]. The optimizations can be used to predict the actual behavior of the aircraft pitching rate. A note to remember for damping the motion and response dynamics of the aircraft from pitch axis, C mq and C mα̇a re the two variables that can be used for fine adjustments according to pilot's requirement. Before, it is necessary to have a correct initializing FDM model. The file that was being used with JSBSim-based simulator was introduced with the specific tables of C mq and C mα̇f or assisting for specific cases of landing and takeoff phase. This was implemented as per pilot's observations.
The roll rate variation attained from FDR, FlightGear, and JSBSim data shows high degree of resemblance from the maxima and minima values. The roll rate data of actual aircraft, i.e., ROLLR, from Figure 4 depicts that during the cruise phase, aircraft rolling rate was within À4 to 4°/s. However, from the flight dynamic model perspective, i.e., for JSBSim and FG, it shows rolling motion rate minima and maxima between À1.25 to 6.25°/s and À2.25 to 2.15°/s, respectively. The depiction gives a clear understanding that the aircraft modeled using FDM is replicating the motion dynamics of the actual aircraft. In addition, the simulative results of FG are in high degree of agreement to actual aircraft roll rate performance. However, JSBSim demonstrates greater roll rates, specifically in positive direction, than the actual aircraft, and this can be because the control dynamics on the JSBSim were being operated using a feedback-based control loading system; however, in an actual aircraft, relative wind component affects the feedback felt by pilot on the stick, hence adding an extra variability to aircraft control. Nevertheless, the roll rate performance was better than the pitching rate. For controlling the oscillation and damping effects in roll axis, C lp (i.e., coefficient of rolling moment due to damping) and C lr (i.e., effect of yaw rate on coefficient of rolling moment) are the major variables. Mostly, C lp is used for altering the damping effect. The effectiveness in rolling moment of the control surface is catered by C lda where da denotes the change in aileron deflection angle. In addition to C lda , the FDM file was constrained in such a way that the code for right and left ailerons was separately designed for giving correct effect of aerodynamic deflections. Nevertheless, in pitch similar steps were also followed for defining aerodynamics of positive elevator separately to that of the negative elevator deflection using C mde variable; however this was causing oscillatory modes at maximum deflection of positive elevators, for this two major variables were given full consideration during the pilot phase optimization which were C mq and C mα̇. These two variables assisted in changing the damping and oscillations caused due to pitching motion.
The yaw rate variation attained from FDR, FlightGear, and JSBSim data shows high degree of resemblance from the maxima and minima values. The yaw rate data of actual aircraft, i.e., YAWR, from Figure 5 depicts that during the cruise phase, aircraft yaw rate was within À2 to 2°/s. However, from the flight dynamic model perspective, i.e., for JSBSim and FG, it shows yaw rate minima and maxima between À2.5 to 1.05°/s and À6.45 to 2.45°/s, respectively. The depiction gives a clear understanding that the aircraft modeled using FDM is replicating the motion dynamics of the actual aircraft. In addition, the simulative results of JSBSim are in high degree of agreement to the actual aircraft yaw rate performance. However, FG demonstrates greater yaw rates, specifically in negative direction, than the actual aircraft, and this can be because the control loading used on the FG was being operated using a joystick control. Logitech extreme 3d edition flight joystick has a small moment in yawing direction, i.e., the maxima and minima is reached in just a slight deflection, causing it hard to deflect gradually, hence adding an extra difficulty during yawing moment while controlling through FlightGear. Therefore, it was found that the yaw rate demonstrated better performance using JSBSim than the FG. The stability derivatives that are used for optimizing the yaw performance of the aircraft were C nr , C ndr , and C nβ , where C nr is coefficient of yaw moment due to yaw rate, C ndr is the yaw moment caused due to rudder deflection, and C nβ is yaw moment caused due to sideslip angle. C nβ and C nr were adjusted using specific values; however, for C ndr table was defined as per to assist for landing and ground run effectiveness of the yaw moment due to rudder deflection.

Data quantification of the results
For the steady-state cases, the linear model with required changes according to pilot's input, satisfactory estimation of the aerodynamic response is achievable. However, if the pilot induces large-amplitude maneuvers or rapid divergences from the steady-state conditions, then nonlinear parameters need to be considered with the basic aerodynamic model. Klein and Morelli in their research state two ways of doing this: (a) using Taylor series expansion for defining nonlinear stability derivatives and (b) combining static terms and treating stability and control derivatives of aircraft as a function of explanatory variables, i.e., angle of attack, angle of sideslip, and Mach number [12]. Moreover, system identification technique is also good for validating the required results of the aerodynamic coefficients from FDR data by reverse engineering the states using observation matrix with specific input and outputs [13][14][15][16]. However, in this study as we are considering the steady-state flight dynamics, some mathematical techniques like standard deviation can be used for evaluating performance of the CFD attained variables. This quantification was carried out using Excel.
Quantification of Figures 5-13 was conducted using Excel through which the local maxima and minima were depicted on different charts for observing the exact value from specific plots. By the help of the local maxima and minima, standard deviation (SD) was also calculated using Eq. (8): where σ is the population standard deviation, N is the size of the population, and x i and μ are the population mean.
It was found that the SD lied on a discrepancy of about 1°/s from mean values, giving us a generic idea that the results deviation was nominal as seen from Table 6. It is found that the max deviation is mostly found in the JSBSim-based results. However, FG results deviation is similar to the deviation of the actual aircraft. From JSBSim perspective, it can be depicted that the atmospheric model was not completely integrated to give the effects as it was found in FlightGear. Figure 6 shows the local maxima and minima of the actual aircraft roll rate during a steady-state flight condition. The figure gives an idea of the local maxima and minima in red and green color, respectively. Using this SD was calculated to understand the difference between different operating platforms. Figure 7 shows the local maxima and minima of the JSBSim aircraft roll rate during a steady-state flight condition. The figure gives an idea of the local maxima and minima in red and green color, respectively. Using this SD was calculated to understand the difference between different operating platforms. Figure 8 shows the local maxima and minima of the FlightGear aircraft roll rate during a steady-state flight condition. The figure gives an idea of the local maxima and minima in red and green color, respectively. Using this SD was calculated to understand the difference between different operating platforms.  Table 6. Standard deviation of the angular rates.

Figure 6.
Local maxima and minima on the actual (P) roll rate in°/s chart, where green are the minima and red are the maxima. Figure 9 shows the local maxima and minima of the actual aircraft pitch rate during a steady-state flight condition. The figure gives an idea of the local maxima and minima in red and green color, respectively. Using this SD was calculated to understand the difference between different operating platforms. Figure 10 shows the local maxima and minima of the JSBSim aircraft pitch rate during a steady-state flight condition. The figure gives an idea of the local maxima and minima in red and green color, respectively. Using this SD was calculated to understand the difference between different operating platforms. Figure 11 shows the local maxima and minima of the FlightGear aircraft pitch rate during a steady-state flight condition. The figure gives an idea of the local maxima and minima in red and green color, respectively. Using this SD was calculated to understand the difference between different operating platforms. Figure 12 shows the local maxima and minima of the actual aircraft yaw rate during a steady-state flight condition. The figure gives an idea of the local maxima  and minima in red and green color, respectively. Using this SD was calculated to understand the difference between different operating platforms. Figure 13 shows the local maxima and minima of the JSBSim aircraft yaw rate during a steady-state flight condition. The figure gives an idea of the local maxima and minima in red and green color, respectively. Using this SD was calculated to understand the difference between different operating platforms. Figure 14 shows the local maxima and minima of the FlightGear aircraft yaw rate during a steady-state flight condition. The figure gives an idea of the local maxima and minima in red and green color, respectively. Using this SD was calculated to understand the difference between different operating platforms.
After observing the SD, standard error measure was also calculated. Table 7 demonstrates values of the standard error measure of angular rates of actual, JSBSim, and FlightGear responses. Standard error measure shows how varied data is acquired from the actual responses. It can be seen that the response variation from  mean is in similar constraints to each other. The standard error deviation for pitch rate (i.e., Q) is higher for JSB and FG to that of the actual model because the values chosen for the aerodynamic coefficients are generating greater pitching moment; nevertheless, here the control column perspective should not be neglected as this has also a greater impact on variability of the results.

Conclusions
The study has compiled a way for designing FDM tables using CFD obtained results and pilots response from FDR data. The performance of CFD-designed FDM was tested and compared with the FDR data with similar steady-state conditions as prevailed during actual flight scenario. It was found that the response  characteristics of the simulated angular rates were in correspondence with the actual rates. In addition to this, a standard deviation error measure was below 0.1 for all the rates from the mean position, giving us a confidence on the values of the aerodynamic coefficients plugged in the FDM. Moreover, 1°/s SD of the angular rates explained that the CFD data is useful for initially designing the FDM;  Table 7.
Standard error measure of the angular rates. however, further modifications using system identification with least square method (LSM), Taylor series expansion, and filtering methods can be employed for increasing the accuracies of the results.
ω Eddy dissipation rate u forward velocity in x-direction v lateral velocity in y-direction w normal velocity in z-direction p roll rate q pitch rate r yaw rate C L lift coeffeicient C D drag coefficient C Y side force coeffeceint C l roll moment coefficient C lp roll moment coefficient due to roll rate C lr roll moment coefficient due to yaw rate C lda roll moment coefficient due to aileron input C ldr roll moment coefficient due to rudder input C lβ roll moment coefficient due to rudder input C mq pitch moment coefficient C mά pitch moment coefficient C mde pitch moment coefficient C n yaw moment coefficient C np roll moment coefficient due to roll rate C nr roll moment coefficient due to yaw rate C nda roll moment coefficient due to aileron input C ndr roll moment coefficient due to rudder input C nβ roll moment coefficient due to rudder input C yβ roll moment coefficient due to rudder input SD standard deviation σ the population standard deviation N the size of the population x i and μ the population mean SDE standard deviation error