HB-2彈體的空氣動力及熱傳模擬

全文

(1)國立交通大學機. 械工程學碩士論文. 系. HB-2 彈體的空氣動力及熱傳模擬. Aerodynamics and Heat Transfer Simulation of HB-2 Model. 研究生：柳志良指導教授：吳宗信. 博士. 中華民國九十七年七月.

(2)

(3) HB-2 彈體的空氣動力及熱傳模擬. Aerodynamics and Heat Transfer Simulation of HB-2 Model 研究生：柳志良. Student：Chih-Liang Liu. 指導教授：吳宗信博士. Advisor：Dr. Jong-Shinn Wu. 國立交通大學機械工程學系碩士論文 A Thesis Submitted to Institute of Mechanical Engineering College of Engineering National Chiao Tung University in Partial Fulfillment of the Requirements for the Degree of Master of Science In Mechanical Engineering July 2008 Hisnchu, Taiwan 2008 年七月.

(4) HB-2 彈體的空氣動力及熱傳模擬學生：柳志良. 指導教授：吳宗信國立交通大學機械工程學系. 摘要. 火箭飛行歷經次音速、穿音速、超音速至極超音速等不同的飛行階段。同時隨著高度的增加，大氣密度迅速降低，火箭也同時飛越連續體流域、過渡流域至稀薄流域。對火箭氣動力設計系統來說，克服各種不同的飛行環境所面對的氣動力問題是很重要的。本文研究重點為使用 UNIC-UNS 模擬 HB-2 彈體飛行在不同馬赫數和攻角下的空氣動力和熱傳情形。首先我們藉由網格測試去選出適當的網格大小，再使用該網格去模擬 M = 3.01 和 M = 5.10 時不同攻角的氣動力情形和 M = 9.59 時的熱傳情形。在氣動力係數方面， M = 3.01， Re = 2.2 E 06 ，α = 0 0 ~ 15 0 時的氣動力結果： Ca = 0.7 ， Cn = 0 ~ 1.6 ， Cm = 0 ~ −1.4 ， Xcp / L = 0.5 ~ 0.56 。 M = 5.10 ， Re = 2.5 E 06 ， α = 0 0 ~ 15 0 時的氣動力結果： Ca = 0.7 ~ 0.9 ， Cn = 0 ~ 1.4 ， Cm = 0 ~ −1.4 ， Xcp / L = 0.52 ~ 0.56 。在熱傳情形方面，我們使用幾種不同的. 壁面溫度和流場模式去模擬，最終是較接近室溫的層流模式的熱傳分佈情形最為接近實驗的結果。而所得到結果的物理現象也都符合我們預估的。將結果和實驗的數據作比較，得出大部份的氣動力參數誤差都在 4% i.

(5) 以下，變化的趨勢也都和實驗的結果很吻合。這些結果都可以提供為後續研究人員的參考資料。. ii.

(6) Aerodynamics and Heat Transfer Simulation of HB-2 Model Student: C. L. Liu. Advisor: Dr. J. S. Wu. Department of Mechanical Engineering National Chiao-Tung University. Abstract. The rocket flies through various stages such as subsonic, sonic, supersonic, and hypersonic. Atmospheric density reduced rapidly with increasing of height, the rocket flies over the continuous flow, transitional flow, and rare flow at the same time too. To the aerodynamics force design system of the rocket, it is very important to overcome the aerodynamics problem that some kinds of flight environment cause. In this thesis, we apply a parallelized Navier-Stokes equation solver, named UNIC-UNS, to do the heat transfer and aerodynamics simulation of HB-2 model at different Mach numbers and attack angles of the flight. We make the grid convergence first to choose appropriate quality of grids. We do the aerodynamics simulation with the cases ( M = 3.01 , M = 5.10 ) with different angles of attack and the heat transfer simulation with the cases ( M = 9.59 , Re = 1.97 E 05 and Re = 1.87 E 05 ) with the grid file. The results of aerodynamics coefficients at M = 3.01 and Re = 2.2 E 06 are Ca = 0.7 , Cn = 0 ~ 1.6 , Cm = 0 ~ −1.4 , and Xcp / L = 0.5 ~ 0.56 .. The results of aerodynamics coefficients at M = 5.10 and Re = 2.5E 06 are Ca = 0.7 ~ 0.9 , Cn = 0 ~ 1.4 , Cm = 0 ~ −1.4 , and Xcp / L = 0.52 ~ 0.56 . In heat. iii.

(7) transfer simulation, we use some kinds of wall temperature and flow model. The results of comparison of non-dimensional heat flux distributions at Tw = 353.6k and laminar flow are more accuracy then other cases. We get the. results with physical phenomena as we know. We compare the results with experimental data. The error of most aerodynamics coefficients that we get is less than 4%. The tendency of coefficient’s change is also very similar to experimental data. The results and mesh can be reference material to supply the people of follow-up study.. iv.

(8) 誌謝. 首先感謝我的家人，支持我在這條不順遂的求學道路上行走下去，現在才得以完成碩士班的學業。在碩士班的兩年中，特別感謝吳宗信老師不辭辛勞的教誨，除了在研究上，也教導我許多身為台灣人應有的自覺；另外，也要特別感謝黃俊誠老師和陳彥升博士，感謝他們在百忙之中仍抽空幫我解答一些學業上的疑惑。感謝鄭凱文、李允民、洪捷粲、邱老大、江明鴻、胡孟樺、林雅茹、林昆模、盧勁全、林宗漢、謝昇汎、洪維呈等諸位學長姐，感謝他們不厭其煩的為我在研究以及其他方面解惑。感謝蘇正勤、吳玟琪、劉育宗、鄭丞志、呂政霖諸位同學，感謝有你們和我一起在研究上互相勉勵和進步。感謝必任、穎志、逸民、丁丁、諸位學弟以及于恩，感謝你們在各方面的協助。感謝 Hadley、Aziz、Matt，感謝你們使我的英文溝通能力和國際觀有所提升。感謝 APPL 實驗室，感謝這個巨大的寶庫，讓我能在研究所兩年的時間裡任意拿取各種學識上的寶物。最後，感謝台灣這個國家，感謝這片孕育我成長茁壯的土地，願我在將來的數十年也能和你繼續成長。. 柳志良謹誌 2008 年 8 月于台灣風城. v.

