Analysis of uncertainty factors in the prediction of missile pitch damping derivatives by unstructured CFD solver

Based on the transient planar pitching method, the dynamic derivative calculation function is implemented on an unstructured computational fluid (CFD) software NNW-FlowStar. Then the pitch-damping dynamic derivatives of the standard model of the ANF missile at zero angles of attack and different Mach numbers (0.5∼4.5) are calculated. The influence of various numerical and physical parameters on the calculation results of dynamic derivatives is studied in detail. It is found that the dynamic derivatives of pitch damping are independent of the values of amplitude (0.125° ≤ A ≤ 1°) and reduction frequency (0.025 ≤ k ≤ 0.2). When the number of calculation steps of the sine oscillation period is 200 and the number of inner iterations is 40, the prediction accuracy and computing efficiency can meet the engineering requirements simultaneously in this case. The influence of different grid resolutions on the dynamic derivative calculation results is studied. Results show that the grid requirements of dynamic derivative calculation are consistent with those of steady-state calculations. Finally, the calculation results are compared with the free flight test results, which verifies the reliability of the methods and parameters suggested in this paper.


Introduction
The dynamic derivative is the derivative of the aerodynamic force or moment coefficient concerning the rate of change of some parameters with time (such as pitch angle velocity, yaw angle velocity, and roll angle velocity), reflecting the response of aerodynamic force and moment to the change of flight state like speed or angle of attack, including damping derivatives, cross derivatives, inertia derivatives, Magnus coefficients, etc.It is an essential original aerodynamic parameter for the design of aircraft navigation and control systems and flight dynamic quality analysis and is also crucial data for the research and analysis of aircraft stability, controllability, and flutter phenomena.So it is very important to predict it accurately for ballistic and missile weapons.
CFD is a very economical and efficient tool for predicting aerodynamic forces and moments of aircraft and can be used as a supplement to free-flight trajectory tests, wind tunnel tests, and semiempirical analyses and estimations.Linear mechanics theory shows that the pitching force and moment acting on the aircraft is related to the pitch damping force and moment [1,2].The pitching damping moment (PDM) is the damping moment that occurs in the process of aircraft pitching and swinging, which forces the aircraft to gradually weaken and disappear, thus ensuring the aircraft pitching stability.In this study, we use the transient planar pitching method [3][4][5][6][7][8][9][10][11][12][13][14] (also known as the forced oscillation approach) to calculate the PDM derivatives of aircraft.Transient planar pitching refers to the simple harmonic oscillation motion of the aircraft around its center of gravity during straight flight.In numerical methods, forced sinusoidal motion is usually used to simulate the process.This motion is time-dependent, therefore time-accurate CFD methods were used to compute the flow field solution.The dynamic aerodynamic force and moment in the flight process were obtained from the flow field solution at each time step to identify the dynamic derivative of the aircraft.
The objective of this paper is to study various factors affecting the calculation of dynamic derivatives by CFD and find the best scheme for accurately predicting dynamic derivatives.We take the Army-Navy Basic Finner (ANF) missile standard model as the research object, use the transient planar pitching method to calculate the dynamic derivatives of ANF, and study various influencing factors.An industrial CFD software NNW-FlowStar [15] is used to simulate the unsteady flow field of the aircraft in sinusoidal motion.As a standard missile model, ANF has been extensively tested in ballistic free flight and wind tunnels in foreign countries for many years [16][17][18][19][20] .Due to the limited information, the test data used to confirm the CFD calculation results in this paper are directly collected from literature [3,5].In this paper, we analyzed the characteristics of PDM dynamic derivatives of ANF at zero angles of attack at various Mach numbers and studied in detail the effects of various simulation parameters (including numerical simulation parameters, aerodynamic model parameters, and mesh density) on the calculation results of PDM dynamic derivatives.
The rest of the paper is arranged as follows: Section 2 presents some theoretical basic of calculating PDM derivatives.Section 3 introduces the numerical methods used in this paper.The details of ANF geometry and meshes are described in section 4. Section 5 illustrates the test results and associated discussions.Section 6 is the concluding remarks.

