Submit manuscript...
eISSN: 2576-4500

Aeronautics and Aerospace Open Access Journal

Research Article Volume 1 Issue 3

Simulation of the supersonic starting jet flow interactions with a rigid body

Pooya Azizian, Milad Azarmanesh

Department of Mechanical Engineering, Babol Noshirvani University of Technology, Iran

Correspondence: Pooya Azizian, Department of Mechanical Engineering, Babol Noshirvani University of Technology, Shariati St, Babol, Iran, Tel +989124578951

Received: September 01, 2017 | Published: November 3, 2017

Citation: Azizian P, Azarmanesh M. Simulation of the supersonic starting jet flow interactions with a rigid body. Aeron Aero Open Access J. 2017;1(3):129-134. DOI: 10.15406/aaoaj.2017.01.00017

Download PDF

Abstract

The interaction of starting flow fields of an axisymmetric supersonic under-expanded jet with the cross section of a cylinder is simulated. To achieve this goal, an implicit pressure correction algorithm based on a finite volume method on a collocated grid for simulation of compressible flows is presented. A two dimensional test problem, the supersonic flow through a wind tunnel with a step, is simulated and the numerical outcomes agree very well with previous results. Then, the method is extended to axisymmetric flows and validated for a supersonic circular jet flow. Finally, the proposed method is used to analyze the concerned problem. To describe the evolution of the starting flow fields, the Mach number, pressure, density and temperature contours evolutions through time are investigated. Also, the evolution of total force and temperature distribution on the cylinder cross section surface are studied during the collision of mushroom head of supersonic jet with the surface of cylinder.

Keywords: finite volume method, starting flow, supersonic flow, axisymmetric flow, gauss theorem, mach disk, prandtl-meyer expansion

Introduction

Understanding of the jet aerodynamics is vital to several areas of implementation of the fluid dynamic science. Due to the significant applications, supersonic under-expanded compressible jets have been investigated widely.1-3 Also, because of many important physical phenomena, the starting flow fields of supersonic jets have been studied previously.4,5 Moreover, interaction of the starting flow fields of supersonic jets with rigid bodies are found in many areas of aerospace and mechanical engineering; however there are not many works reported in this field.

This paper focuses on unsteady supersonic flows and is concerned with the development of a simple implicit pressure correction method for computing compressible fluids and extending it to the ability of handle supersonic axisymmetric flows, which contain strong shocks. We aim to simulate and analyze the interaction of starting flow fields of an axisymmetric supersonic under-expanded jet with the cross section of a cylinder.

Supersonic jet flow characteristics can be very complex. Due to the supersonic nozzle exit Mach number, jet exit pressure can be different from ambient pressure. This pressure difference between the jet and the ambient fluid must be resolved locally either across an oblique shock, by a prominent streamline curvature at the jet boundary, or by a Mach disk inside the jet.6 In the starting flow of a supersonic jet, the shock diffraction leads to the misalignment of the pressure and density gradients and makes a vortex ring which moves downstream like a mushroom head.5,7 The interactions of the acoustic wave and the mushroom head with the rigid body lead to propagate the acoustic waves from the impingement region out and toward the nozzle exit. Also, overexpansion of the flow occurs at the cylinder lip.

This paper is organized as follows. In Section 2, the numerical method, which consists of the governing equations and the algorithm, is described. In Section 3, the method is validated for the supersonic flow through a wind tunnel with a step and the under-expanded circular jet. Also, the simulation results of the interaction of starting flow fields of an axisymmetric supersonic under-expanded jet with the cross section of a cylinder are provided. Finally, a brief summary and conclusions are given in the last Section.

Numerical method

Governing equations

In high speed flow, the Reynolds number is extremely high and turbulence effect is confined to thin boundary layer that leads to small frictional drag compared to pressure drag and wave drag. If frictional drag is ignored, high speed flow can be computed using the inviscid Euler equations.8 Continuity, momentum and energy equations for an ideal gas and inviscid flow can be written as follow:

ρ d t d + ( ρ d u j d ) x j d =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaqcLbsapeGaeyOaIyRaeqyWdiNcpaWaaWbaaSqabKqa GeaajugWa8qacaWGKbaaaaGcpaqaaKqzGeWdbiabgkGi2kaadshak8 aadaahaaWcbeqcbasaaKqzadWdbiaadsgaaaaaaKqzGeGaey4kaSIc daWcaaWdaeaajugib8qacqGHciITkmaabmaapaqaaKqzGeWdbiabeg 8aYPWdamaaCaaaleqajeaibaqcLbmapeGaamizaaaajugibiaadwha l8aadaqhaaqcbasaaKqzadWdbiaadQgaaKqaG8aabaqcLbmapeGaam izaaaaaOGaayjkaiaawMcaaaWdaeaajugib8qacqGHciITcaWG4bWc paWaa0baaKqaGeaajugWa8qacaWGQbaajeaipaqaaKqzadWdbiaads gaaaaaaKqzGeGaeyypa0JaaGimaaaa@5DA0@                               (1)

( ρ d u i d ) t d + ( ρ d u i d u j d ) x j d = P d x i d MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaqcLbsapeGaeyOaIyRcdaqadaWdaeaajugib8qacqaH bpGCk8aadaahaaWcbeqcbasaaKqzadWdbiaadsgaaaqcLbsacaWG1b WcpaWaa0baaKqaGeaajugWa8qacaWGPbaajeaipaqaaKqzadWdbiaa dsgaaaaakiaawIcacaGLPaaaa8aabaqcLbsapeGaeyOaIyRaamiDaO WdamaaCaaaleqajeaibaqcLbmapeGaamizaaaaaaqcLbsacqGHRaWk kmaalaaapaqaaKqzGeWdbiabgkGi2QWaaeWaa8aabaqcLbsapeGaeq yWdiNcpaWaaWbaaSqabKqaGeaajugWa8qacaWGKbaaaKqzGeGaamyD aSWdamaaDaaajeaibaqcLbmapeGaamyAaaqcbaYdaeaajugWa8qaca WGKbaaaKqzGeGaamyDaSWdamaaDaaajeaibaqcLbmapeGaamOAaaqc baYdaeaajugWa8qacaWGKbaaaaGccaGLOaGaayzkaaaapaqaaKqzGe WdbiabgkGi2kaadIhal8aadaqhaaqcbasaaKqzadWdbiaadQgaaKqa G8aabaqcLbmapeGaamizaaaaaaqcLbsacqGH9aqpcqGHsislkmaala aapaqaaKqzGeWdbiabgkGi2kaadcfak8aadaahaaWcbeqcbasaaKqz adWdbiaadsgaaaaak8aabaqcLbsapeGaeyOaIyRaamiEaSWdamaaDa aajeaibaqcLbmapeGaamyAaaqcbaYdaeaajugWa8qacaWGKbaaaaaa aaa@7AB1@          (2)

