AL-Khwarizmi Engineering Journal, 4, 2, (2009) Al-Khwarizmi Engineering Journal Al-Khwarizmi Engineering Journal, Vol. 5, No. 1, PP 42-52 (2009) Numerical Computations of Transonic Critical Aerodynamic Behavior of a Realistic Artillery Projectile Ahmed F. M. Kridi Department of Manufacturing Operation Eng. /Al–Khwarizmi College of Eng. /University of Baghdad (Received 31 January 2008; accepted 25 March 2009) Abstract The determination of aerodynamic coefficients by shell designers is a critical step in the development of any projectile design. Of particular interest is the determination of the aerodynamic coefficients at transonic speeds. It is in this speed regime that the critical aerodynamic behavior occurs and a rapid change in the aerodynamic coefficients is observed. Two-dimensional, transonic, flow field computations over projectiles have been made using Euler equations which were used for solution with no special treatment required. In this work a solution algorithm is based on finite difference MacCormack’s technique for solving mixed subsonic-supersonic flow problem. Details of the asymmetrically located shock waves on the projectiles have been determined. Computed surface pressures have been compared with experimental data and are found to be in good agreement. The pitching moment coefficient, determined from the computed flow fields, shows the critical aerodynamic behavior observed in free flights. Keywords: CFD, Euler Equation, Artillery Projectile, MacCormack’s technique. 1. Introduction The flight of projectiles covers a wide range of speeds. The accurate prediction of projectile aerodynamics at these speeds is of significant importance in the early design stage of a projectile. The critical aerodynamic behavior occurs in the transonic speed regime, 0 .9< M < 1.2 where the aerodynamic coefficients have been found to change by as much as 100%. Of particular interest is the determination of the pitching-moment coefficient since it is used to determine the static stability of the projectile. The critical behavior in this case is usually characterized by a rapid increase in the coefficient followed by a sharp drop. This rapid change in the pitching moment coefficient can be attributed in part to the complex flow structure and, in particular, to the asymmetrically located shock waves, which exist on projectiles flying at transonic speeds at angle of attack. Computations of two-dimensional flow fields at transonic speeds are thus needed to predict the critical aerodynamic behavior. A considerable research effort has been focused on the development of modern predictive capabilities for determining projectile aerodynamic 1-5 . Numerical capabilities have been developed primarily using Euler equation computational techniques and have been used to compute flow over slender bodies of revolution at transonic speeds. Flow field computations have included both axisymmetric 4 and two-dimensional situations. 1, 2, 3, 5 initial computations 1-3 did not include the wake or base region of a projectile and, thus, ignored the upstream effect of the baser region flow on the afterbody. Benek et. el. 6 show the development of a chimera grid scheme. This scheme provides multiple regions where communications between grids are done by interpolating in regions of overlap. A blocked grid approach reported by Belk and Whitfield 7 does not require interpolations at the interfaces and has been successfully used to obtain Euler solutions over a wing. The scheme used by Benek et. el. 6 is generally complicated since it allows for embedding a block or zone into another. Recently, a simple composite grid scheme has been Ahmed F. M. Kridi Al-Khwarizmi Engineering Journal, Vol. 5, No. 1, PP 42-52 (2009) 43 developed here a large single grid was partitioned into smaller grids. Each of the smaller problems was solved separately with simple data transfers at the interfaces. The initial results obtained were very promising. The present effort extends the use of this simple algebraic grid scheme to include the correct modeling of a projectile. The aim of the present work is to study numerically the two-dimensional projectile flow field for air missiles under the conditions of free stream Mach number 0.9, 0.92, 0.94, 0.96, 0.98, 1.0, 1.1, 1.2 in transonic speed case at atmospheric conditions, pressure 101325 N/m 2 and temperature 15 o C. 2. Theoretical Analysis The calculation of the projectile flow field is of considerable importance to the efficient design of projectile. These flow fields are very complex due to their mixed hyperbolic-elliptic nature, the influence the forbody and viscous effects, as well as the three-dimensionality. The complexities of three-dimensional viscous inlet flow make their numerical prediction a very difficult task; therefore, the calculation of two dimensional inlets is an step toward that direction. The projectile flow fields calculated by a two- dimensional computational method, the problem of employing an explicit, time-marching, finite difference procedure to solve the Euler equation formulated in body-fitted coordinates. The method can be used for a flow field in both supersonic and subsonic regions. 2.1. Model and Computational Grids The model used for the computational study presented here is an idealization of a realistic artillery projectile geometry. The experimental model shown in Fig. 1 is a secant-ogive cylinder- boattail (SOCBT) projectile. It consists of a three caliber (one-caliber=maximum body diameter), sharp, secant-ogive nose, a two-caliber, cylindrical midsection, and One-caliber boattail. The computational grid used for this computation is shown in Fig. 2 shows the longitudinal cross section of the two-dimensional grid. The algebraic equation is used to relate the grid points in the computational domain to those of the physical domain. This objective is met by using an interpolation scheme between the specified boundary grid points to generate the interior grid points. Clearly, many algebraic equations can be introduced for this purpose. 2.2. Governing equations For high Reynolds number flows, viscous effects are confined to the vicinity of the surface, where large velocity gradients exist. This region is known as the boundary layer. Outside the boundary layer, the velocity gradients are negligible resulting in zero shear stresses. This region is called the inviscid region, and solution procedures for the inviscid flow region are governed by the Euler equations and the solution of this research depends on it, which is written in conservation-law form for two-dimensional flows of a perfect gas 8 . The general compact vector form is given as:- 0         y F x E t U , Where:                                                )( , )( , 2 2 pe p u F peu uv pu u E e v u U ttt            …(1) Fig.1. Model Geometry of SOCBT Projectile. Calibe r Cylinder Boattail Fig.2. Two Dimension Projectile (General Shape). m m Lower outer flow boundary Upper outer flow boundary Out flow Out flow x y Ahmed F. M. Kridi Al-Khwarizmi Engineering Journal, Vol. 5, No. 1, PP 42-52 (2009) 44 u and v are the velocities along the x and y coordinates, respectively, p is the pressure,  is the density, and et is the total energy per unit volume. And U, E, F, are the fluxes vectors. To transform the Euler Equation (1) into curvilinear coordinates (,), an independent variable may be written as follows:-          FE t U …(2) Where:                1 2222 4 13 12 1 ) 22 () 22 ( U vu e vu e J U vU J v U uU J u U J U J U U     …(2a)                     v J p Uu J p UE J p vUvUE uU J p uUE UUE FE J E yx yx yx yx yx ).().( ).().( ).().( )( 1 444 323 322 321      …(2b)                     v J p Uu J p UF J p vUvUF uU J p uUF UUF FE J F yx yx yx yx yx ).().( ).().( ).().( )( 1 444 323 322 321      …(2c) where: , u, v, p and e are a (primitive variables) non dimensional density, velocity in x-direction, velocity in y-direction, pressure and internal energy respectively. 2.3. MacCormack's Technique: The MacCormack's time marching method is an explicit finite-difference technique. It is second-order-accurate in both space and time. This method will be used to solve the Euler Equation itemized in Equations (2a) to (2c) with march in time to steady state solution by solving the flow properties at every (i,j) spatial location, assuming that the flow field at each node is known at time t. Consider the U flow field variable at grid point (i,j) at time t+t. In MacCormack's method, this is obtained from t t U av t ji tt ji UU           ,, …(3) Where (U/t)av is a represented mean value of (U/t) between time t and t+t. The value of (U/t)av is calculated as a second order accuracy, and once again, U is a flow field variable known at time t. Either from initial condition or as a result from the previous iteration in time. (U/t)av is defined as: av tt ji c t jiav t U t U t U                                  ,, 2 1 …(4) To obtain a value of (U/t)av , produce initial fluxes from primitive variables by using Equations (2a) to (2c) and then there are two major steps taken as:- 1. Predictor step: (U/t) t i,j is calculated using forward spatial difference on the right side of the governing equation (2) from the known flow field at time t. The predicted value of the flow field variable can be obtained at t+t, as follows:- )()( ,1,,,1, , FFEEUU t ji t ji t ji t ji t ji tt ji tt p          ...(5) For the interior nodes a new value of parameters (, u, v, p and e) will be found from the new fluxes using equation (2a). After that, the updating boundaries must be done. From the new value of parameters that have been derived, the influence of the boundary updating flux must be done using equations (2a) through (2c). 2. Corrector step: Using backward spatial differences, the predicted value (from step 1) is inserted into the governing equation such that a Ahmed F. M. Kridi Al-Khwarizmi Engineering Journal, Vol. 5, No. 1, PP 42-52 (2009) 45 predicted time derivative tt jip tU   , )/( can be obtained. The equation of backward space is illustrated as:- )()( 1,,,1,,, FFEEUU tt P tt P tt P tt P tt P tt C jijijijijiji tt           ...(6) Then, substitute tt jic tU   , )/( and (U/t) t i,j by equation (4) to obtain the average value, to find the corrected second order accurate value of U at time t+t, combining equations (4) and (6) and substituting them by equation (3) yields:-                    )()( 2 1 1,,,1,, ,, pp t pp t p FFEEUUU tt ji tt ji tt ji tt ji tt ji t ji tt ji  …(7) And from the new corrected fluxes it is possible to obtain the correct value of parameters (, u, v, p and e) for all interior nodes by using equation (2a) then updating boundaries. The above steps are repeated until the flow field variable approaches a steady state value; this is the desired steady state solution. 2.4. Time Step Calculation: The value of t cannot be arbitrary, rather it must be less than some maximum values for stability, it was stated that t must obey the Courant-Friedriches-Lowry criterion CFL. The CFL criterion states that physically the explicit time step must be not greater than the time required for a sound wave to propagate from one grid to next. The maximum allowable value of CFL factor for stability in explicitly time dependent finite difference calculation can vary from approximately 0.5 to 0.1. To determine the value of time step, the following version of the CFL criterion 9 is used. Where ai,j is the local speed of sound in meters per second, and C is a constant. 1 22ji, ji,ji, ji,CFL Δη 1 Δξ 1 a Δη v Δξ u )(          t …(8) And, t = min [C (t CFL) i ,j] . 2.5. Boundary Condition: The Euler equation has an unlimited number of solutions. What makes a solution unique is the proper specification of initial and boundary conditions for a given PDE (Euler equation). A set of boundary conditions must be specified, it referred to as the “analytical boundary condition” Once the PDE is approximated by a FDE, Thus the FDE will require additional boundary conditions. This boundary condition will be referred to as “numerical boundary condition”. As for the problem under consideration, there are four types of boundaries: solid, inflow, outer and outflow. 2.5.1. Solid boundary Condition For the two solid boundary conditions (projectile upper surface and lower surface), the tangency grid body surface must be satisfied for inviscid flow. The components of the momentum equation for the two-dimension flow may be expressed with some mathematical steps, as 10 :-            0) J Pη (η) J Pξ (η) J Pη (η) J pξ (η ) J Vρu (η) J Vρv (η) J Uρu (η) J Uρu (η η y yξ y yη x xξ x x ηxηxξyξx …(9) a. Upper projectile surface:- A finite difference equation for the upper equation is obtained, as a second order central difference approximation for the  derivatives and a second order forward difference approximation for  derivatives are used. b. Lower projectile surface:- A second order central difference approximation for  derivatives and second-order backward difference approximation for  derivatives are used. 2.5.2. Outer Flow Boundary The upper outer flow boundary is the air flow out from the numerical simulation of two- dimensional projectile at 1.4 meter in the x- direction. To calculate the properties at this boundary, first order backward transformation derivatives are used. The lower outer flow boundary is the air flow out from the numerical simulation of two-dimensional projectile at zero meter in the x-direction, the first order forward transformation derivatives are used to calculate the properties as shown in Fig. 2. Ahmed F. M. Kridi Al-Khwarizmi Engineering Journal, Vol. 5, No. 1, PP 42-52 (2009) 46 2.5.3. Out flow Boundary The outflow boundary illustrated in Fig. 2 represents the airflow above and below the projectile. This airflow is got out from the numerical simulation and is two meters far from the original point. The air flow properties are calculated by using backward transformation derivatives. 3. Results and Discussions The projectile characteristics at transonic speed are dominated by the shock-wave systems that go into their design. In the following results we put aside temporarily the problems of boundary layer and flow separation and consider the simple nature and properties of shock wave. The implicit time marching procedure was used to obtain the desired steady-state result. Initial conditions were free-stream everywhere, and the boundary conditions were up-dated explicitly at each time step. For (SOCBT) Projectile, 0.9