Theoretical basic of calculating PDM derivatives
By solving the unsteady RANS equation, we can obtain the aerodynamic forces and moments under the body-axis coordinate system and then use the free flow density, velocity, reference length, and reference area as the reference variables to nondimensionalize these aerodynamic forces and moments to obtain the aerodynamic force coefficients and moment coefficients.Then these data are processed by linear flight mechanics theory to calculate the PDM derivatives.The PDM derivatives is defined as follows [1], [2]: Where q and   are the pitching angular rate and the rate of change of the angle of attack, respectively.mq C is the derivative of the pitching moment coefficient to the dimensionless pitching angular rate, and m C   is the derivative of the pitching moment coefficient to the dimensionless rate of change of the angle of attack.  is the free flow density, V  is the free flow velocity, S is the reference area, D is the reference length.
It should note that since it is usually difficult for experimental and numerical methods to calculate individual components mq C and m C   , the sum of dynamic derivative coefficients of PDM is generally calculated in practical applications, i.e. [ ] In the following contents of this paper, PDM represents the sum of the two dynamic derivative coefficients, which will not be explained again.All forces and moments are defined according to the missile body-axis coordinate system, the x-axis is positive from the head to the tail along the missile body axis, the z-axis is positive upward in the symmetry plane, the y-axis meets the right-handed spiral rule, and the right-handed is positive.
Transient planar pitching motion, also known as forced planar pitching motion or forced oscillation, is the motion whereby harmonically oscillates about the center of gravity of the aircraft in the process Where 0 ( ) m C  is the static pitching moment coefficient at α 0 .m C  is the static derivative of the pitching moment.Note that the unit of attack angle is the radian.At the angle of attack α=α 0 , the PDM derivative is obtained by a forced transient planar pitching motion of sinusoidal oscillation around α 0 .The sine function of the angle of attack is defined as: Where A represents the amplitude of forced oscillation (usually less than 1 o ), ω is the angular velocity, (t)


is the instantaneous angle of attack at time t.For forced planar pitching motion, the pitching angular velocity and the rate of change of the angle of attack are equal, i.e In order to simplify and facilitate, and also to be consistent with the definition of literature, the reduction frequency k is defined as: Now we calculate the PDM derivative at the point of average angular displacement.This method is often applied to the calculation of dynamic derivatives of general missile shapes [4,5,11,14].Fig. 1 is an example of the pitching moment coefficient changing with the angle of attack in the forced pitching oscillation motion of ANF [5].The forced oscillation causes the hysteresis effect of the pitching moment coefficient with the change of the angle of attack, so the loop in the figure is usually called a "hysteresis loop".For most aircraft, this change is periodic, and the hysteresis loop is symmetry about the angle of attack α 0 .

Figure 1. Sample pitching moment history as α0=30 o
From equation (3), when α=α 0 , 0 0 , ( 1 ) , 0 , 1 ,2 ,3 , Since the missile passes through every pitch cycle α 0 twice, so n represents every half cycle.Combine equations ( 2), ( 4), ( 5) and ( 6), omit higher order terms, then get, 0 ( ) According to the hysteresis curve, it is noted that the missile starts to pitch up when t=0, C m is defined as follows: Therefore, the sum of PDM derivatives can be obtained by the Pitch-up part or Pitch-down part of equation (8).Due to the symmetry of pitching oscillation motion, add the two parts of equation ( 8), Then the PDM derivative only needs two points at α 0 during the oscillation period.This calculation is very simple.
In order to simulate the forced sinusoidal motion defined by equation ( 3), it is necessary to determine the values of amplitude A and reduced frequency k.According to the research results of ref [3], we select A=0.25 o and k=0.1 as the benchmark parameters in this paper.In the research scope of this paper, PDM is independent of the selection of A and k values.According to equation( 9), the pitching moment C m and C m change linearly with k and A.
The period T and frequency f are defined as: Because the planar pitching simulation is time-accurate calculation, it is necessary to determine the time step of numerical simulation.Assuming that each pitch cycle is iterated in N steps, the time step is, Research in literature [3] shows that N ≥ 100 can meet the accurate pitch damping prediction of ANF missiles.Here we select N=200 iterations per pitch period as the benchmark parameter.