(9) Table of Contents 摘要…………………………………………………………………………...………………………………..………….i Abstract………………………………………………………………………………………………………...……….iii 誌謝…………………………………………………………………………………….…………………………………v Table of Contents .................................................................................................................. vi List of Tables....................................................................................................................... viii List of Figures....................................................................................................................... ix Nomenclature........................................................................................................................ xi Chapter 1 Introduction .........................................................................................................1 1.1 Motivation and Background .....................................................................................1 1.2 Literatures Survey.....................................................................................................2 1.3 Specific Objectives of the Thesis..............................................................................3 Chapter 2 Numerical Method ..............................................................................................4 2.1 Governing Equations ................................................................................................4 2.2 Spatial Discretization................................................................................................4 2.3 Time Integration........................................................................................................6 2.4 Pressure-Velocity-Density Coupling.........................................................................7 2.5 Linear Matrix Solver.................................................................................................8 2.6 Parallelization ...........................................................................................................9 Chapter 3 Results and Discussions ....................................................................................10 3.1 Overview.................................................................................................................10 3.2 Grid Convergence Test............................................................................................10 3.2.1 Grid Configuration.......................................................................................10 3.2.2 Flow Conditions and Simulation Conditions...............................................11 3.2.3 Validation .....................................................................................................12 3.2.3.1 Density, Pressure and Mach Number Distributions..........................12 3.2.3.2 Coefficients of Normal Force ...........................................................13 3.2.3.3 Coefficients of Pitching Moment......................................................13 3.3 Aerodynamics Simulation with Different Mach numbers and Angles of Attack 14 3.3.1 Results in Different Mach Numbers and Attack Angles..............................14 3.3.1.1 Flow Conditions and Simulation Conditions....................................14 3.3.1.2 Density, Pressure and Mach Number Distributions..........................15 3.3.1.3 Coefficients of Axial-Force...............................................................16 3.3.1.4 Coefficients of Normal-Force ...........................................................17 3.3.1.5 Coefficients of Pitching-Moment .....................................................17 3.3.1.6 Normal-Force Curve Slope...............................................................18 3.3.1.7 Center of Pressure Location..............................................................18 vi.

(10) 3.3.1.8 Comparison of Non-Dimensional Pressure Distributions ................19 3.4 Heat Transfer Simulation ........................................................................................20 3.4.1 Flow Conditions and Simulation Conditions...............................................20 3.4.2 Temperature Distributions ...........................................................................21 3.4.3 Comparison of Non-Dimensional Heat Flux Distributions .........................21 Chapter 4 Conclusions .......................................................................................................23 4.1 Summary.................................................................................................................23 References ………………………………………………………………………………….…………………………24. vii.

(11) List of Tables Table.I Flow conditions for grid convergence test............................................................25 Table.II Simulation conditions for grid convergence test.................................................25 Table.III Flow conditions for simulation with different Mach numbers and angles of attack...................................................................................................................25 Table.IV Simulation conditions for simulation with different Mach numbers and angles of attack...............................................................................................................26 Table.V Flow conditions for heat transfer simulation.......................................................26 Table.VI Simulation conditions for heat transfer simulation............................................26. viii.

(12) List of Figures Fig.1.1 The flying diagram of Formosat-2 ........................................................................27 Fig.1.2 Enlarged overset grid around the nodes of the objects cut view at X = 0.22 and symmetric boundary............................................................................................27 Fig.1.3. Intergrid boundaries of the booster subgrid: angles of attack of the airplane at 2o and the booster at 0o and ΔX = 0 : ΔZ = 0.4 , ΔZ = 1.0 , and ΔZ = 3.0 . ...........28. Fig.1.4. Comparison between computed Mach contours and schlieren photographs of supersonic airplane/rocket booster separation at M ∞ = 2.5 , angles of attack of the airplane at 2o and the booster at 0o, and ΔX = 0 : ΔZ = 0.4 , ΔZ = 2.4 , and ΔZ = 5.0 m. .........................................................................................................28. Fig.1.5 Fig.1.6. Standard hypervelocity ballistic correlation configurations..................................29 Generic supersonic transport configuration SYN87-MB Grid Structure: 180 Blocks .................................................................................................................29. Fig.1.7. Generic supersonic transport configuration SYN87-MB solution pressure coefficient (Mach=2.2 α = 3.15 0 C L = 0.105 ) .................................................30. Fig.2.1 Fig.3.1 Fig.3.2 Fig.3.3. Unstructured control volume.................................................................................30 Typical grid distribution of aerodynamics simulation of HB-2 model. ................32 Typical time history of residuals of mass, momentum and energy equations.......32 The distributions of (a) density; (b) pressure; (c) Mach number at various slices of the computational domain. .............................................................................35. Fig.3.4 Fig.3.5 Fig.3.6. The diagram of normal-force and axial-force. ......................................................36 The comparison of normal-force coefficients of grid convergence test................36 The comparison of pitching-moment coefficients of grid convergence test. ........37. The distributions of (a) density at α = 0 0 ; (b) density at α = 15 0 ; (c) 0 0 0 pressure at α = 0 ; (d) pressure at α = 15 ; (e) Mach number at α = 0 ; (f) Mach number at α = 15 0 at various slices of the computational domain...........43 Fig.3.8 The comparison of axial-force coefficients at (a) M = 3.01 ; (b) M = 5.10 ..........44 Fig.3.9 The comparison of normal-force coefficients at (a) M = 3.01 ; (b) M = 5.10 .......45 Fig.3.10 The comparison of pitching-moment coefficients at (a) M = 3.01 ; (b) M = 5.10 .. .....................................................................................................46 Fig.3.7. Fig.3.11. The comparison of normal-force curve slope at different Mach number and Reynolds number. ...............................................................................................47. Fig.3.12 The center of pressure location at different angles of attack...............................47 Fig.3.13 The comparison of non-dimensional pressure distributions at (a) M = 3.01 ,. α = 5 0 ; (b) M = 3.01 , α = 15 0 ; (c) M = 5.10 , α = 5 0 ; (d) M = 5.10 , α = 15 0 ....49 Fig.3.14 The 3-D pressure distribution at (a) M = 3.01 , α = 15 0 ; (b) M = 5.10 , α = 15 0 … ...........................................................................................................50 ix.