( ρ d E d ) t d + ( ( ρ d E d + P d ) u j d ) x j d =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaqcLbsapeGaeyOaIyRaaiikaiabeg8aYPWdamaaCaaa leqajeaibaqcLbmapeGaamizaaaajugibiaadweak8aadaahaaWcbe qcbasaaKqzadWdbiaadsgaaaqcLbsacaGGPaaak8aabaqcLbsapeGa eyOaIyRaamiDaOWdamaaCaaaleqajeaibaqcLbmapeGaamizaaaaaa qcLbsacqGHRaWkkmaalaaapaqaaKqzGeWdbiabgkGi2QWaaeWaa8aa baWdbmaabmaapaqaaKqzGeWdbiabeg8aYPWdamaaCaaaleqajeaiba qcLbmapeGaamizaaaajugibiaadweak8aadaahaaWcbeqcbasaaKqz adWdbiaadsgaaaqcLbsacqGHRaWkcaWGqbGcpaWaaWbaaSqabKqaGe aajugWa8qacaWGKbaaaaGccaGLOaGaayzkaaqcLbsacaWG1bWcpaWa a0baaKqaGeaajugWa8qacaWGQbaajeaipaqaaKqzadWdbiaadsgaaa aakiaawIcacaGLPaaaa8aabaqcLbsapeGaeyOaIyRaamiEaSWdamaa DaaajeaibaqcLbmapeGaamOAaaqcbaYdaeaajugWa8qacaWGKbaaaa aajugibiabg2da9iaaicdaaaa@6DF7@       (3)

where the superscript ‘d’ represents dimensional values. Also, the variables ρ, u i MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaKqzGeaeaaaaaa aaa8qacqaHbpGCcaGGSaGaaGPaVlaadwhak8aadaWgaaqcbasaaKqz adWdbiaadMgaaSWdaeqaaaaa@3E95@  and E show the density, velocity component in the direction of i and the total energy per unit mass, which can be defined as below:

E=e+ u i u i 2 = c v T+ u i u i 2 = RT γ1 + u i u i 2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaKqzGeaeaaaaaa aaa8qacaWGfbGaeyypa0JaamyzaiabgUcaROWaaSaaa8aabaqcLbsa peGaamyDaOWdamaaBaaajeaibaqcLbmapeGaamyAaaWcpaqabaqcLb sapeGaamyDaOWdamaaBaaajeaibaqcLbmapeGaamyAaaWcpaqabaaa keaajugib8qacaaIYaaaaiabg2da9iaadogak8aadaWgaaqcbasaaK qzadWdbiaadAhaaSWdaeqaaKqzGeWdbiaadsfacqGHRaWkkmaalaaa paqaaKqzGeWdbiaadwhak8aadaWgaaqcbasaaKqzadWdbiaadMgaaS WdaeqaaKqzGeWdbiaadwhak8aadaWgaaqcbasaaKqzadWdbiaadMga aSWdaeqaaaGcbaqcLbsapeGaaGOmaaaacqGH9aqpkmaalaaapaqaaK qzGeWdbiaadkfacaWGubaak8aabaqcLbsapeGaeq4SdCMaeyOeI0Ia aGymaaaacqGHRaWkkmaalaaapaqaaKqzGeWdbiaadwhak8aadaWgaa qcbasaaKqzadWdbiaadMgaaSWdaeqaaKqzGeWdbiaadwhak8aadaWg aaqcbasaaKqzadWdbiaadMgaaSWdaeqaaaGcbaqcLbsapeGaaGOmaa aaaaa@6834@                 (4)