Numerical method
In this study, NNW-FlowStar [15] (FlowStar for short) is used to numerically solve the threedimensional compressible unsteady RANS equation to simulate the flow field of the missile, calculate the aerodynamic force and moment of the missile, and then solve the static derivative and pitching damping dynamic derivative of the missile.Double-precision format is used in the calculation.
FlowStar is an industrial CFD software developed by China Aerodynamics Research & Development Center (CARDC) based on arbitrary unstructured hybrid grid solution technology and large-scale parallel computing technology.The software applies to aircraft, helicopters, missiles, and other aircraft with subsonic, transonic, supersonic at high angle of attack steady and unsteady aerodynamic calculation, weapon delivery trajectory calculation, inlet internal flow calculation, jet flow field calculation, etc.
For the steady state calculation, the implicit LUSGS method is used for time discretization, and the local time step and multigrid technology are used to accelerate the convergence of solution.The spatial discrete scheme adopts the second-order upwind Roe scheme, the gradient calculation adopts the nodal Gauss method, and the flux limiter uses the Venkatakrishnan limiter.The SST two-equation model is used for the turbulence model, and the SA one-equation model is used to calculate some states for comparison.Unsteady RANS calculation adopts the dual-time stepping method with the outer integration using the physical global time step and the inner iteration using the pseudo-time step subiteration.The internal iteration is solved by the same method as the steady calculation.Generally, if the sub-iteration number is taken as 40, the residual of internal iteration can be reduced by 3 to 4 orders of magnitude, which is sufficient to meet the requirements of convergence.
The transient planar pitching motion method is implemented firstly by obtaining the hysteretic aerodynamic force directly related to the vibration frequency under the condition of periodic motion by solving the unsteady fluid dynamic equation of the moving grid system and then obtaining the aerodynamic derivative.The principle of the transient planar pitching motion method is completely consistent with the forced oscillation test.The specific calculation method is described as follows.Firstly, we calculate the solution of steady flow field at a given Mach number and an average angle of attack α=α 0 , which is taken as the initial condition for unsteady planar pitching motion simulation.The transient planar pitching motion is defined by the sine function of equation ( 3).The initial transient change of unsteady calculation usually only affects the oscillation of the first period.When N=200 iterations are selected for each pitching period, three pitching periods are usually enough to obtain the convergence hysteresis loop of the aerodynamic moment coefficient.Using the m C history of the last cycle and using equation ( 8), PDM can be obtained.
The viscous non-slip adiabatic wall condition is adopted for the wall boundary, and the characteristic boundary condition based on local one-dimensional Riemann invariants is adopted for the far field boundary condition.
The calculation parameters are as follows: freestream pressure P=101325.0Pa,temperature T=293.15K,Mach number is 0.5~4.5, angle of attack is 0 o , and Reynolds number based on body diameter is 0.342e6~3.08e6.

ANF's Geometry and meshes
The test case is the ANF missile, which is a standard missile model.Fig. 2 shows the geometric details of the model [5].The ANF standard model has a diameter of 30.0mm, a 10° conical head of 2.84 times caliber, a column body of 7.16 times caliber, and four symmetrical tails.The tail is a 1x1 square tail with a sharp leading edge and a bottom thickness of 0.08 times caliber.There is no deflection angle for the tail in this calculation.The center of gravity is 5.5 times the caliber of the nose.The reference length of aerodynamic force and moment coefficient is 0.03m, the reference area is 7.06858347e-4m 2 , and the moment reference point is the centroid position (x=0.165m,y=0.0m, z=0.0m).

Figure 2. Computational geometry
According to the geometric parameters in Fig. 2, we generated the CAD digital model of ANF and generated the unstructured hybrid mesh containing prism, pyramid, and tetrahedron elements of CGNS format.
Taking Ma=4.5, Re=3.08e6 (based on the reference length of 0.03m) as the standard, based on y + =1, the spacing of the first layer of the grid in the normal direction is determined as d=0.00075mm, and the boundary layer grid grows with a ratio of 1.2.The far field is taken as 20 times the length of the projectile.The grid quantity is 17.27 million, including 9.92 million tetrahedral cells, 7.29 million prism cells, and 60 thousand pyramid cells, as shown in Fig. 3.

Figure 3. The computational grid of ANF missile model
To study the influence of the grid resolution on the results, four sets of grids with different sizes (the coarsest, coarse, medium, and fine grid) are generated based on the above grid (named the medium grid).The scale ratio of each set of grids in the three coordinate directions of x, y, and z is 2 1/3 , so the total grid size between two adjacent sets of grids differs by about 2 times.Based on the domain decomposition strategy, the generated grid is divided into different sub-domains and loaded into different processors for parallel computing through pre-processing.But when the grid size is too large, serial preprocessing is impractical due to the memory bottleneck and is too timeconsuming.A parallel preprocessing technique for unstructured-grid is used here to deal with this problem [20].As shown in Fig. 4 all the preprocessing flow is in parallel.

Results and discussion
Section 5.1 briefly introduces the computational results of the steady state, and analyzes the influence of the grid resolution on the calculation results.Section 5.2 and section 5.3 study the different simulation parameters on the precision of dynamic derivatives simulation.Section 5.4 compares the CFD results and the free flight test data under different Mach numbers to further verify the proposed method.