(13) Fig.3.15 The distributions of temperature at (a) Tw = 53.6k ; (b) Tw = 353.6k .................52 Fig.3.16 The comparison of non-dimensional heat flux distributions at (a) Re = 197,000 , α = 15 ; (b) Re = 187,000 , α = 0 .............................................53. x.

(14) Nomenclature P P∞ P0 Pb. ： pressure ： pressure of freestream far field stagnation pressure base pressure temperature temperature of freestream far field velocity velocity of freestream far field density. ρ. ：：：：：：：. ρ∞ μ. ： density of freestream far field ： viscosity. α φ. ：：：：：：：：：：：：：. T T∞ V V∞. Re. M Cn Ca Cm C nα X CP L ΔX q q0 S l. angle of attack angle of roll Reynolds number Mach number normal-force coefficient axial-foecr coefficient pitching-moment coefficient normal-force curve slope center of pressure location model length character length of HB-2 model heat flux stagnation heat flux. ： area of cross section ： diameter of cylinder. xi.

(15) Chapter 1. Introduction. 1.1 Motivation and Background In recent years the rocket’s development has become one of the focal point which a lot of countries pay close attention to. With different payload, the rocket to be allowed to carry out the different task like the satellite launch, air sounding, space probe, various types of experiment and so on. The rocket design is one kind of the conformity technologies that contains much knowledge. It has contained the aerodynamics, structure analysis, control system, propulsion system and so on. In the design process, one but had decided the mission and the rocket flight path, then have often decided the majority of designs like the rocket outlook size, propelling power and so on. We may see the Fig.1.1, the rocket flies through various stages such as subsonic, sonic, supersonic, and hypersonic. Atmospheric density reduced rapidly with increasing of height, the rocket flies over the continuous flow, transitional flow, and rare flow at the same time too. The rocket under the different flow, the different speed and the different shape also can have the different aerodynamics forces influence. These aerodynamics forces influences all needs to go overcomes when designing the rocket. Because of nowadays computer computational ability is great strength, these data which aerodynamics forces influences to the rocket already did not need the affiliation to obtain again by experiments. By computer simulation, we may obtain the data that error in the scope which may accept. Compare. 1.

(16) with doing the experiment to survey, computer simulation may save hundred time of above the expenditure. This is a very big progress in the rocket design. In the simulation, the accuracy and the times which simulates is very important. In the same time of simulation, the different method of simulation can obtain the data with different accuracy. Therefore in order to obtain a more accurate data, we need to begin to improve from the method of simulation and production of grid. In the same method of simulation, the different time of simulation also can obtain the data with different accuracy. Basically the more time spent, the more accurate data that we can obtain. But we need to weigh the increase of accuracy of data and time used. The first step is simulating successfully to aerodynamics designing of rocket. Then the next step is improving accuracy and time of simulation.. 1.2 Literatures Survey Because the computational ability of computer has progressed, aerodynamics simulation is applied in many science and technology like car, airplane, rocket and so on. In 2001 Fumiya Togashi [1] uses overset unstructured grids to simulate supersonic airplane/booster separation. An unstructured grid around the rocket booster is overset on the stationary grid around the airplane and moves with time to simulate the separation process. Some results are shown in Fig.1.2, Fig.1.3, and Fig.1.4. In 1963 J. Don Gray [2] do the force tests of standard hypervelocity ballistic models. 2.

(17) HB-1 and HB-2, that is shown in Fig.1.5. They used the two models to test the accuracy of some wind tunnels. In 1999 J. Reuther [3] do the application of a control theory-based aerodynamic shape optimization method do the problem of supersonic aircraft design. A high fidelity computational fluid dynamics (CFD) algorithm modeling the Euler equations is used to calculate the aerodynamic properties of complex three-dimensional aircraft configurations. Some results are shown in Fig.1.6 and Fig.1.7.. 1.3 Specific Objectives of the Thesis Based on previous reviews, the current objectives of the thesis are summarized as follows: 1. We do the grid convergence test of HB-2 model with some different quality of grids to choose suitable mesh file with UNIC-UNS code [4]. 2. We do the heat transfer simulation and aerodynamics simulation of HB-2 model with different Mach number and different attack angles with the suitable mesh file. 3. We compare the results with experimental data and discuss it to verify the ability of simulation of UNIC-UNS code at supersonic flow.. 3.

(18) Chapter 2. Numerical Method. In this thesis, we use the UNIC-UNS code, developed by Y.S. Chen et al, to simulate an unsteady compressible flow. It uses Navier-Stokes solver with finite volume method. The governing equation, boundary condition, numerical methods, algorithm and so on will be discussed below.. 2.1 Governing Equations The general form of mass conservation, energy conservation, Navier-Stokes equation and other transport equations can be written in Cartesian tensor form:. ⎛ ∂ (ρφ ) ∂ (ρU jφ ) = ∂ ⎜⎜ μφ ∂φ + ∂x j ⎝ ∂x j ∂x j ∂t. ⎞ ⎟ + Sφ ⎟ ⎠. (1). where μφ is an effective diffusion coefficient, S φ denotes the source term, ρ is the fluid density and φ = (1, u, v, w, h, k, ε ) stands for the variables for the mass, momentum, total energy and turbulence equation, respectively.. 2.2 Spatial Discretization The cell-centered scheme is employed here then the control volume surface can be represented by the cell surfaces and the coding structure can be much simplified. The transport equations can also be written in integral form as:. r r ∂ ρφ d Ω + F ∫Γ ⋅ ndΓ = Ω∫ SΩ dΩ ∂t Ω∫. (2) 4.