that by using the equation of state ( P d = ρ d R d T d ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaKqzGeaeaaaaaa aaa8qacaGGOaGaamiuaOWdamaaCaaaleqajeaibaqcLbmapeGaamiz aaaajugibiabg2da9iabeg8aYPWdamaaCaaaleqajeaibaqcLbmape Gaamizaaaajugibiaadkfak8aadaahaaWcbeqcbasaaKqzadWdbiaa dsgaaaqcLbsacaWGubGcpaWaaWbaaSqabKqaGeaajugWa8qacaWGKb aaaOWdaiaacMcaaaa@49BC@  the above equation leads to:

ρE= P γ1 + 1 2 ρ u i u i MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaKqzGeaeaaaaaa aaa8qacqaHbpGCcaWGfbGaeyypa0JcdaWcaaWdaeaajugib8qacaWG qbaak8aabaqcLbsapeGaeq4SdCMaeyOeI0IaaGymaaaacqGHRaWkkm aalaaapaqaaKqzGeWdbiaaigdaaOWdaeaajugib8qacaaIYaaaaiab eg8aYjaadwhak8aadaWgaaqcbasaaKqzadWdbiaadMgaaSWdaeqaaK qzGeWdbiaadwhak8aadaWgaaqcbasaaKqzadWdbiaadMgaaSWdaeqa aaaa@4DAA@                                               (5)

where γ MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape Gaeq4SdCgaaa@380D@  is the ratio of specific heats, and by using this equation in energy equation (eq. 3), the number of variables are reduced.

In the aforementioned governing equations, the conservation forms are used that the fluxes of mass, momentum and energy into and out of the control volume, themselves become important dependent variables rather than just the primitive variables. This is due to the fact that there are discontinuities in primitive variables across the shock; however, the fluxes are constant across the shock wave.9

To non-dimensionalization, the reference variables are used as follow:

  u i = u i d u r ,     ρ= ρ d ρ r ,     t= t d L/ u r ,     P= P d ρ r u r 2 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaKqzGeaeaaaaaa aaa8qacaGGGcGaamyDaOWdamaaBaaajeaibaqcLbmapeGaamyAaaWc paqabaqcLbsapeGaeyypa0JcdaWcaaWdaeaajugib8qacaWG1bWcpa Waa0baaKqaGeaajugWa8qacaWGPbaajeaipaqaaKqzadWdbiaadsga aaaak8aabaqcLbsapeGaamyDaOWdamaaBaaajeaibaqcLbmapeGaam OCaaWcpaqabaaaaKqzGeWdbiaacYcacaGGGcGaaiiOaiaacckacaGG GcGaaiiOaiabeg8aYjabg2da9OWaaSaaa8aabaqcLbsapeGaeqyWdi NcpaWaaWbaaSqabKqaGeaajugWa8qacaWGKbaaaaGcpaqaaKqzGeWd biabeg8aYPWdamaaBaaajeaibaqcLbmapeGaamOCaaWcpaqabaaaaK qzGeWdbiaacYcacaGGGcGaaiiOaiaacckacaGGGcGaaiiOaiaadsha cqGH9aqpkmaalaaapaqaaKqzGeWdbiaadshak8aadaahaaWcbeqcba saaKqzadWdbiaadsgaaaaak8aabaqcLbsapeGaamitaiaac+cacaWG 1bGcpaWaaSbaaKqaGeaajugWa8qacaWGYbaal8aabeaaaaqcLbsape GaaiilaiaacckacaGGGcGaaiiOaiaacckacaGGGcGaamiuaiabg2da 9OWaaSaaa8aabaqcLbsapeGaamiuaOWdamaaCaaaleqajeaibaqcLb mapeGaamizaaaaaOWdaeaajugib8qacqaHbpGCk8aadaWgaaqcbasa aKqzadWdbiaadkhaaSWdaeqaaKqzGeWdbiaadwhal8aadaqhaaqcba saaKqzadWdbiaadkhaaKqaG8aabaqcLbmapeGaaGOmaaaaaaaaaa@8875@        (6)

Where u r , ρ r MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaKqzGeaeaaaaaa aaa8qacaWG1bGcpaWaaSbaaKqaGeaajugWa8qacaWGYbaal8aabeaa jugibiaacYcapeGaeqyWdiNcpaWaaSbaaKqaGeaajugWa8qacaWGYb aal8aabeaaaaa@4064@ and L are reference velocity, density and length, respectively. In this way, the non-dimensional governing equations are:

ρ t + ( ρ u j ) x j =0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaqcLbsapeGaeyOaIyRaeqyWdihak8aabaqcLbsapeGa eyOaIyRaamiDaaaacqGHRaWkkmaalaaapaqaaKqzGeWdbiabgkGi2Q WaaeWaa8aabaqcLbsapeGaeqyWdiNaamyDaOWdamaaBaaajeaibaqc LbmapeGaamOAaaWcpaqabaaak8qacaGLOaGaayzkaaaapaqaaKqzGe WdbiabgkGi2kaadIhak8aadaWgaaqcbasaaKqzadWdbiaadQgaaSWd aeqaaaaajugib8qacqGH9aqpcaaIWaaaaa@504C@                                                                    (7)

( ρ u i ) t + ( ρ u i u j ) x j = P x i MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaqcLbsapeGaeyOaIyRcdaqadaWdaeaajugib8qacqaH bpGCcaWG1bGcpaWaaSbaaKqaGeaajugWa8qacaWGPbaal8aabeaaaO WdbiaawIcacaGLPaaaa8aabaqcLbsapeGaeyOaIyRaamiDaaaacqGH RaWkkmaalaaapaqaaKqzGeWdbiabgkGi2QWaaeWaa8aabaqcLbsape GaeqyWdiNaamyDaOWdamaaBaaajeaibaqcLbmapeGaamyAaaWcpaqa baqcLbsapeGaamyDaOWdamaaBaaajeaibaqcLbmapeGaamOAaaWcpa qabaaak8qacaGLOaGaayzkaaaapaqaaKqzGeWdbiabgkGi2kaadIha k8aadaWgaaqcbasaaKqzadWdbiaadQgaaSWdaeqaaaaajugib8qacq GH9aqpcqGHsislkmaalaaapaqaaKqzGeWdbiabgkGi2kaadcfaaOWd aeaajugib8qacqGHciITcaWG4bGcpaWaaSbaaKqaGeaajugWa8qaca WGPbaal8aabeaaaaaaaa@637F@                                                    (8)

t ( P+ ( γ1 ) 2 ρ u i u i )+ x j ( ( γP+ ( γ1 ) 2 ρ u i u i ) u j )=0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaqcLbsapeGaeyOaIylak8aabaqcLbsapeGaeyOaIyRa amiDaaaakmaabmaapaqaaKqzGeWdbiaadcfacqGHRaWkkmaalaaapa qaa8qadaqadaWdaeaajugib8qacqaHZoWzcqGHsislcaaIXaaakiaa wIcacaGLPaaaa8aabaqcLbsapeGaaGOmaaaacqaHbpGCcaWG1bGcpa WaaSbaaKqaGeaajugWa8qacaWGPbaal8aabeaajugib8qacaWG1bGc paWaaSbaaKqaGeaajugWa8qacaWGPbaal8aabeaaaOWdbiaawIcaca GLPaaajugibiabgUcaROWaaSaaa8aabaqcLbsapeGaeyOaIylak8aa baqcLbsapeGaeyOaIyRaamiEaOWdamaaBaaajeaibaqcLbmapeGaam OAaaWcpaqabaaaaOWdbmaabmaapaqaa8qadaqadaWdaeaajugib8qa cqaHZoWzcaWGqbGaey4kaSIcdaWcaaWdaeaapeWaaeWaa8aabaqcLb sapeGaeq4SdCMaeyOeI0IaaGymaaGccaGLOaGaayzkaaaapaqaaKqz GeWdbiaaikdaaaGaeqyWdiNaamyDaOWdamaaBaaajeaibaqcLbmape GaamyAaaWcpaqabaqcLbsapeGaamyDaOWdamaaBaaajeaibaqcLbma peGaamyAaaWcpaqabaaak8qacaGLOaGaayzkaaqcLbsacaWG1bGcpa WaaSbaaKqaGeaajugWa8qacaWGQbaal8aabeaaaOWdbiaawIcacaGL Paaajugibiabg2da9iaaicdaaaa@78DA@    (9)

Algorithm

 The variables are collocated in space at the center of control volumes. Moreover, the numerical method is implicit in time. The governing equations can be discretized by integrating over a control volume and using the Gauss theorem as follows:

ρ c.v. t+1 ρ c.v. t Δt + 1 V faces ρ face t+1 u N t+1 A face =0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaqcLbsapeGaeqyWdi3cpaWaa0baaKqaGeaajugWa8qa caWGJbGaaiOlaiaadAhacaGGUaaajeaipaqaaKqzadWdbiaadshacq GHRaWkcaaIXaaaaKqzGeGaeyOeI0IaeqyWdi3cpaWaa0baaKqaGeaa jugWa8qacaWGJbGaaiOlaiaadAhacaGGUaaajeaipaqaaKqzadWdbi aadshaaaaak8aabaqcLbsapeGaeyiLdqKaamiDaaaacqGHRaWkkmaa laaapaqaaKqzGeWdbiaaigdaaOWdaeaajugib8qacaWGwbaaaOWaay buaeqajeaipaqaaKqzadWdbiaadAgacaWGHbGaam4yaiaadwgacaWG Zbaaleqan8aabaqcLbsapeGaeyyeIuoaaiabeg8aYTWdamacObyhaa qcbasaiGgGjugWa8qacGaoKoOzaiac4q6GHbGaiGdPdogacGaoKoyz aaqcbaYdaeacObycLbmapeGaiGgGdshacWaAaA4kaSIaiGgGigdaaa qcLbsacaWG1bWcpaWaa0baaKqaGeaajugWa8qacaWGobaajeaipaqa aKqzadWdbiaadshacqGHRaWkcaaIXaaaaKqzGeGaamyqaOWdamaaBa aajeaibaqcLbmapeGaamOzaiaadggacaWGJbGaamyzaaWcpaqabaqc LbsapeGaeyypa0JaaGimaaaa@81DC@                     (10)

g i,c.v. t+1 g i,c.v. t Δt + 1 V faces g i,face t+1 u N t+1 A face = P c.v. t+1 x i MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaqcLbsapeGaam4zaSWdamaaDaaajeaibaqcLbmapeGa amyAaiaacYcacaWGJbGaaiOlaiaadAhacaGGUaaajeaipaqaaKqzad WdbiaadshacqGHRaWkcaaIXaaaaKqzGeGaeyOeI0Iaam4zaSWdamaa DaaajeaibaqcLbmapeGaamyAaiaacYcacaWGJbGaaiOlaiaadAhaca GGUaaajeaipaqaaKqzadWdbiaadshaaaaak8aabaqcLbsapeGaeyiL dqKaamiDaaaacqGHRaWkkmaalaaapaqaaKqzGeWdbiaaigdaaOWdae aajugib8qacaWGwbaaaOWaaybuaeqajeaipaqaaKqzadWdbiaadAga caWGHbGaam4yaiaadwgacaWGZbaaleqan8aabaqcLbsapeGaeyyeIu oaaiaadEgal8aadaqhaaqcbasaaKqzadWdbiaadMgacaGGSaGaamOz aiaadggacaWGJbGaamyzaaqcbaYdaeaajugWa8qacaWG0bGaey4kaS IaaGymaaaajugibiaadwhal8aadaqhaaqcbasaaKqzadWdbiaad6ea aKqaG8aabaqcLbmapeGaamiDaiabgUcaRiaaigdaaaqcLbsacaWGbb GcpaWaaSbaaKqaGeaajugWa8qacaWGMbGaamyyaiaadogacaWGLbaa l8aabeaajugib8qacqGH9aqpcqGHsislkmaalaaapaqaaKqzGeWdbi abgkGi2kaadcfak8aadaqhaaqcbasaaKqzadWdbiaadogacaGGUaGa amODaiaac6caaKqaG8aabaqcLbmapeGaamiDaiabgUcaRiaaigdaaa aak8aabaqcLbsapeGaeyOaIyRaamiEaOWdamaaBaaajeaibaqcLbma peGaamyAaaWcpaqabaaaaaaa@8CFC@      (11)

( P+ ( γ1 ) 2 ρ u i u i ) c.v. t+1 ( P+ ( γ1 ) 2 ρ u i u i ) c.v. t Δt + 1 V faces ( γP+ ( γ1 ) 2 ρ u i u i ) face t+1 u N t+1 A face =0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOabaeqabaaeaaaaaa aaa8qadaWcaaWdaeaapeWaaeWaa8aabaqcLbsapeGaamiuaiabgUca ROWaaSaaa8aabaWdbmaabmaapaqaaKqzGeWdbiabeo7aNjabgkHiTi aaigdaaOGaayjkaiaawMcaaaWdaeaajugib8qacaaIYaaaaiabeg8a Yjaadwhak8aadaWgaaqcbasaaKqzadWdbiaadMgaaSWdaeqaaKqzGe Wdbiaadwhak8aadaWgaaqcbasaaKqzadWdbiaadMgaaSWdaeqaaaGc peGaayjkaiaawMcaa8aadaqhaaqcbasaaKqzadWdbiaadogacaGGUa GaamODaiaac6caaKqaG8aabaqcLbmapeGaamiDaiabgUcaRiaaigda aaqcLbsacqGHsislkmaabmaapaqaaKqzGeWdbiaadcfacqGHRaWkkm aalaaapaqaa8qadaqadaWdaeaajugib8qacqaHZoWzcqGHsislcaaI XaaakiaawIcacaGLPaaaa8aabaqcLbsapeGaaGOmaaaacqaHbpGCca WG1bGcpaWaaSbaaKqaGeaajugWa8qacaWGPbaal8aabeaajugib8qa caWG1bGcpaWaaSbaaKqaGeaajugWa8qacaWGPbaal8aabeaaaOWdbi aawIcacaGLPaaapaWaa0baaKqaGeaajugWa8qacaWGJbGaaiOlaiaa dAhacaGGUaaajeaipaqaaKqzadWdbiaadshaaaaak8aabaqcLbsape GaeyiLdqKaamiDaaaaaOqaaKqzGeGaaCzcaiaaxMaacqGHRaWkkmaa laaapaqaaKqzGeWdbiaaigdaaOWdaeaajugib8qacaWGwbaaaOWaay buaeqajeaipaqaaKqzadWdbiaadAgacaWGHbGaam4yaiaadwgacaWG Zbaaleqan8aabaqcLbsapeGaeyyeIuoaaOWaaeWaa8aabaqcLbsape Gaeq4SdCMaamiuaiabgUcaROWaaSaaa8aabaWdbmaabmaapaqaaKqz GeWdbiabeo7aNjabgkHiTiaaigdaaOGaayjkaiaawMcaaaWdaeaaju gib8qacaaIYaaaaiabeg8aYjaadwhak8aadaWgaaqcbasaaKqzadWd biaadMgaaSWdaeqaaKqzGeWdbiaadwhak8aadaWgaaqcbasaaKqzad WdbiaadMgaaSWdaeqaaaGcpeGaayjkaiaawMcaaSWdamaaDaaajeai baqcLbmapeGaamOzaiaadggacaWGJbGaamyzaaqcbaYdaeaajugWa8 qacaWG0bGaey4kaSIaaGymaaaajugibiaadwhal8aadaqhaaqcbasa aKqzadWdbiaad6eaaKqaG8aabaqcLbmapeGaamiDaiabgUcaRiaaig daaaqcLbsacaWGbbGcpaWaaSbaaKqaGeaajugWa8qacaWGMbGaamyy aiaadogacaWGLbaal8aabeaajugib8qacqGH9aqpcaaIWaaaaaa@B50C@  (12)

where the superscripts denote the time steps. Also, the subscripts c,v and face represent the values at the center of control volumes and the values that located on the control volume surfaces, respectively. Here, V and A face MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaKqzGeaeaaaaaa aaa8qacaWGbbGcpaWaaSbaaKqaGeaajugWa8qacaWGMbGaamyyaiaa dogacaWGLbaal8aabeaaaaa@3D1A@  determine the volume and surface area of the control volumes. Furthermore, g i MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaKqzGeaeaaaaaa aaa8qacaWGNbGcpaWaaSbaaKqaGeaajugWa8qacaWGPbaal8aabeaa aaa@3A8B@  is the momentum in the direction of i and   u N MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaKqzGeaeaaaaaa aaa8qacaGGGcGaamyDaOWdamaaBaaajeaibaqcLbmapeGaamOtaaWc paqabaaaaa@3BA2@  is the face normal velocity, which is acquired by projection and nor that interpolation.10

To find the values on the control volume surfaces, which there are in the convection term of the governing equations, the following scheme is used:

Q face=i+ 1 2 = ( 1+ w i ) Q i +( 1 w i ) Q i+1 2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaKqzGeaeaaaaaa aaa8qacaWGrbGcpaWaaSbaaKqaGeaajugWa8qacaWGMbGaamyyaiaa dogacaWGLbGaeyypa0JaamyAaiabgUcaRSWaaSaaaKqaG8aabaqcLb mapeGaaGymaaqcbaYdaeaajugWa8qacaaIYaaaaaWcpaqabaqcLbsa peGaeyypa0JcdaWcaaWdaeaapeWaaeWaa8aabaqcLbsapeGaaGymai abgUcaRiaadEhak8aadaWgaaqcbasaaKqzadWdbiaadMgaaSWdaeqa aaGcpeGaayjkaiaawMcaaKqzGeGaamyuaOWdamaaBaaajeaibaqcLb mapeGaamyAaaWcpaqabaqcLbsapeGaey4kaSIcdaqadaWdaeaajugi b8qacaaIXaGaeyOeI0Iaam4DaOWdamaaBaaajeaibaqcLbmapeGaam yAaaWcpaqabaaak8qacaGLOaGaayzkaaqcLbsacaWGrbGcpaWaaSba aKqaGeaajugWa8qacaWGPbGaey4kaSIaaGymaaWcpaqabaaakeaaju gib8qacaaIYaaaaaaa@623D@                     (13)

Where Q can be each of the terms on the control volume surfaces. At the vicinity of discontinuity wi takes the value of +1 and -1 depends on the flow direction, but in the rest of the domain it takes the value of zero.

Adapted from the researches of Hou & Mahesh,10,11 the algorithm is based on a pressure-correction method, which the face normal velocities are projected to satisfy a constraint on the divergence that is determined by the energy equation. This method, from the outermost layer to inside, has a time loop to advance from the time step t to t+1 an outer loop to advance from the iteration k to k+1 and three inner loops to solve the governing equations by an iterative approach. The iterative governing equations are as below:

ρ c.v. t+1,k+1 ρ c.v. t Δt + 1 V faces ρ face t+1,k+1 u N t+1,k A face =0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaqcLbsapeGaeqyWdi3cpaWaa0baaKqaGeaajugWa8qa caWGJbGaaiOlaiaadAhacaGGUaaajeaipaqaaKqzadWdbiaadshacq GHRaWkcaaIXaGaaiilaiaadUgacqGHRaWkcaaIXaaaaKqzGeGaeyOe I0IaeqyWdi3cpaWaa0baaKqaGeaajugWa8qacaWGJbGaaiOlaiaadA hacaGGUaaajeaipaqaaKqzadWdbiaadshaaaaak8aabaqcLbsapeGa eyiLdqKaamiDaaaacqGHRaWkkmaalaaapaqaaKqzGeWdbiaaigdaaO Wdaeaajugib8qacaWGwbaaaOWaaybuaeqajeaipaqaaKqzadWdbiaa dAgacaWGHbGaam4yaiaadwgacaWGZbaaleqan8aabaqcLbsapeGaey yeIuoaaiabeg8aYTWdamacObyhaaqcbasaiGgGjugWa8qacGaoKoOz aiac4q6GHbGaiGdPdogacGaoKoyzaaqcbaYdaeacObycLbmapeGaiG gGdshacWaAaA4kaSIaiGgGigdacGaAakilaiacOb4GRbGamGgGgUca RiacObiIXaaaaKqzGeGaamyDaSWdamaaDaaajeaibaqcLbmapeGaam OtaaqcbaYdaeaajugWa8qacaWG0bGaey4kaSIaaGymaiaacYcacaWG RbaaaKqzGeGaamyqaOWdamaaBaaajeaibaqcLbmapeGaamOzaiaadg gacaWGJbGaamyzaaWcpaqabaqcLbsapeGaeyypa0JaaGimaaaa@8DB6@                                          (14)

g i,c.v. t+1,* g i,c.v. t Δt + 1 V faces g i,face t+1,* u N t+1,k A face = P c.v. t+1,k x i MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaqcLbsapeGaam4zaSWdamaaDaaajeaibaqcLbmapeGa amyAaiaacYcacaWGJbGaaiOlaiaadAhacaGGUaaajeaipaqaaKqzad WdbiaadshacqGHRaWkcaaIXaGaaiilaiaacQcaaaqcLbsacqGHsisl caWGNbWcpaWaa0baaKqaGeaajugWa8qacaWGPbGaaiilaiaadogaca GGUaGaamODaiaac6caaKqaG8aabaqcLbmapeGaamiDaaaaaOWdaeaa jugib8qacqGHuoarcaWG0baaaiabgUcaROWaaSaaa8aabaqcLbsape GaaGymaaGcpaqaaKqzGeWdbiaadAfaaaGcdaGfqbqabKqaG8aabaqc LbmapeGaamOzaiaadggacaWGJbGaamyzaiaadohaaSqab0Wdaeaaju gib8qacqGHris5aaGaam4zaSWdamaaDaaajeaibaqcLbmapeGaamyA aiaacYcacaWGMbGaamyyaiaadogacaWGLbaajeaipaqaaKqzadWdbi aadshacqGHRaWkcaaIXaGaaiilaiaacQcaaaqcLbsacaWG1bWcpaWa a0baaKqaGeaajugWa8qacaWGobaajeaipaqaaKqzadWdbiaadshacq GHRaWkcaaIXaGaaiilaiaadUgaaaqcLbsacaWGbbGcpaWaaSbaaKqa GeaajugWa8qacaWGMbGaamyyaiaadogacaWGLbaal8aabeaajugib8 qacqGH9aqpcqGHsislkmaalaaapaqaaKqzGeWdbiabgkGi2kaadcfa k8aadaqhaaqcbasaaKqzadWdbiaadogacaGGUaGaamODaiaac6caaK qaG8aabaqcLbmapeGaamiDaiabgUcaRiaaigdacaGGSaGaam4Aaaaa aOWdaeaajugib8qacqGHciITcaWG4bGcpaWaaSbaaKqaGeaajugWa8 qacaWGPbaal8aabeaaaaaaaa@92F8@                                  (15)

( P+ ( γ1 ) 2 ρ u i u i ) c.v. t+1,k+1 ( P+ ( γ1 ) 2 ρ u i u i ) c.v. t Δt + 1 V faces ( γP+ ( γ1 ) 2 ρ u i u i ) face t+1,k+1 u N t+1,k+1 A face =0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOabceqabaaC=haaqa aaaaaaaaWdbmaalaaapaqaa8qadaqadaWdaeaajugib8qacaWGqbGa ey4kaSIcdaWcaaWdaeaapeWaaeWaa8aabaqcLbsapeGaeq4SdCMaey OeI0IaaGymaaGccaGLOaGaayzkaaaapaqaaKqzGeWdbiaaikdaaaGa eqyWdiNaamyDaOWdamaaBaaajeaibaqcLbmapeGaamyAaaWcpaqaba qcLbsapeGaamyDaOWdamaaBaaajeaibaqcLbmapeGaamyAaaWcpaqa baaak8qacaGLOaGaayzkaaWdamaaDaaajeaibaqcLbmapeGaam4yai aac6cacaWG2bGaaiOlaaqcbaYdaeaajugWa8qacaWG0bGaey4kaSIa aGymaiaacYcacaWGRbGaey4kaSIaaGymaaaajugibiabgkHiTOWaae Waa8aabaqcLbsapeGaamiuaiabgUcaROWaaSaaa8aabaWdbmaabmaa paqaaKqzGeWdbiabeo7aNjabgkHiTiaaigdaaOGaayjkaiaawMcaaa Wdaeaajugib8qacaaIYaaaaiabeg8aYjaadwhak8aadaWgaaqcbasa aKqzadWdbiaadMgaaSWdaeqaaKqzGeWdbiaadwhak8aadaWgaaqcba saaKqzadWdbiaadMgaaSWdaeqaaaGcpeGaayjkaiaawMcaa8aadaqh aaqcbasaaKqzadWdbiaadogacaGGUaGaamODaiaac6caaKqaG8aaba qcLbmapeGaamiDaaaaaOWdaeaajugib8qacqGHuoarcaWG0baaaaGc baqcLbsacaWLjaGaaCzcaiabgUcaROWaaSaaa8aabaqcLbsapeGaaG ymaaGcpaqaaKqzGeWdbiaadAfaaaGcdaGfqbqabKqaG8aabaqcLbma peGaamOzaiaadggacaWGJbGaamyzaiaadohaaSqab0Wdaeaajugib8 qacqGHris5aaGcdaqadaWdaeaajugib8qacqaHZoWzcaWGqbGaey4k aSIcdaWcaaWdaeaapeWaaeWaa8aabaqcLbsapeGaeq4SdCMaeyOeI0 IaaGymaaGccaGLOaGaayzkaaaapaqaaKqzGeWdbiaaikdaaaGaeqyW diNaamyDaOWdamaaBaaajeaibaqcLbmapeGaamyAaaWcpaqabaqcLb sapeGaamyDaOWdamaaBaaajeaibaqcLbmapeGaamyAaaWcpaqabaaa k8qacaGLOaGaayzkaaWcpaWaa0baaKqaGeaajugWa8qacaWGMbGaam yyaiaadogacaWGLbaajeaipaqaaKqzadWdbiaadshacqGHRaWkcaaI XaGaaiilaiaadUgacqGHRaWkcaaIXaaaaKqzGeGaamyDaSWdamaaDa aajeaibaqcLbmapeGaamOtaaqcbaYdaeaajugWa8qacaWG0bGaey4k aSIaaGymaiaacYcacaWGRbGaey4kaSIaaGymaaaajugibiaadgeak8 aadaWgaaqcbasaaKqzadWdbiaadAgacaWGHbGaam4yaiaadwgaaSWd aeqaaKqzGeWdbiabg2da9iaaicdaaaaa@C035@  (16)

where the superscript * shows the predicted values at the iteration k+1. By using equation (13) and extracting ρ c.v. t+1,k+1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaKqzGeaeaaaaaa aaa8qacqaHbpGCl8aadaqhaaqcbasaaKqzadWdbiaadogacaGGUaGa amODaiaac6caaKqaG8aabaqcLbmapeGaamiDaiabgUcaRiaaigdaca GGSaGaam4AaiabgUcaRiaaigdaaaaaaa@44EB@  from equation (14) and g i,c.v. t+1,* MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaKqzGeaeaaaaaa aaa8qacaWGNbWcpaWaa0baaKqaGeaajugWa8qacaWGPbGaaiilaiaa dogacaGGUaGaamODaiaac6caaKqaG8aabaqcLbmapeGaamiDaiabgU caRiaaigdacaGGSaGaaiOkaaaaaaa@43D6@  from equation (15), the iterative equations for density and momentum predictor will be produced. The following relations can be found by subtracting equation (15) from (11):

g i,c.v. t+1,k+1 = g i,c.v. t+1,* Δt ( ( δP ) x i ) c.v. MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaKqzGeaeaaaaaa aaa8qacaWGNbWcpaWaa0baaeaajugWa8qacaWGPbGaaiilaiaadoga caGGUaGaamODaiaac6caaSWdaeaajugWa8qacaWG0bGaey4kaSIaaG ymaiaacYcacaWGRbGaey4kaSIaaGymaaaajugibiabg2da9iaadEga l8aadaqhaaqaaKqzadWdbiaadMgacaGGSaGaam4yaiaac6cacaWG2b GaaiOlaaWcpaqaaKqzadWdbiaadshacqGHRaWkcaaIXaGaaiilaiaa cQcaaaqcLbsacqGHsislcqGHuoarcaWG0bGcdaqadaWdaeaapeWaaS aaa8aabaqcLbsapeGaeyOaIyRcdaqadaWdaeaajugib8qacqaH0oaz caWGqbaakiaawIcacaGLPaaaa8aabaqcLbsapeGaeyOaIyRaamiEaO WdamaaBaaaleaajugib8qacaWGPbaal8aabeaaaaaak8qacaGLOaGa ayzkaaWdamaaBaaaleaajugWa8qacaWGJbGaaiOlaiaadAhacaGGUa aal8aabeaaaaa@6A00@                             (17)

u i,c.v. t+1,k+1 = u i,c.v. t+1,* Δt ρ c.v. t+1 ( ( δP ) x i ) c.v. MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaceaaW9peaaaaaa aaa8qacaWG1bWdamaaDaaaleaapeGaamyAaiaacYcacaWGJbGaaiOl aiaadAhacaGGUaaapaqaa8qacaWG0bGaey4kaSIaaGymaiaacYcaca WGRbGaey4kaSIaaGymaaaakiabg2da9iaadwhapaWaa0baaSqaa8qa caWGPbGaaiilaiaadogacaGGUaGaamODaiaac6caa8aabaWdbiaads hacqGHRaWkcaaIXaGaaiilaiaacQcaaaGccqGHsisldaWcaaWdaeaa peGaeyiLdqKaamiDaaWdaeaapeGaeqyWdi3damaaDaaaleaapeGaam 4yaiaac6cacaWG2bGaaiOlaaWdaeaapeGaamiDaiabgUcaRiaaigda aaaaaOWaaeWaa8aabaWdbmaalaaapaqaa8qacqGHciITdaqadaWdae aapeGaeqiTdqMaamiuaaGaayjkaiaawMcaaaWdaeaapeGaeyOaIyRa amiEa8aadaWgaaWcbaWdbiaadMgaa8aabeaaaaaak8qacaGLOaGaay zkaaWdamaaBaaaleaapeGaam4yaiaac6cacaWG2bGaaiOlaaWdaeqa aaaa@69DF@                          (18)

u N t+1,k+1 = u N t+1,* Δt ρ face t+1 ( ( δP ) N ) face MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaceaaW9peaaaaaa aaa8qacaWG1bWdamaaDaaaleaapeGaamOtaaWdaeaapeGaamiDaiab gUcaRiaaigdacaGGSaGaam4AaiabgUcaRiaaigdaaaGccqGH9aqpca WG1bWdamaaDaaaleaapeGaamOtaaWdaeaapeGaamiDaiabgUcaRiaa igdacaGGSaGaaiOkaaaakiabgkHiTmaalaaapaqaa8qacqGHuoarca WG0baapaqaa8qacqaHbpGCpaWaa0baaSqaa8qacaWGMbGaamyyaiaa dogacaWGLbaapaqaa8qacaWG0bGaey4kaSIaaGymaaaaaaGcdaqada WdaeaapeWaaSaaa8aabaWdbiabgkGi2oaabmaapaqaa8qacqaH0oaz caWGqbaacaGLOaGaayzkaaaapaqaa8qacqGHciITcaWGobaaaaGaay jkaiaawMcaa8aadaWgaaWcbaWdbiaadAgacaWGHbGaam4yaiaadwga a8aabeaaaaa@60E7@                      (19)

Where N denotes the face normal direction. To find the pressure at the iteration k+1 the following equation must be used:

P c.v. t+1,k+1 = P c.v. t+1,k +δ P c.v. MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaceaaW9peaaaaaa aaa8qacaWGqbWdamaaDaaaleaapeGaam4yaiaac6cacaWG2bGaaiOl aaWdaeaapeGaamiDaiabgUcaRiaaigdacaGGSaGaam4AaiabgUcaRi aaigdaaaGccqGH9aqpcaWGqbWdamaaDaaaleaapeGaam4yaiaac6ca caWG2bGaaiOlaaWdaeaapeGaamiDaiabgUcaRiaaigdacaGGSaGaam 4AaaaakiabgUcaRiabes7aKjaadcfapaWaaSbaaSqaa8qacaWGJbGa aiOlaiaadAhacaGGUaaapaqabaaaaa@5307@ (20)

By substituting equations (17)-(20) into equation (16) and ignoring the higher-order terms of δP MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape GaeqiTdqMaamiuaaaa@38E1@ , the pressure-correction equation can be produced as follow:

a P δ P c.v. + nb a nb δ P nb =RHS MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaceaaW9peaaaaaa aaa8qacaWGHbWdamaaBaaaleaapeGaamiuaaWdaeqaaOWdbiaabs7a caWGqbWdamaaBaaaleaapeGaam4yaiaac6cacaWG2bGaaiOlaaWdae qaaOWdbiabgUcaRmaawafabeWcpaqaa8qacaWGUbGaamOyaaqab0Wd aeaapeGaeyyeIuoaaOGaamyya8aadaWgaaWcbaWdbiaad6gacaWGIb aapaqabaGcpeGaeqiTdqMaamiua8aadaWgaaWcbaWdbiaad6gacaWG IbaapaqabaGcpeGaeyypa0JaamOuaiaadIeacaWGtbaaaa@507F@                                    (21)

Where the subscript nb  represents the neighboring control volumes. The flowchart of the implemented numerical method is shown in Figure 1.

Figure 1 Flowchart of the numerical method.

 Applying boundary conditions and geometry

Ghost cells are used to impose boundary conditions. These cells, which extend the computational domain, are a few additional cells on either end that their values are set at the beginning of each iteration based on the boundary conditions and the interior solution.12

For creating solid geometries inside the domain, a criterion is used to discriminate between the fluid cells and the solid cells. During marching of the computational domain, calculation of the governing equations for the solid cells must be ignored. Also, to calculate the adjacent cells of the solid cells, a sub-program imposes the required values of the solid cells, as boundary ghost cells, based on the fluid cells. Figure 2 shows the flowchart of creating solid geometries.

Figure 2 Flowchart of creating solid geometries.

Results and discussion

Code validation

For validation, a two dimensional test problem, the supersonic flow through a wind tunnel with a step, is simulated. This problem was introduced in the paper by Emery13 and completely investigated by Woodward & Colella14 later. The tunnel heights and lengths are 1 length unit and 3 length units, respectively. The initial condition is a uniform flow with Mach 3, density 1.4 kg/ m 3 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaceaaW9peaaaaaa aaa8qacaaIXaGaaiOlaiaaisdacaGGGcGaam4AaiaadEgacaGGVaGa amyBa8aadaahaaWcbeqaa8qacaaIZaaaaaaa@3FB2@ , pressure 1atm and γ=1.4 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaceaaW9peaaaaaa aaa8qacaqGZoGaeyypa0JaaGymaiaac6cacaaI0aaaaa@3C43@ , which leads to the strong shocks. Also, the step height is 0.2 length units and is located 0.6 length units from the tunnel entrance. The additional boundary conditions near the corner of the step are applied as same as reference14 to minimize numerical errors generated at the corner of the step.

Figure 3 shows the comparison of the density contours for different times with reference14 which indicates remarkable agreement. The dimensionless time is obtained using the length unit and the sound speed.This Figure shows the overexpansion of the flow at the corner of the step and the formation of the lambda shaped shock within the time evolution.

Figure 3 The comparison of the density contours with Reference 14, The CFL number = 1.02, Grid number = 80X240, Dimensionless time = a: 0.5, b: 1.0, c: 2.0, d: 4.0.

Axisymmetric flow

The numerical method for an axisymmetric flow is as same as Section 2 with the exception that the governing equations are as bellow that it ignores circulation ( V θ =0) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaceaaW9peaaaaaa aaa8qacaGGOaGaamOvaSWdamaaBaaabaqcLbmapeGaeqiUdehal8aa beaak8qacqGH9aqpcaaIWaGaaiykaaaa@3F30@ :

ρ t + ( ρu z ) z + 1 r ( ρru r ) r =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaWdbiabgkGi2kaabg8aa8aabaWdbiabgkGi2kaabsha aaGaey4kaSYaaSaaa8aabaWdbiabgkGi2oaabmaapaqaa8qacaqGbp GaaeyDa8aadaWgaaWcbaWdbiaabQhaa8aabeaaaOWdbiaawIcacaGL Paaaa8aabaWdbiabgkGi2kaabQhaaaGaey4kaSYaaSaaa8aabaWdbi aaigdaa8aabaWdbiaabkhaaaWaaSaaa8aabaWdbiabgkGi2oaabmaa paqaa8qacaqGbpGaaeOCaKqzGeGaaeyDaOWdamaaBaaaleaapeGaae OCaaWdaeqaaaGcpeGaayjkaiaawMcaaaWdaeaapeGaeyOaIyRaaeOC aaaacqGH9aqpcaaIWaaaaa@559D@ (22)

( ρu z ) t + ( ρu z u z ) z + 1 r ( ρru r u z ) r = P z MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaWdbiabgkGi2oaabmaapaqaa8qacaqGbpGaaeyDa8aa daWgaaWcbaWdbiaabQhaa8aabeaaaOWdbiaawIcacaGLPaaaa8aaba WdbiabgkGi2kaabshaaaGaey4kaSYaaSaaa8aabaWdbiabgkGi2oaa bmaapaqaa8qacaqGbpGaaeyDa8aadaWgaaWcbaqcLbmapeGaaeOEaa WcpaqabaGcpeGaaeyDa8aadaWgaaWcbaWdbiaabQhaa8aabeaaaOWd biaawIcacaGLPaaaa8aabaWdbiabgkGi2kaabQhaaaGaey4kaSYaaS aaa8aabaWdbiaaigdaa8aabaWdbiaabkhaaaWaaSaaa8aabaWdbiab gkGi2oaabmaapaqaa8qacaqGbpGaaeOCaiaabwhal8aadaWgaaqaaK qzadWdbiaabkhaaSWdaeqaaOWdbiaabwhapaWaaSbaaSqaa8qacaqG 6baapaqabaaak8qacaGLOaGaayzkaaaapaqaa8qacqGHciITcaqGYb aaaiabg2da9iabgkHiTmaalaaapaqaa8qacqGHciITcaqGqbaapaqa a8qacqGHciITcaqG6baaaaaa@6576@   (23)

( ρu r ) t + ( ρu z u r ) z + 1 r ( ρru r u r ) r = P r MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaWdbiabgkGi2oaabmaapaqaa8qacaqGbpGaaeyDa8aa daWgaaWcbaWdbiaabkhaa8aabeaaaOWdbiaawIcacaGLPaaaa8aaba WdbiabgkGi2kaabshaaaGaey4kaSYaaSaaa8aabaWdbiabgkGi2oaa bmaapaqaa8qacaqGbpGaaeyDa8aadaWgaaWcbaWdbiaabQhaa8aabe aak8qacaqG1bWdamaaBaaaleaapeGaaeOCaaWdaeqaaaGcpeGaayjk aiaawMcaaaWdaeaapeGaeyOaIyRaaeOEaaaacqGHRaWkdaWcaaWdae aapeGaaGymaaWdaeaapeGaaeOCaaaadaWcaaWdaeaapeGaeyOaIy7a aeWaa8aabaWdbiaabg8acaqGYbqcLbsacaqG1bGcpaWaaSbaaSqaa8 qacaqGYbaapaqabaGcpeGaaeyDa8aadaWgaaWcbaWdbiaabkhaa8aa beaaaOWdbiaawIcacaGLPaaaa8aabaWdbiabgkGi2kaabkhaaaGaey ypa0ZaaSaaaeaacqGHciITcaWGqbaabaGaeyOaIyRaamOCaaaaaaa@6255@ (24)

( ρE ) t + ( (ρE+P)u Z ) z + ( r(ρE+P)u r ) r =0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape WaaSaaa8aabaWdbiabgkGi2oaabmaapaqaa8qacaqGbpGaaeyraaGa ayjkaiaawMcaaaWdaeaapeGaeyOaIyRaaeiDaaaacqGHRaWkdaWcaa WdaeaapeGaeyOaIy7aaeWaa8aabaWdbiaabIcacaqGbpGaaeyraiaa bUcacaqGqbGaaeykaiaabwhapaWaaSbaaSqaa8qacaqGAbaapaqaba aak8qacaGLOaGaayzkaaaapaqaa8qacqGHciITcaqG6baaaiabgUca Rmaalaaapaqaa8qacqGHciITdaqadaWdaeaapeGaaeOCaiaabIcaca qGbpGaaeyraiaabUcacaqGqbGaaeykaiaabwhapaWaaSbaaSqaa8qa caqGYbaapaqabaaak8qacaGLOaGaayzkaaaapaqaa8qacqGHciITca qGYbaaaiabg2da9iaaicdaaaa@5C95@                  (25)

ρE= P γ1 + 1 2 ρ( u z u z + u r u r ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape GaaeyWdiaabweacqGH9aqpdaWcaaWdaeaapeGaaeiuaaWdaeaapeGa ae4SdiabgkHiTiaaigdaaaGaey4kaSYaaSaaa8aabaWdbiaaigdaa8 aabaWdbiaaikdaaaGaaeyWdmaabmaapaqaa8qacaqG1bWdamaaBaaa leaapeGaaeOEaaWdaeqaaOWdbiaabwhapaWaaSbaaSqaa8qacaqG6b aapaqabaGcpeGaey4kaSIaaeyDa8aadaWgaaWcbaWdbiaabkhaa8aa beaak8qacaqG1bWdamaaBaaaleaapeGaaeOCaaWdaeqaaaGcpeGaay jkaiaawMcaaaaa@4D8A@     (26)

For validation of the written axisymmetric code, the under-expanded circular jet, which containing Mach disk, is simulated. The jet exit Mach number and the exit static pressure ratio ( P exit / P ambient ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqkY=gjVeeu0dXdPqFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape GaaiikaiaadcfapaWaaSbaaSqaa8qacaWGLbGaamiEaiaadMgacaWG 0baapaqabaGcpeGaai4laiaadcfapaWaaSbaaSqaa8qacaWGHbGaam yBaiaadkgacaWGPbGaamyzaiaad6gacaWG0baapaqabaGccaGGPaaa aa@4546@  are 1.5 and 3.15, respectively. The contour in Figure 4 shows the comparison of the obtained Mach number distribution with reference,6 which indicates desirable agreement. A Mach disk is placed at the location of about 4.5 times greater than the jet exit radius from the jet exit location. Also, the radius of the Mach disk is approximately 0.7 times smaller than the jet exit radius. The contour value shows that the Mach number upstream of the Mach disk has accelerated to values about 4.0, but the Mach number downstream of the Mach disk is reduced to values below 0.5.

Figure 4 The comparison of the Mach number contour with Reference 6.

The supersonic starting flow interactions with the cross section of a cylinder

The supersonic starting flow jet interactions with a cross section of a cylinder is simulated in which the jet exit Mach number and the exit static pressure ratio are as same as the last Section equal to 1.5 and 3.15, respectively. The radius of the cylinder is two times greater than the radius of the jet exit plane. The cylinder cross section is located at the 4 length units from the jet exit plane. To describe the evolution of the starting flow fields, the Mach number, pressure, density and temperature evolutions are considered.

Figure 5 shows the Mach number evolution. When the jet ejects, Prandtl-Meyer expansion wave is created at the exit. Because of the expansion, gas accelerates to higher Mach number. At first, the flow structure is as like as a typical supersonic under-expanded planar jet, which moves downstream like a mushroom head. However, this mushroom head collides with the cylinder cross section around8ms. The interactions of the acoustic wave and the mushroom head with the rigid body lead to propagate the acoustic waves from the impingement region out and toward the nozzle exit that disrupts the vortex rings structures. In following, it is observed that the Mach disk is formed at the location of about 3 length units from the jet exit plane. Also, another shock wave is commenced to propagate at the cylinder lip.

Figure 5 The Mach number evolution.

Where the supersonic gas is accelerated, pressure, density and temperature are decreased (Figures 6-8). Interaction of the mushroom head with the cylinder cross section increased pressure distribution significantly at the time of 8ms. It creates a pressure wave in the opposite direction that reduces the pressure distribution on the cylinder surface at following times. Moreover, the temperature contour shows that by crossing the shock width, temperature increases about two times. The highest temperature on the cylinder cross section surface happens about collision of mushroom head.

Figure 6 The pressure contour evolution.

Figure 7 The density contour evolution.

Figure 8 The temperature contour evolution.

Figure 9 indicates the evolution of total force due to jet flow on the cylinder cross section. The highest total force on the cylinder cross section surface occurs just after collision of mushroom head with the surface which is about 2500 kpa. Also, oscillation in the magnitude of total force on the surface is happened in the following times. Figure 10 shows the pressure distribution along the radius of the cylinder cross section surface for different time. It demonstrates reducing in the pressure value in the vicinity of cylinder lip which is evident because of increasing in velocity magnitude adjacent to cylinder lip. In an opposite manner with time of 12ms, depression in the pressure values occurred in time of 8 ms and 16 ms along the surface, which are due to the impacts of vortex rings.

Figure 9 The evolution of total force on the cylinder cross section surface.

Figure 10 The pressure evolution along the radius of the cylinder cross section.

Summary and conclusion

In this study, an implicit pressure correction method for simulation of compressible flows is presented. The presented method is based on a finite volume approach on a collocated grid. The algorithm is explained in the numerical method Section completely and the flowchart of the implemented procedure is provided. The method is validated with remarkable agreement for two-dimensional supersonic flows with strong shock. Moreover, a simple procedure is used to discriminate between the fluid cells and the solid cells, and then it finds the ghost cells location immersed in the computational domain. The no-slip condition along the solid surface is applied on the ghost cells.The method is extended to axisymmetric flows and validated for an axisymmetric under-expanded circular jet. Finally, the interaction of starting flow fields of an axisymmetric supersonic under-expanded jet with the cross section of a cylinder is simulated and analyzed. To describe the evolution of the starting flow fields the Mach number, pressure, density and temperature evolutions are provided. It is observed that the highest force and temperature distribution on the cylinder cross section occur just after collision of mushroom head with the surface.

Acknowledgements

None.

Conflicts of Interest

Author declares that there is no conflict of interest.

References

  1. Bulent K, Otugen MV. Scaling parameters of under-expanded supersonic jets. Phys Fluids. 2002;14(12):4206‒4215.
  2. Vuorinen V, Yu J, Tirunagari S, et al. Large-eddy simulation of highly under-expanded transient gas jets. Phys Fluids. 2013;25(1):016101.
  3. Zhang H, Chen Z, Li B, et al. The secondary vortex rings of a supersonic under-expanded circular jet with low pressure ratio. European Journal of Mechanics-B/Fluids. 2014;46:172‒180.
  4. Golub VV. Development of shock wave and vortex ring structures in unsteady jets. Shock Waves. 1994;3(4):279‒285.
  5. Zhang H, Chen Z, Jiang X, et al. The starting flow structures and evolution of a supersonic planar jet. Computers & Fluids. 2015;114:98‒109.
  6. Pao SP, Abdol-Hamid KS. Numerical Simulation of Jet Aerodynamics Using the Three-Dimensional Navier-Stokes Code PAB3D. NASA Technical Report. 1996. p. 1‒48.
  7. Takayama K, Inoue O. Shock wave diffraction over a 90 degree sharp corner-posters presented at the 18th ISSW. Shock Waves. 1991;1(4):301‒312.
  8. Ferziger JH, Peric M. Compressible Flow. In: Computational methods for fluid dynamics. Chapter 10, Springer, Germany; 2002. p. 309‒328.
  9. Anderson JR, John D. Introduction to Numerical Methods. In: Computational fluid dynamics. Chapter 2, Springer, Germany; 1995. p. 21‒37.
  10. Hou Y. A novel algorithm for DNS/LES of compressible turbulent flows. A dissertation for the degree of doctor of philosophy at university of Minnesota, USA; 2007. p. 88.
  11. Hou Y, Mahesh K. A robust, colocated, implicit algorithm for direct numerical simulation of compressible, turbulent flows. Journal of Computational Physics. 2005;205(1):205‒221.
  12. Leveque RJ. Boundary Conditions and Ghost Cells. In: Finite volume methods for hyperbolic problems. Chapter 7, Cambridge university press, UK; 2002. p. 129‒138.
  13. Emery AE. An evaluation of several differencing methods for inviscid fluid flow problems. Journal of Computational Physics. 1968;2(3):306‒331.
  14. Woodward P, Colella P. The numerical simulation of two-dimensional fluid flow with strong shocks. Journal of Computational Physics. 1984;54(1):115‒173.
Creative Commons Attribution License

©2017 Azizian, et al. This is an open access article distributed under the terms of the, which permits unrestricted use, distribution, and build upon your work non-commercially.