Steady state result
Steady-state simulation can be used to test the accuracy of the CFD solver, and the results can be used as the initial flow field for unsteady calculations of dynamic derivatives.Fig. 5 shows the typical subsonic and supersonic flowfield ( contour of Mach) at a steady state.From the figure, we can clearly see the shock wave near the nose and the wing, and the low-speed flow at the bottom.Fig. 6 compares the calculation results and test results of the steady state at zero angles of attack when the Mach number is between 0.5 and 4.5, including the axial force coefficient C x , the normal force static derivative C Nα , and the static derivative of pitching moment C mα .Since the aerodynamic coefficient changes linearly with the angle of attack near zero, C Nα and C mα are obtained using the result of the angle of attack 0 o and 1 o by linear fitting.
It can be seen from Fig. 6 that the calculated results agree well with the test results within the calculated Mach number range.In the case of supersonic, compared with the wind tunnel test, the calculated results are in better agreement with the free flight test results.
At the same time, the calculation results of four sets of grids with different sizes are compared when the Mach number is 0.5, 0.9, 1.1, 2.5, and 4.5.It can be seen from Fig. 6 that under the subsonic and transonic conditions, the grid resolution has a certain influence on C x .but little influence for C xα and C mα .At Ma=0.5, the difference between the C x calculated by the coarsest grid and the fine grid is 4.4%.At supersonic, the results of the four sets of grids are basically consistent.In summary, the medium grid that was selected as the reference grid meets the requirements of grid independence.
A comparison of the difference between the SA turbulence and the SST turbulence model is also made in Fig. 6.Compared to the SST model, the C x calculated by the SA model results is significantly higher at subsonic and transonic, about 19.6% at Ma=0.5.While for C Nα and C mα , the difference between the two models is not obvious.At supersonic, all the results calculated by the two turbulence models are basically the same.Therefore, we use the medium grid and SST turbulence model for the following calculations.

Impact analysis of numerical parameters
As mentioned earlier, since the planar pitching motion is an unsteady phenomenon, the dual-time stepping method is used for calculation.Physical time step ΔT is calculated based on N=200 steps for each selected period, and the sub-iteration number (represented by i) is based on 40.Fig. 7 ~ Fig. 10 show the influence of changing i and N on the PDM derivative respectively.Here, two incoming flow Mach numbers 0.9 and 4.5 are computed, and the amplitude and reduction frequency taken as the reference values A=0.25 o , k=0.1.When i changes from 10 to 80, it can be seen from Fig. 7 that for these two Mach numbers, whether N=100 or N=200, except for i=10, PDM is a gradual convergence process with the increase of i.When i=40 or i=50, the difference between PDM and the result of i=80 is less than 2%, and when N=200, the difference is smaller.Comprehensively considering the calculation accuracy and efficiency, i=40 is selected as the optimal parameter for calculation.Fig. 7 shows the impact of inner iteration numbers on the hysteresis loop at Ma=4.5.It can be seen that all the hysteresis loops are very closed except i=10.When i increases gradually, the hysteresis loops gradually converge to a fixed loop.When i=40, N changes from 100 to 500, as shown in Fig. 8.For these two Mach numbers, PDM is a process of gradual convergence.When N increases from 200 to 500, PDM only increases by 0.5% for these two Mach numbers.Considering the calculation efficiency, when N is taken as 200, a good calculation efficiency can be obtained at the expense of small calculation accuracy.Therefore, N=200 is selected as the optimal parameter for calculation in this paper.Fig. 9 shows the hysteresis loop comparison of different N values with Ma=4.5.It can be seen that the hysteresis loop changes very little when N changes from 100 to 500.From the enlarged figure, with N increases, the hysteresis loop will gradually converge into a fixed loop.The difference between N=200 and N=500 can be ignored.