(19) r where Ω is the domain of interest, Γ the surrounding surface, n the unit normal in. r outward direction. The flux function F consists of the inviscid and the viscous parts: r v F = ρVφ − μφ ∇φ. (3). The finite volume formulation of flux integral can be evaluated by the summation of the flux vectors over each face, r r F ∫ ⋅ ndΓ = Γ. ∑F. j = k (i ). i, j. ΔΓ j. (4). where k(i) is a list of faces of cell i, Fi,j represents convection and diffusion fluxes through the interface between cell i and j, ΔΓ j is the cell-face area. The viscous flux for the face e between control volumes P and E as shown in Fig.2.1 can be approximated as:. r r ⎛ r rE − rP φE − φP r ⎜ (∇φ ⋅ n )e ≈ r r + ∇φe ⋅ ⎜ n − r r rE − rP rE − r ⎝. ⎞ ⎟ ⎟ ⎠. (5). That is based on the consideration that. r. r. φ E − φ P ≈ ∇φe ⋅ (rE − rP ). (6). where ∇φ is interpolated from the neighbor cells E and P. The inviscid flux is evaluated through the values at the upwind cell and a linear reconstruction procedure to achieve second order accuracy. r. r. φe = φu + Ψe ∇φu ⋅ (re − ru ). (7). where the subscript u represents the upwind cell and Ψe is a flux limiter used to prevent from local extrema introduced by the data reconstruction. The flux limiter proposed by 5.

(20) Barth [5] is employed in this work. Defining φmax = max (φu , φ j ), φmin = min (φu , φ j ) , the. scalar Ψe associated with the gradient at cell u due to edge e is ⎧ ⎛ φ max − φu ⎞ 0 ⎟ifφ e − φ > 0 ⎪min⎜⎜1, 0 ⎟ ⎝ φ e − φu ⎠ ⎪ ⎪ ⎛ φ − φu ⎞ 0 Ψe = ⎨ min⎜1, min ⎟ ⎜ φ 0 − φ ⎟ifφ e − φ < 0 ⎪ e u ⎠ ⎝ ⎪ ⎪ ⎩1. (8). where φ e0 is computed without the limiting condition (i.e. Ψe =1). 2.3 Time Integration A general implicit discretized time-marching scheme for the transport equations can be written as: ⎛ ρn ⎞ n+1 NB ( ρφ P )n n +1 ⎜⎜ + AP ⎟⎟φ P = ∑ Amφm + + Sφ Δt m =1 ⎝ Δt ⎠. (9). where NB means the neighbor cells of cell P. The high order differencing terms and cross diffusion terms are treated using known quantities and retained in the source term and updated explicitly. The Δ-form used for time-marching in this work can be written as: NB ⎛ ρn ⎞ ⎜⎜ + AP ⎟⎟Δφ P = ∑ Am Δφm + SU φ m =1 ⎝ Δt ⎠. (10). NB ⎛ ⎞ ⎜ Sφ + ∑ Am Δφmn − APφ n ⎟ m =1 ⎠ SU φ = ⎝. (11). θ. where θ is a time-marching control parameter which needs to specify. θ = 1 and θ = 0.5 are 6.

(21) for implicit first-order Euler time-marching and second-order time-centered time-marching schemes. The above derivation is good for non-reacting flows. For general applications, a dual-time sub-iteration method is now used in UNIC-UNS for time-accurate time-marching computations.. 2.4 Pressure-Velocity-Density Coupling In an extended SIMPLE [6] family pressure-correction algorithm, the pressure correction equation for all-speed flow is formulated using the perturbed equation of state, momentum and continuity equations. The simplified formulation can be written as:. ρ′ =. ρ′ γRT. r r r r ; u ′ = − Du ∇p′; u n+1 = u n + u ′; p n+1 = p n + p′. (12). n. r r rn ∂ρ ′ ⎛ ∂ρ ⎞ + ∇(u ρ ′) + ∇(ρu ′) = −⎜ ⎟ − ∇(ρu ) ∂t ⎝ ∂t ⎠. (13). where Du is the pressure-velocity coupling coefficient. Substituting Eq. (12) into Eq. (13), the following all-speed pressure-correction equation is obtained, n. rn 1 p′ ⎛ Δρ ′ ⎞ ⋅ + ∇ ⋅ (ρDu ∇p ′) = −⎜ ⎟ − ∇ ⋅ (ρu ) γRT Δt ⎝ Δt ⎠. (14). For the cell-centered scheme, the flux integration is conducted along each face and its contribution is sent to the two cells on either side of the interface. Once the integration loop is performed along the face index, the discretization of the governing equations is completed. First, the momentum equation (9) is solved implicitly at the predictor step.. 7.

(22) Once the solution of pressure-correction equation (14) is obtained, the velocity, pressure and density fields are updated using Eq. (12). The entire corrector step is repeated 2 and 3 times so that the mass conservation is enforced. The scalar equations such as turbulence transport equations, species equations etc. are then solved sequentially. Then, the solution procedure marches to the next time level for transient calculations or global iteration for steady-state calculations. Unlike for incompressible flow, the pressure-correction equation, which contains both convective and diffusive terms is essentially transport-like. All treatments for inviscid and the viscous fluxes described above are applied to the corresponding parts in Eq. (14).. 2.5 Linear Matrix Solver The discretized finite-volume equations can be represented by a set of linear algebra equations, which are non-symmetric matrix system with arbitrary sparsity patterns. Due to the diagonal dominant for the matrixes of the transport equations, they can converge even through the classical iterative methods. However, the coefficient matrix for the pressure-correction equation may be ill conditioned and the classical iterative methods may break down or converge slowly. Because satisfaction of the continuity equation is of crucial importance to guarantee the overall convergence, most of the computing time in fluid flow calculation is spent on solving the pressure-correction equation by which the. 8.

(23) continuity-satisfying flow field is enforced. Therefore the preconditioned Bi-CGSTAB [7] and GMRES [8] matrix solvers are used to efficiently solve, respectively, transports equation and pressure-correction equation.. 2.6 Parallelization Compared with a structured grid approach, the unstructured grid algorithm is more memory and CPU intensive because “links” between nodes, faces, cells, needs to be established explicitly, and many efficient solution methods developed for structured grids such as approximate factorization, line relaxation, SIS, etc. cannot be used for unstructured methods. As a result, numerical simulation of three-dimensional flow fields remains very expensive even with today’s high-speed computers. As it is becoming more and more difficult to increase the speed and storage of conventional supercomputers, a parallel architecture wherein many processors are put together to work on the same problem seems to be the only alternative. In theory, the power of parallel computing is unlimited. It is reasonable to claim that parallel computing can provide the ultimate throughput for large-scale scientific and engineering applications. It has been demonstrated that performance that rivals or even surpasses supercomputers can be achieved on parallel computers.. 9.

(24) Chapter 3 Results and Discussions 3.1 Overview In the thesis, we want to simulate the flow outside HB-2 model with UNIC-UNS. It is very important to choose suitable quality of grids. Therefore first we must do grid convergence test to choose suitable quality of grids. Then we use the grids to do the aerodynamics simulation of HB-2model at different velocity and angle of attack. Last we compare the results with experimental data and discuss it.. 3.2 Grid Convergence Test The quality of grids will affect the time and the accuracy of simulation. Finer grids will obtain more accurate results but cost more time of simulation. Therefore we simulate the flow at the same flow conditions with different quality of grids. We compare the results with experimental data to choose the most suitable grids. Then we do the following research with the mesh file.. 3.2.1 Grid Configuration The grid configuration is shown in Fig.3.1. The finer grid is developed near the boundary of wall of HB-2 model in order to result more complex behaviors. Relatively, the coarser grid is used at the boundary of freestream far-field to reduced computational cost. The construction of grid needs to be balanced between computational cost and accurate of 10.

(25) solutions, and most of the grid configuration depends on trial and error. In chapter 3.1 of the thesis, a series of grid testing demonstrates to show how to optimize the mesh. The thesis focus on the aerodynamics physic on the rocket, hence, the feedback of baseflow can be reasonable neglected. In that sense, the baseflow zone is not included in the grid because of its heavy load of calculation and less contribution of the aerodynamics properties.. 3.2.2 Flow Conditions and Simulation Conditions Table.I lists the flow conditions of all the simulations of Mach number, Reynolds number, angle of attack, temperature, pressure, density, and viscosity. The case with M = 3.01 we choose is in order to observe the shock wave. The physical phenomenon is. just happen in supersonic flow. The form of shock wave is related to Mach number. The case with α = 8 0 we choose is in order to observe the location that the shock wave happen in. If we choose the case at α = 0 0 , the shock wave is happen in the nose of HB-2 model as we know. If we choose the case at α ≠ 0 0 , the position that the shock wave is happen in will change slightly. The other flow conditions pressure, density, and viscosity are in the cause of getting the same Reynolds number as experimental data. So we can do the following comparison. Table.II gives the simulation conditions for grid testing, within the cell numbers and computation time. The minimum grid size is different in the three. 11.

(26) cases. They are 200 times ΔX , 50 times ΔX , and 30 times ΔX separately. There is wall function in UNIC-UNS. It makes the result accurate even if the minimum grid size of the case is dozens of times of ΔX . The residual is shown in Fig.3.2. Grid testing criteria is the residual being reduced at least one order and variation less than 1% . It shows that the result has already convergent to that we request.. 3.2.3 Validation 3.2.3.1 Density, Pressure and Mach Number Distributions The density, pressure and Mach number distribution at various slices of the computational domain are shown in Fig.3.3. We can observe that there is shock wave in the nose of the HB-2 model. As the flow pass the shock wave region, the pressure increase, the density increase and the velocity decrease suddenly. For the attack angle of these case is not zero, the stagnation point leave the top of nose of HB-2 model slightly. After the flow leave the shack wave region, the temperature, density, and velocity of flow recover original value gradually. Because the attack angle is not zero in the case, there are higher temperature and higher density and lower velocity in the windward side of wall then the leeward side of wall. The temperature, pressure, and density of the flow near the corner of wall have larger change.. 12.

(27) 3.2.3.2 Coefficients of Normal Force The schematic drawing of aerodynamics coefficients is shown in Fig.3.4. The Normal-force coefficient is one of the aerodynamics coefficients. The comparison of normal-force coefficient between the results and experimental data is shown in Fig.3.5. Cn is defined as,. Cn =. Fn 1 ρ ∞V∞ 2 S 2. , where Fn is normal-force, ρ ∞ is density of freestream far-field, V∞ is velocity of freestream far-field, S is the area of cross section. By the comparison, we can observe that the error in case1 is less than 2%. The errors in case2 and case3 are all less than 0.5%, which is better than the error in case1. The remainder of the normal-force coefficient between case1 and case2 is about 30 times of the remainder of the normal-force coefficient between case2 and case3.. 3.2.3.3 Coefficients of Pitching Moment The pitching-moment coefficient is also one of the aerodynamics coefficients. The pitching-moment coefficient’s comparison between the results and experimental data is shown in Fig.3.6. Cm is defined as,. Cm =. Mp 1 ρ ∞V∞ 2 Sl 2. 13.

(28) , where Mp is pitching moment, ρ ∞ is density of freestream far-field, V∞ is velocity of freestream far-field, S is the area of cross section, l is the diameter of cylinder. We take the moment reference as original point to calculate Cm . By the comparison, we can observe that the error in case1 is less than 5%. And the errors in case2 and case3 are all less than 2%, which is better than case1. The remainder of the normal-force coefficient between case1 and case2 is about 35 times of the remainder of the normal-force coefficient between case2 and case3.. 3.3 Aerodynamics Simulation with Different Mach numbers and Angles of Attack 3.3.1 Results in Different Mach Numbers and Attack Angles By the result of the grid convergence test, we choose the mesh file with 300,000 grids to simulate the cases with two kinds of Mach numbers M = 3.01 , M = 5.10 . And each kind of case contains five kinds of attack angles α = 0 0 , α = 2 0 , α = 5 0 , α = 10 0 , and. α = 15 0 . We also compare the results with experimental data and discuss the results with different Mach numbers and attack angles. We also explain the physical meaning of each coefficient.. 3.3.1.1 Flow Conditions and Simulation Conditions The flow conditions are shown in Table.III. The flows with two kinds of Mach numbers 14.