Impact analysis of physical parameters
To study the impact of amplitude A, the calculation results of PDM are compared when A is taken as 0.125 o , 0.25 o , 0.5 o and 1 o under the condition of Ma=0.9 and Ma=4.5 respectively (Fig. 11).In order to study the impact of reduced frequency k, the calculation results of PDM are also compared when k is 0.025, 0.05, 0.1, and 0.2 at Ma=0.9 and Ma=4.5 respectively (Fig. 12).
In the above two figures, the pitch moment difference is the ordinate, and A and k are the abscissa.According to equation ( 9), for the same flight attitude, i.e. the same PDM, (C m+ -C m-) is proportional to A and k.The curve is very close to a straight line.'r' in the figure is a linear correlation coefficient, and the value is very close to 1.This indicates that the PDM does not change with amplitude A and the reduction frequency k within the range tested in this paper.Therefore, the standard parameters A=0.25 o , k=0.1 selected in this paper meet the requirements for solving the PDM.  13 shows the effect of amplitude A change on pitch moment when Ma=4.5.With an angle of attack α as the X-axis, as the amplitude A increases, the hysteresis loop curve also increases, and the center of the hysteresis loop is located at(α, C m )=(0.0,0.0).With time t as the X-axis, as the increase of amplitude A, the amplitude of the sine wave increases, but the frequency does not change.It can be seen that the initial start of the planar pitching motion basically only affects the first quarter of the period and then becomes the periodic state.Fig. 14 shows the effect of reduction frequency k change on pitching moment when Ma=4.5.With an angle of attack α as the X-axis, as the increase of the reduction frequency k, the eccentricity of the hysteresis loop decreases, and the center of the hysteresis loop is located at(α, C m )=(0.0,0.0).With time t as the X-axis, as the increase of amplitude A, the amplitude of the sine wave and the frequency will also increase.It is found again that the initial start of the planar pitching motion basically only affects the first half of the period, and then becomes periodic.
We can conclude that the reference parameters A=0.25 o and k=0.1 selected in this paper provide sufficient disturbance for the calculation of pitch damping derivatives and meet the linear assumptions of the transient planar pitching method.It can be seen that in the subsonic region, there is a very big difference between the CFD result and the test data.At Ma=0.75, the calculation result is 101% larger than the test result.In the transonic region of 0.9<Ma<1.3, the dispersion of the test results is large, and the calculated results are approximately near the average of the test results.According to the analysis of document [3], it is believed that the free flight test data used in this paper in the subsonic and transonic region is not reliable and needs more test data to confirm.In the supersonic region (Ma>1.3), the calculated results are in good agreement with the test results.
'CFD_document' in the figure is the calculation result of the same dynamic derivative calculation method used in ref [3].The two CFD calculation results are in good agreement, which further verifies the correctness of the calculation results in this paper.

Conclusion
In this paper, the pitching damping moment of the ANF standard missile model in the subsonic and transonic range was calculated based on the transient planar pitching method.The influence of different numerical parameters and physical parameters on the calculation accuracy of dynamic derivatives was studied.Comparisons of the numerical results with the free flight test results and the data from the literature were made.The main conclusions are as follows: In the transient planar pitching method, the dynamic derivative of pitching damping is related to the pitching moment in the forced sinusoidal oscillation motion.In this paper, the pitching moment changes linearly with amplitude (0.125 o ≤ A ≤ 1 o ) and reduced frequency (0.025 ≤ k ≤ 0.2), so the prediction of pitching damping derivative is independent of amplitude and reduced frequency.The range of this linear change depends on the geometric shape and flight conditions of the aircraft, such as Mach number, flight altitude, etc.
The calculation accuracy of the dynamic derivative of pitch damping is related to the physical time step and the iteration number of the dual-time stepping method in RANS calculation.The smaller the time step or the more internal iteration steps, the more accurate the result, but it will increase the amount of calculation and affect the efficiency.On the premise of meeting the engineering requirements, accuracy can be appropriately sacrificed to improve calculation efficiency.After a comprehensive analysis, in this paper, the calculation steps of each sinusoidal oscillation period N=200 and the subiteration number i=40 meet the calculation accuracy requirements.
Compared with the free flight test results, the calculated results are in good agreement in the supersonic region, but there is some difference in the subsonic and transonic ranges.According to the analysis of the reference paper, in the subsonic and transonic range, the accuracy of the free-flight test data used in this paper maybe be uncertain, which needs more test results to confirm.
The research results of this paper can provide significate guidance in the prediction of missile pitch damping derivatives by unstructured CFD solver, including the selection of turbulence models, grids, numerical parameters, and physical parameters.Further work will concentrate on the confirmation work under the condition of a large angle of attack and research dynamic derivatives of aircraft shape.

Figure 4 .
Figure 4. Parallel preprocessing flow of unstructured grid

Figure 9 .Figure 10 .
Figure 9. Impact of global iterations per cycle on PDM

Figure 11 .Figure 12 .
Figure 11.Effect of amplitude A on pitch moment difference

Figure 13 .Figure 14 .
Figure 13.Effect of amplitude A on pitch moment

Figure 15 .
Figure 15.Comparison of CFD and experiment on PDM