(29) M = 3.01 , M = 5.10 are all supersonic flow. So we can compare the shock wave between. different Mach number. We also can compare the shock wave between different angles of attack. The simulation conditions are shown in Table.IV. Because we have already done the grid convergence test, the simulation conditions are the same as the case of grid convergence test. The residual of case’s criteria is also equal.. 3.3.1.2 Density, Pressure and Mach Number Distributions The density, pressure and Mach number distribution are shown in Fig.3.7. We compare the cases that involve two kinds of attack angles α = 0 0 and α = 15 0 . Their remainder of attack angle is the biggest. It makes us more convenient to observe the difference between them. We first observe the location of shock wave. The location of shock wave of the case at α = 0 0 is just in the top of nose, and the location of shock wave of the case at α = 15 0 leave the top of nose slightly to the windward side of wall. In the case at α = 0 0 , all changes of flow near the wall are axial symmetry include temperature, density, and velocity. In the case at α = 15 0 , there are different change between the windward side of wall and the leeward side of wall--there are higher temperature, higher density, and lower velocity in the windward side of wall then the leeward side of wall for example.. 15.

(30) 3.3.1.3 Coefficients of Axial-Force The results of axial-force coefficient are shown in Fig.3.8. Ca is defined as, Ca =. Fa 1 ρ ∞V∞ 2 S 2. , where Fa is axial-force, ρ ∞ is density of freestream far-field, V∞ is velocity of freestream far-field, S is the area of cross section, l is the diameter of cylinder. When the flow conditions are fixed, Ca is affected by Fa . In the same geometry form and flow conditions, higher Ca means that the axial-force of rocket is higher and the rock need higher thrust to fly. We can observe that the change of axial-force coefficient is very small between different attack angles. It is because the cases with less attack angle have less influence on axial-force. There is almost no change of axial-force coefficient. The results of simulation doesn’t include base flow region. In order to do some comparison, we set Pb = P∞ to correct the axial-force. The error of the correctional axial-force coefficients is less than 10%. In fact, Pb is lower than P∞ . The correct change line of C a should be between the line at Pb = 0 and the line at Pb = P∞ . We do the comparison between two kinds cases of Mach numbers. Although the case at M = 5.10 has higher velocity than the case at M = 3.01 , the case at M = 5.10 has lower density than the case at M = 3.01 . The pressure of wall in the cases at M = 5.10 is always higher than the pressure of wall in the cases at M = 3.01 . It makes the higher normal-force in the cases at M = 3.01 , and also makes the higher normal-force coefficient in the cases at M = 3.01 . 16.

(31) 3.3.1.4 Coefficients of Normal-Force The results of axial-force coefficient are shown in Fig.3.9. Cn is affected by Fn . The HB-2 model is axial symmetry. If the attack angle of HB-2 model is zero, the normal-force on wall of HB-2 model is also zero. When the angle of attack is larger, the normal-force on wall of HB-2 model is also larger if other conditions like pressure, density, velocity of flow are the same. We can observe that when the attack angle is bigger, the coefficient is also bigger. The error in the case at α = 2 0 is less then 9%. The errors in other cases are all less than 5%. We do the comparison between two kinds cases of Mach numbers. Although the case at M = 5.10 has higher velocity than the case at M = 3.01 , the case at M = 5.10 has lower density than the case at M = 3.01 . The pressure of wall in the cases at M = 5.10 is always higher than the pressure of wall in the cases at M = 3.01 . It makes the higher normal-force in the cases at M = 3.01 , and also makes the. higher normal-force coefficient in the cases at M = 3.01 .. 3.3.1.5 Coefficients of Pitching-Moment The results of pitching-moment coefficient are shown in Fig.3.10. The minus sign means the counterclockwise direction. Cm is affected by Mp and Mp is affected by Fa and Fn . The influence of Fn on Cm is larger than Fa . When the angle of attack. is larger, the total force that on wall of HB-2 model is also larger. It makes 17.

(32) Pitching-Moment that on the wall of HB-2 model larger. It is also that we can observe by the figure. The error in the case of M = 5.10 and α = 2 0 is about 40%. The error in the case of M = 5.10 and α = 5 0 is about 20%. The other errors in these cases are less than 8%. Perhaps the reason that makes the larger error is separation happened at the leeward side of wall.. 3.3.1.6 Normal-Force Curve Slope The results of the normal-force curve slope are shown in Fig.3.11. C nα is defined as the slope of normal-force coefficient at α = 0 0 . When the rocket flies actually, the normal-force curve slope affects the sensitivity of control. Precise control makes the trajectory that we want to achieve the mission of flying. C nα in the case at lower Mach number is larger. When the Mach number is higher and the Reynolds number is almost no change, the density is lower. The influence of density on normal-force is larger than velocity on normal-force. So C nα in the case at higher Mach number is lower than C nα in the case at lower Mach number. The errors in these cases are less than 5%.. 3.3.1.7 Center of Pressure Location The results of center of pressure location are shown in Fig.3.12. X CP is defined as, X CP = Cm Cn. 18.

(33) We revise X CP become the value that we get by taking the top of nose of HB-2 model as the original point. The center of pressure location determines the stability of the rocket’s flying directly. If the center of pressure location is more close to nose of rocket than the center of gravity, the rocket is unstable while it is flying. If the center of pressure location is farther from nose of rocket than the center of gravity, the rocket is stable while it is flying. We can observe that when the attack angle is bigger and the center of pressure location is also bigger. The center of pressure of the case at M = 3.01 is large than it of the case at M = 5.10 . It is because the two kinds cases of Mach number have almost the same Cm , but Cn in the case at M = 3.01 is higher than Cn in the case at M = 5.10 . It mean that the HB-2 model at M = 5.10 is more stable than the HB-2 model at M = 3.01 .. 3.3.1.8 Comparison of Non-Dimensional Pressure Distributions The results of pressure distribution are shown in Fig.3.13. The pressure is higher in the nearby nose, it is probably P P0 = 0.8 ~ 0.9 . After X / L pass 0.1, the pressure will lower to P P0 = 0.1 . It is because the flow region near the wall just breaks away the shock region and the velocity is higher to supersonic once again. The pressure rises a little bit after X / L pass 0.6. It is because the diameter of cylinder of HB-2 model is becomes large and the velocity that perpendicular to the wall is also become large. So the pressure 19.

(34) rises a little bit after X / L pass 0.6. The pressure in the case at φ = 0 0 are higher than the pressure in the case at φ = 90 0 and φ = 180 0 . It is because there is cross flow effect. The 3-D pressure distribution are shown in Fig.3.14.We can observe the cross flow effect. The most errors in these cases are about 5%. A part of errors in the case at X / L = 0.4 ~ 1 and φ = 180 0 are higher than 30%. It is because there is separation point, and the separation point affects the flow, therefore the deviation is higher.. 3.4 Heat Transfer Simulation There is aerodynamic influence when the rocket flies, besides there is influence of heat. When the rocket flies with the high speed like to hypersonic, there is high temperature of air that near the body of rocket. The temperature is so high that can melt the wall of rocket. If we can get the heat flux data of rocket while it is flying with simulation, we can strengthen the structure of body to prevent the rocket being destroyed by high temperature.. 3.4.1 Flow Conditions and Simulation Conditions The flow conditions are shown in Table.V. The flow conditions of the experimental data of the paper [9] are not complete. It doesn’t include wall’s temperature. The wall temperature is very important to heat transfer simulation. It wall make different temperature difference between wall and flow. The result of heat flux of simulation is also different. So we set some wall’s temperature to simulate and compare it. The wall’s 20.

(35) temperature is fixed that is different from the aerodynamics cases. This is for approaching experiment situation to get more accurate results. We also simulate with laminar flow and turbulence in order to observe the difference between them. The simulation conditions are shown in Table.VI.. 3.4.2 Temperature Distributions The temperature distribution is shown in Fig.3.15. There highest temperature in the stagnation point nearby the nose of HB-2 model. The temperature boundary layer in the windward side of wall is less than it in the leeward side of wall. It wall make the heat flux in the windward side of wall higher than the heat flux in the leeward side of wall. By the comparison between two kinds of wall temperature, there is higher difference between the wall temperature and flow temperature near the wall in the case at Tw = 53.6k . It wall make the higher heat flux in the wall in the case at Tw = 53.6k .. 3.4.3 Comparison of Non-Dimensional Heat Flux Distributions The results of heat transfer distribution are shown in Fig.3.16. The heat flux is higher in the nearby nose, it is probably q q0 = 0.8 ~ 0.9 . After X / L pass 0.05, the heat flux will lower to q q 0 = 0.08 . The heat flux rises a little bit after X / L pass 0.6. By the comparison we can observe that the result of case at Tw = 353.6k and laminar flow is. 21.

(36) more accurate then other cases. It is because most flow region is laminar flow in the case at Re = 197,000 . The turbulence will raise the heat flux, so that the heat flux at turbulence is higher than the heat flux at laminar flow. The wall temperature of case4 is relatively close to experimental situation then other cases. Just the change of wall temperature wall makes such a big change of heat flux. We can know again the wall temperature is very important to heat transfer simulation. If we want to get more accurate results, we can set about from the wall temperature and grid file.. 22.

(37) Chapter 4. Conclusions. 4.1 Summary The current study can be briefly summarized as follow: 1. The influence of change of attack angle on Cm and Cn are greater than the influence of change of attack angle on Ca . 2. The higher velocity makes the higher stability of HB-2 model at the similar Re . 3. The error in most aerodynamics coefficients that we get is less than 4% and it prove that UNIC-UNS code has very good accuracy in aerodynamics simulation at supersonic continuity flow. 4. Although there are no very accurate results in heat transfer simulations, the tendency of change and physical phenomena are all rational. It is a main reason that the flow conditions of experimental data are not intact enough. 5. The results and mesh can be reference material to supply the people of follow-up study. Using finer grids to focus on the flow region with complicated change is also a good topic that is worth probing into for the people of follow-up study.. 23.

(38) References [1] Fumiya Togashi, “Flow simulation of NAL experimental supersonic airplane/booster separation using overset unstructured grids,” Computers & Fluids 30 (2001) 673-688 [2] J. Don Gray, “Force Tests of Standard Hypervelocity Ballistic Models HB-1 and HB-2 at Mach 1.5 to 10,” AEDC-TDR-63-137 [3] J. Reuther, “Aerodynamic Shape Optimization of Supersonic Aircraft Configurations via an Adjoint Formulation on Distributed Memory Parallel Computers,” Computer & Fluids 28 (1999) 675-700 [4] Chen, Y.S.,”UNIC-UNS USER’S MANUAL” [5] Barth, T.J., “Recent Development in High Order K-Exact Reconstruction on Unstructured Meshes,” AIAA Paper 93-0668, 1993 [6] Chen, Y.S., “An Unstructured Finite Volume Method for Viscous Flow Computations,” 7th International Conference on Finite Element Methods in Flow Problems, Feb. 3-7, 1989, University of Alabama in Huntsville, Huntsville, Alabama [7] Van Der Vorst, H.A., “Bi-CFSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the solution of Nonsymmetric Linear Systems,” SIAM J. Sci. Stat. Comput., Vol. 13(2), pp. 631-644, 1992 [8] Saad, Y. and Schultz, M.H., “GMRES: A Generalized Minimal Residual Algorithm for solving Nonsymmeric Linear Systems,” SIAM J. Sci. Stat. Comput., Vol. 7(3), pp. 856-869, 1986 [9] Shigeru KUCHI-ISHI, “Compartive Force/Heat Flux Measurements between JAXA Hypersonic Test Facilities Using Standard Model HB-2,” ISSN 1349-1113 JAXA-RR-04-035E. 24.

(39) Tables Table.I Flow conditions for grid convergence test Case. Cells. 1. 100,000. 2. 300,000. 3. 500,000. M. α(deg). 5.1. 8. Re. T(k) P(atm) ρ(kg/m^3) μ(kg/m*s). 2.20E+06 198.9 0.055. 0.0889. Wall Temperature. 2.98E-05. Adiabatic. Table.II Simulation conditions for grid convergence test Case. Cells. cpus Computation Time(hrs). 1. 100,000. 2. 300,000. 3. 500,000. Simulation Time(s). 3.3 12. 7.5. Number of Time Step (Time Step Size) 1~2000(E-8), 2000~12000(2E-7). 2.02E-03. 13.3. 1~2000(E-8), 2000~12000(2E-7) 1~2000(E-8), 2000~12000(2E-7). Table.III Flow conditions for simulation with different Mach numbers and angles of attack M α(deg) Re T(k) P(atm) ρ(kg/m^3) μ(kg/m*s) 3.01. 5.10. 0 2 5 10 15. 2.50E+06. 199.3. 0.05. 0.0819. 1.43E-05. 2.20E+06. 198.9. 0.055. 0.0889. 2.98E-05. 25.

(40) Table.IV Simulation conditions for simulation with different Mach numbers and angles of attack Cells cpus Computation Simulation Number of Time Step Time(hrs) Time(s) (Time Step Size) 300,000. 12. 7.5. 2.02E-03. 1~2000(E-8), 2000~12000(2E-7). Table.V Flow conditions for heat transfer simulation ρ α Case M 1. (deg). T∞ (k) P(atm). (kg/m^3). μ (kg/m*s) TW (k) 53.6. 2 3. Re. 15. 1.97E+05. 53.6. 7.39E-04. 4.87E-03. 1.70E-05. 4. Turbulence. 353.6 53.6. 9.59. Flow model. Laminar flow. 353.6. 5. 0. 1.87E+05. 55.1. 7.36E-04. 4.72E-03. 1.76E-05. 6. 55.1. Laminar flow. 300. Table.VI Simulation conditions for heat transfer simulation Case Cells cpus Computation Time Simulation Time Number of Time Step (hrs) (s) (Time Step Size) 1~6. 300,000. 12. 7.5. 2.02E-03. 26. 1~2000(E-8), 2000~12000(2E-7).

(41) Figures. Fig.1.1 The flying diagram of Formosat-2. Fig.1.2 Enlarged overset grid around the nodes of the objects cut view at X = 0.22 and. symmetric boundary.. 27.

(42) Fig.1.3 Intergrid boundaries of the booster subgrid: angles of attack of the airplane at 2o and the booster at 0o and ΔX = 0 : ΔZ = 0.4 , ΔZ = 1.0 , and ΔZ = 3.0 .. Fig.1.4 Comparison between computed Mach contours and schlieren photographs of supersonic airplane/rocket booster separation at M ∞ = 2.5 , angles of attack of. the airplane at 2o and the booster at 0o, and ΔX = 0 : ΔZ = 0.4 , ΔZ = 2.4 , and ΔZ = 5.0 m. 28.

(43) Fig.1.5 Standard hypervelocity ballistic correlation configurations.. Fig.1.6 Generic supersonic transport configuration SYN87-MB Grid Structure: 180 Blocks. 29.

(44) Fig.1.7 Generic supersonic transport configuration SYN87-MB solution pressure coefficient (Mach=2.2 α = 3.15 0 C L = 0.105 ). Fig.2.1 Unstructured control volume.. 30.

(45) 31.

(46) Fig.3.1 Typical grid distribution of aerodynamics simulation of HB-2 model.. Fig.3.2 Typical time history of residuals of mass, momentum and energy equations.. 32.

(47) (a). 33.

(48) (b). 34.

(49) (c) Fig.3.3 The distributions of (a) density; (b) pressure; (c) Mach number at various slices of the computational domain.. 35.

(50) Fig.3.4 The diagram of normal-force and axial-force.. Fig.3.5 The comparison of normal-force coefficients of grid convergence test.. 36.

(51) Fig.3.6 The comparison of pitching-moment coefficients of grid convergence test.. 37.

(52) (a). 38.

(53) (b). 39.

(54) (c). 40.

(55) (d). 41.

(56) (e). 42.

(57) (f) Fig.3.7 The distributions of (a) density at α = 0 0 ; (b) density at α = 15 0 ; (c) pressure at α = 0 0 ; (d) pressure at α = 15 0 ; (e) Mach number at α = 0 0 ; (f) Mach number at α = 15 0 at various slices of the computational domain. 43.

(58) (a). (b) Fig.3.8 The comparison of axial-force coefficients at (a) M = 3.01 ; (b) M = 5.10. 44.

(59) (a). (b) Fig.3.9 The comparison of normal-force coefficients at (a) M = 3.01 ; (b) M = 5.10. 45.

(60) (a). (b) Fig.3.10 The comparison of pitching-moment coefficients at (a) M = 3.01 ; (b) M = 5.10. 46.

(61) Fig.3.11 The comparison of normal-force curve slope at different Mach number and Reynolds number.. Fig.3.12 The center of pressure location at different angles of attack.. 47.

(62) (a). (b). 48.

(63) (c). (d) Fig.3.13 The comparison of non-dimensional pressure distributions at (a) M = 3.01 ,. α = 5 0 ; (b) M = 3.01 , α = 15 0 ; (c) M = 5.10 , α = 5 0 ; (d) M = 5.10 , α = 15 0 49.

(64) (a). (b) Fig.3.14 The 3-D pressure distribution at (a) M = 3.01 , α = 15 0 ; (b) M = 5.10 , α = 15 0. 50.

(65) (a). 51.

(66) (b) Fig.3.15 The distributions of temperature at (a) Tw = 53.6k ; (b) Tw = 353.6k. 52.

(67) (a). (b) Fig.3.16 The comparison of non-dimensional heat flux distributions at (a) Re = 197,000 , α = 15 ; (b) Re = 187,000 , α = 0 53.

(68)