Visualizing 2D Prandtl-Meyer Shock Flow Emanating Over A Sharp Convex Corner using Converge CFD & Post-processing It In ParaView.

Abstract:- The work is focused on visualizing supersonic expansion fan technically known as Prandtl-Meyer shock waves, by allowing high-speed fluid to enter the geometry through an inlet duct followed by a sudden sharp corner, this sharp corner would create an infinite number of expansion waves or Mach lines which would fan out creating a certain flow pattern. We are interested in capturing these flow patterns by using the Adaptive Mesh Refinement (AMR) technique for 3 distinct Sub Grid-Scale (SGS). Pictorial representation of the geometry is shown below in Fig 1.1. The geometry is designed in CAD software & imported to Converge as an STL file which is then edited using the various geometric editing tools available in Converge Studio. Newtonian fluid Air is considered as the working fluent. We will also be making use of boundary flagging techniques for assigning each side of the geometry to a particular boundary & these boundaries will later be assigned to a volumetric region. Next, we set up the simulation parameters using a steady-state solver. Once the setup is complete, we export these inputs files generated by Converge CFD & run it using CYGWIN to generate 3D post output files. These 3D output files are first converted into either ParaView vtk inline binary format or EnSite format for ParaView to read & generate the flow field simulation. The paper will also cover explanations on Mach numbers, basics of shockwaves, different types of shock waves, boundary conditions & it's applications. The simulation will be run for 3 different subgrid-scale :- 0.05, 0.04 & 0.002. We will also be comparing the results of the simulation obtained.

 

Fig 1.1:- When a supersonic flow encounters a convex corner, it forms an expansion fan, which consists of an infinite number of expansion waves centered at the corner. The figure shows one such ideal expansion fan.

 

Theory                                                                                                                                                                                                                                       

In aerodynamics, the regimes of flight speeds are classified into subsonic, transonic, supersonic & hypersonic as shown in Fig 2.1. They are classified on the basis of their speed when compared to the speed of sound & is known as Mach number. Mach number can be defined as the ratio of the speed of the object by the speed of the sound (speed of sound a = γRT, where γ is the ratio of specific heat, R is the universal gas constant and T is the temperature of the fluid). Mach number is a dimensionless measure of speed common in aerodynamics. The speed of sound depends on the type of medium & temperature of the medium through which it is propagating. For our simplicity, we consider a constant value for the speed of sound at atmospheric conditions which is 330 m/s or 760 mph. Subsonic conditions occur when this Mach number is less than one i.e. (M < 1). For the lowest subsonic conditions, compressibility can be ignored. As the speed of the object approaches the speed of sound, the flight Mach number is nearly equal to one, M = 1, and the flow is said to be transonic. At some places on the object, the local speed exceeds the speed of sound. Compressibility effects are most important in transonic flows. Supersonic conditions occur for Mach numbers greater than one but less than 3 i.e. 1 < M < 3. Compressibility effects are important for supersonic aircraft as shock waves are generated by the surface of the object. For high supersonic speeds i.e. 3 < M < 5, aerodynamic heating consideration also becomes very important for aircraft designing. These supersonic speeds result in Sonic Boom as shown in Fig 2.2. A sonic boom is a loud, thunder-like noise heard by a person on the ground when an aircraft flies overhead at supersonic speeds. The shock wave forms a cone of pressurized air. A sharp release of pressure after the build-up of a shock wave is heard as a sonic boom (similar to the sharp release of pressure when a pin pops a balloon and makes a loud noise.). For speeds greater than five times the speed of sound, M > 5, the flow is said to be hypersonic. At these speeds, some of the energy of the object goes into exciting the chemical bonds which hold together the nitrogen and oxygen molecules of the air. The Space Shuttle re-enters the atmosphere at high hypersonic speeds of M ~ 25. Under these conditions, the heated air becomes an ionized plasma of gas and the spacecraft must be insulated from the high temperatures.

As an aircraft moves through the air, the air molecules near the aircraft are disturbed and move around the aircraft. If an aircraft passes at a low speed, typically less than 250 mph, the density of the air remains constant. But for higher speeds, some of the energy of the aircraft goes into compressing the air and locally changing the density of the air. This compressible effect alters the amount of resulting force on the aircraft. The effect becomes more important as speed increases, near and beyond the speed of sound. Small disturbances in the flow are transmitted to other locations isentropically or with constant entropy. But a sharp disturbance generates a Shock wave that affects both the lift and drag of an aircraft.

 

                                                        Fig 2.1:- Different flight regimes in aerodynamics.

 

 

                                 Fig 2.2:- Jet fighter planes with conical shock waves made visible by condensation.

 

 

                  In physics, a shock wave or shocks, are essentially non-linear waves that propagate at supersonic speeds. The speed of a shock wave is always greater than the speed of sound in the fluid and decreases as the amplitude of the wave decreases. When the shock wave speed equals the normal speed, the shock wave dies and is reduced to an ordinary sound wave. Like any ordinary wave, a shock wave carries energy and can propagate through a medium. Physically the occurrence of a shock wave is always characterized in a fluid flow by instantaneous changes in pressure, velocity, and temperature. This means that the region between the vehicle and the shock wave usually referred to as the shock layer will be a region of high pressure, temperature, and density compared to the freestream flow conditions. In other words, whenever a fluid streamline crosses a standing shock wave, an instantaneous increase in static pressure, temperature and density are observed with a substantial decrease in the velocity of the flow. This sudden change in the medium of propagation is one of the unique features that characterize the shock wave. Because of the sudden changes in flow properties, shockwaves are often regarded as discontinuities with highly localized reversibility.

 

                              Fig 2.3:- Representation of Wavefronts for subsonic, transonic & supersonic speeds.

 

It is well known that any disturbance in a given flow always travels at local acoustic velocity. For example, let us consider a subsonic flow scenario. If we generate a whistling sound from point A this disturbance reaches the flow always at the local acoustic speed. If we assume that the local temperature of the atmospheric air is 300 K, then the speed of sound will be ~ 340 m/s. Let us say the velocity of flow is ~ 30 m/s (subsonic), then the flow will have prior knowledge of the disturbance wave. On the other hand, when the velocity of flow is supersonic say 400 m/s, then the disturbance waves will never be able to catch up with the flow, and all the disturbance wavefronts will coalesce to form a shock wave. In fact, the angle formed by these compression wavefronts is a unique signature of the flow and in some sense, shock waves are the means of information transfer in supersonic flow. Since the information transfer can take place through acoustic disturbance, there will be no shock waves in the subsonic flow regime.

                                                                       

 Fig 2.4:- Different aerodynamic flow regimes indicating the presence of shock waves in the flow field around generic bodies.

 

Since the flow in subsonic speeds has prior knowledge of the disturbance, the streamlines diverge and go around the body consistent with the continuum principle. However, in the supersonic flow because of the presence of shock waves ahead of the body, the velocity of flow suddenly reduces to subsonic values after the shock wave and this in some sense informs the fluid streamlines about the presence of the body. In supersonic flow, shock waves perform the dual function of information transfer as well as a mechanism by which the flow can turn towards itself.

 

Classification of shock waves are:-

 

A. Normal  Shocks

 

If the shock wave is perpendicular to the flow direction it is called a normal shock. Some equations which describe the change in flow variables for flow across a normal shock are written below. The equations presented here were derived by considering the conservation of mass, momentum & energy for a compressible gas while ignoring viscous effects. The equations have been further specialized for a one-dimensional flow without heat addition.

Fig 2.6

 

 

A normal shock occurs in front of a supersonic object if the flow is turned by a large amount and the shock cannot remain attached to the body. The detached shock occurs for both wedges and cones. A normal shock is also present in most supersonic inlets. Across the normal shock the flow changes from supersonic to subsonic conditions. Since gas turbine engines operate under subsonic conditions, it is necessary to introduce a normal shock in the inlet compression system. Normal shocks also are generated in shock tubes. A shock tube is a high-velocity wind tunnel in which the temperature jump across the normal shock is used to simulate the high heating environment of spacecraft re-entry.

 

Across the normal shock wave, the Mach number decreases to a value specified as M1 & can be expressed mathematically as:

(M1)^2 = [(γ - 1) * M^2 + 2] / [2 * γ * M^2 - (γ - 1)]

The total temperature T(t) across the shock is constant,

T(t1) / T(t0) = 1

The static temperature T increases in zone 1 to become:

T1 / T0 = [2 * γ * M^2 - (γ - 1)] * [(γ - 1) * M^2 + 2] / [(γ + 1)^2 * M^2]

The static pressure p increases to:

p1 / p0 = [2* γ * M^2 - (γ - 1)] / (γ + 1)

And the density r changes by:

r1 / r0 = [(γ + 1) * M^2 ] / [(γ -1 ) * M^2 + 2]

The total pressure P(t) decreases according to:

pt1 / pt0 = {[( γ + 1) * M^2 ] / [(γ - 1) *M^2 + 2]}^[ γ /( γ -1)] * {( γ + 1) /[2 * γ * M^2 - (γ - 1)]}^[1/( γ - 1)]

Where γ or gama is the ratio of specific heats (Cp/Cv) & M is upstream Mach number. For air, γ = 1.4

The right-hand side of all these equations depend only on the free stream Mach number. So, knowing the Mach number, we can determine all the conditions associated with the normal shock

. 

 

 

A.                                                          B.  Oblique Shocks

 

When a shock wave is inclined to the flow direction it is called an oblique shock. The equations which describe the change in flow variables for flow across an oblique shock are expressed below. The equations presented were derived by considering the conservation of mass, momentum, and energy for a compressible gas while ignoring viscous effects. The equations have been further specialized for a two-dimensional flow without heat addition. The equations only apply for those combinations of free-stream Mach number and deflection angle for which an oblique shock occurs. If the deflection is too high, or the Mach too low, a normal shock occurs. For the Mach number change across an oblique shock there are two possible solutions; one supersonic and one subsonic. In nature, the supersonic ("weak shock") solution occurs most often. However, under some conditions the "strong shock", a subsonic solution is possible.

 

 

Fig 2.7

 

Oblique shocks are generated by the nose and by the leading edge of the wing and tail of a supersonic aircraft. Oblique shocks are also generated at the trailing edges of a subsonic aircraft as the flow is brought back to free stream conditions. Oblique shocks also occur downstream of a nozzle if the expanded pressure is different from free stream conditions. In high-speed inlets, oblique shocks are used to compress the air going into the engine. The air pressure is increased without using any rotating machinery. From Fig 2.7 we observe, a supersonic flow at Mach number M approaches a shock wave that is inclined at angle s. The flow is deflected through the shock by an amount specified as the deflection angle - a. The deflection angle is determined by resolving the incoming flow velocity into components parallel and perpendicular to the shock wave. The component parallel to the shock is assumed to remain constant across the shock, the component perpendicular is assumed to decrease by the normal shock relations. Combining the components downstream of the shock determines the deflection angle.

 

Across the shock wave the Mach number decreases to a value specified as M1:

M1^2 * sin^2(s -a) = [(γ-1)M^2 sin^2(s) + 2] / [2 * γ * M^2 * sin^2(s) - (γ - 1)]

The total temperature across the shock is constant, but the static temperature T increases in zone 1 to become:

T1 / T0 = [2 * γ * M^2 * sin^2(s) - (γ - 1)] * [(γ -1) * M^2 * sin^2(s) + 2] / [(γ + 1)^2 * M^2 * sin^2(s)]

The total pressure P(t) decreases according to:

pt1 / pt0 = {[(γ + 1) * M^2 * sin^2(s)]/[(γ-1)*M^2 * sin^2(s) + 2]}^[γ/((γ-1)] * {(γ+1)/[2 * γ * M^2 * sin^2(s)-(γ-1)]}^[1/(γ-1)]

The static pressure p increases to:

p1 / p0 = [2 * γ * M^2 * sin^2(s)-(γ -1)] / (γ + 1)

And the density r changes by:

r1 / r0 = [(γ + 1) * M^2 * sin^2(s)] / [(γ -1) * M^2 * sin^2(s) + 2]

 

The right-hand side of all these equations depends only on the free-stream Mach number and the shock angle. The shock angle depends in a complex way on the free-stream Mach number and the wedge angle. So, knowing the Mach number and the wedge angle, we can determine all the conditions associated with the oblique shock. 

 

               C. Bow Shocks

A bow shock, also called a detached shock or normal shock, is a curved propagating disturbance wave characterized by an abrupt, nearly discontinuous, change in pressure, temperature & density. It occurs when a supersonic flow encounters a body around which the necessary deviation angle of the flow is higher than the maximum achievable deviation angle for an attached oblique shock. Then, oblique shock transforms in a curved detached shock wave. As bow shocks occur for high flow deflection angles, they are often seen forming around blunt bodies, because of the high deflection angle that the body imposes to the flow around it. The thermodynamic transformation across a bow shock is non-isentropic and the shock decreases the flow velocity from supersonic velocity upstream to subsonic velocity downstream. The bow shock significantly increases the drag in a vehicle traveling at a supersonic speed. This property was utilized in the design of the return capsules during space missions such as the Apollo Program, which need a high amount of drag in order to slow down during atmospheric re-entry.

 

For this project, we are going to concentrate only on supersonic expansion fan, technically known as the Prandtl-Meyer expansion fan which is a 2d simple wave in a centered expansion process that occurs when a supersonic flow turns around a convex corner. The fan consists of an infinite number of Mach waves, diverging from a sharp corner. When a flow turns around a smooth and circular corner, these waves can be extended backward to meet at a point. Each wave in the expansion fan turns the flow gradually (in small steps). It is physically impossible for the flow to turn through a single "shock" wave because this would violate the 2^nd law of thermodynamics. Across the expansion fan, the flow accelerates (velocity increases) and the Mach number increases, while the static pressure, temperature, density decrease. Since the process is isentropic, the stagnation properties (e.g. the total pressure and total temperature) remain constant across the fan.

 

 

Boundary conditions are a required component of the mathematical model to solve ODE & PDE’s. When solving the Navier-Stokes equation and continuity equation, appropriate initial conditions and boundary conditions need to be applied. A boundary value problem has conditions specified at the extremes ("boundaries") of the independent variable in the equation whereas an initial value problem has all of the conditions specified at the same value of the independent variable (and that value is at the lower boundary of the domain, thus the term "initial" value). A boundary value is a data value that corresponds to a minimum or maximum input, internal, or output value specified for a system or component.

o   Dirichlet boundary condition specifies the value of a variable at a boundary. It may also be called a fixed condition. Therefore, it is used when a value needs to be imposed at a boundary. For example, a no-slip condition in fluid mechanics is a Dirichlet condition because it sets the value of the velocity to zero. In solid mechanics, prescribing a certain load or displacement, and in heat transfer—setting the temperature at a surface are other examples of this type of condition.

o   Neumann boundary condition specifies the derivative of a variable at a boundary. It is used rather than the actual value when a certain rate of change in value needs to be imposed for a variable. A common example in fluid mechanics is the fully developed condition at an outlet where the gradient of flow variables is set to zero. Traction conditions in solid mechanics and insulated surfaces in heat transfer are other examples.

o   Robin boundary condition specifies a combination of the value and the derivative of a variable at a boundary. Therefore, these conditions are suited to more complex behavior. For example, in heat transfer, Robin conditions are used to model Newton’s law of cooling, where heat flux (derivative) is proportional to temperature (value). However, Robin conditions are not the only class of boundary conditions that utilize value and derivative to specify a condition.

o   Some other boundary conditions besides Dirichlet, Neumann, and Robin that exploit certain flow features. For example, in scenarios where the domain features spatial symmetry, the size of the domain, and therefore the computational effort, is reduced by using symmetry boundaries. Similarly, periodic conditions are used to simulate periodicity effects without having to run computations on the full domain.

Thus, boundary conditions can be used to effectively gain an advantage.

 

 

 

 GOVERNING EQUATIONS                                                                                                                                                                                             

 

 The solution for a CFD simulation is obtained by solving the Navier-Stokes (NS) equations. The NS equations comprise of 5 equations, namely the continuity, momentum (3 equations), and energy equations. The continuity and momentum equations are solved for all flow problems with the energy equation being optional and can be used only when heat transfer is taking place to save computation time.

 

For our purpose, we will be assuming that there is no heat transfer and that the temperature is constant, leaving us with 4 equations now, the continuity and 3 momentum equations, one for each axis direction. 

 

 

 An important property of Converge is that any geometry created or exported is assumed to be made of triangles. In boundary flagging, we group these triangles to a particular ‘Boundary’ & assign these boundaries to a ‘Volumetric region’. Defining these boundary conditions is a fundamental step in any CFD simulation as this helps solve the NS equations. Converge creates and exports ‘input’ file of these complex governing equations which are then solved by CYGWIN for running the simulation.

 

                                                   Prandtl-Meyer Shock Problem Simulation                                                                                           

 

 

The pre-modeled duct was designed using Solidwork’s & imported into Converge Studio as an STL file. After importing the STL file of the duct into converge we convert the dimensions. Since the geometry was created in solid works, the dimensions of the geometry are expressed in millimeters or ‘mm’ but Converge software package uses S.I unit systems i.e. ‘meters’ & thereby assumes that the geometry imported is in meters as shown in Fig 1.1. To fix this we use the geometry editing tools available in converge, we ‘Transform’ the entire geometry by a uniform scale factor of 0.001 as shown in Fig 1.2. Once, the geometry is scaled down we run a ‘Diagnostic’ test to check for any anomaly like intersection errors, nonmanifold problems, open edges, etc within the geometry contour. If no errors, which is denoted by “green checks”, then we proceed to check for the orientations of the “Normal’s” in the geometry.

    Fig 1.1 & 1.2 :- Transforming dimensions

 

 

Every geometry will have a normal vector which is perpendicular to the geometry, for this problem we use the “Normal Toggle” option to check for the direction of the normal’s. If the normal’s are pointing outside of the geometry, then it’s essential to transform these normals to point inside the geometry, where the fluid flow will occur as shown in Fig 1.3 &1.4.

Fig 1.3 & 1.4 :- Transforming Normal’s

 

 In Converge, all the surface information is stored in the form of 3 entities vertices, edges & triangles. Therefore, any geometry created or exported from other CAD software (like in our case) into converge are converted into triangles & that is the fundamental entity in converge, where every part of the geometry is assumed to be a triangle. The process of grouping these triangles into boundaries is known as boundary flagging. Each boundary is assigned with a distinct ID as shown in (Fig 1.5). For our geometry, we have 5 distinct surfaces Inlet, Outlet, Front 2D, Back 2D & Top & Bottom Walls.

Fig 1.5:-  Boundary flagging of the geometry

 

   2. Case Setup

  Having edited the geometry, the next stage of the process is to set up the simulation parameters. This step involves specifying certain properties to capture the physics of the problem perfectly, such as materials used, simulation time parameters, solver parameters, initial & boundary conditions to solve the Navier Stokes equation which are PDE’s, defining body forces, type of fluid flow, species involved & grid size of the mesh.

We start setting up the time-based general flow simulation by first defining the Materials involved in the simulation. Since we are simulation fluid (air) flow inside a channel we choose ‘gas simulation’ and our choice of fluid for this project is ‘air’, hence we choose that as our pre-defined mixture. For gas simulation parameters we will be using Redlich-Kwong equation of state with a critical temperature of 133K & critical pressure of 3770000pa. For global transport parameters, we use the default values of ‘turbulent Prandtl number’ and ‘Schmidt number’ which is 0.90 & 0.78 respectively. Since the fluid flow is air our ‘species’ will be a mixture of Oxygen (O2) & Nitrogen (N2) with a chemical composition of 23% & 77% respectively.

 

The next parameter for our case setup is the Simulation Parameters, which include the ‘run parameters’, ‘simulation time parameters’ & the ‘solver parameters’. For the elbow duct with throttle flow problem, we will be running compressible gas flow using Steady-state solver at full hydrodynamics simulation mode, since our geometry is simple & has no moving part & we also require to solve the NS equations. In addition to this, we will be using ‘density-based’ PISO Navier strokes solver. As for simulation time parameters, we will be running the simulation for 25000 cycles with an initial & minimum time step as (1e-07).

 

As earlier, after designing the geometry, we grouped the triangles and assigned them to a particular boundary ID. Similarly, in Boundary Condition, we group the 5 boundaries & assign each of them to a volumetric region. To do so we add a volumetric region from Initial Conditions & Events & assign the volumetric region (Region 0) to each of the boundaries. These volumetric regions are required to set up the initial conditions for solving PDE’s. Initial conditions are assigned at the volume whereas boundary conditions are given at boundaries. The Inlet & Outlet boundaries are inflow & outflow boundaries respectively where the air will be flowing inside the duct at supersonic speeds of 680 m/s with temperature boundary condition specified as 286.1 K. At the outlet pressure & velocity boundary condition will be following Neumann boundary conditions with a specified temperature of 300k with Air being the species, if backflow occurs. Front2D & Back 2D are just 2D boundaries while Top & Bottom walls are stationary walls with fixed surface movement following the ‘slip’ conditions & a temperature boundary condition of ‘zero normal gradients’ which is Neumann boundary condition.

 

 

The final parameter that we need to setup is defining the Mesh that is the Base grid size. For our case mesh grid size will be kept constant of dx=dy=dz=0.8. In addition, we will also be using temperature-based Adaptive Mesh Refinement (AMR) a technique to monitor the temperature at certain locations where a shock is expected to occur. Adaptive Mesh Refinement was provided with the maximum embedding of 2 layers & a sub-grid scale of 0.05, 0.004 & 0.002. Sub-grid scale (SGS) is a parameter that will be used by converge for refining the mesh at a specific location. Converge, monitors, the curvature of the fluid property provided in the AMR and refines the mesh when the property curvature variation between the consecutive grid is more than the defined SGS value. We will be running the simulation for all 3 different cases of SGS keeping the base mesh size constant.

 

                                                                             3.  Post Processing

 

The function of Converge studio is to set up the simulation & then create several input files that are then exported to a particular folder. To run these input files we use CYGWIN, a command-line interface that reads these inputs files & solves the complex PDE’s of governing equation to generate several output files. These output files are then post-converted into 3D output files which are readable files for ParaView.

 

                                                                                    4. Results

 

Case 1   Base Mesh Grid Size of 0.8 with Sub Grid-Scale temperature of 0.05                    

                                   

A.  Mesh Generated with Adaptive Mesh Refinement

 

 

 

B.  Temperature, Pressure, Velocity Contour with no AMR

       

 

C.  Animation For Temperature, Pressure, Velocity with AMR

 

       

       

 

 

   

C                                 Case 2   Base Mesh Grid Size of 0.8 with Sub Grid-Scale temperature of 0.04                   

                                   

A.  Mesh Generated with Adaptive Mesh Refinement

 

 

 

B.  Temperature, Pressure, Velocity Contour with no AMR

 

       

 

 

C.  Animation For Temperature, Pressure, Velocity with AMR

 

     

     

     

 

 

 

                                  Case 3   Base Mesh Grid Size of 0.8 with Sub Grid-Scale temperature of 0.02                  

                                   

 

A.  Mesh Generated with Adaptive Mesh Refinement

 

 

 

B.  Temperature, Pressure, Velocity Contour with no AMR

 

       

 

 

    Temperature, Pressure, Velocity with AMR

 

     

     

ht 

 

 

 

                                                                                     Observation

 

From the data collected we observe that the number of mesh grid cells generated varies for each case. The finer the mesh the better, as it allows us to capture the physics of the expansion fan-generated due to the sharp convex corner. The total number of mesh generated for an initial condition is equal for each of the 3 cases, however due to variable Adaptive Mesh Refinement used we observe a difference in the total number of cells generated in each case as fluid starts entering the domain. It is quite evident that as we decrease the value of temperature-based Sub Grid-Scale the number of cells generated increases & vice versa equally holds true. Thereby we can conclude that that case 3 has the greatest number of cells generated & capture the most physics while case 1 has the least number of cells generated near the shocked fan. Due to these criteria, we observe a much prominent temperature, pressure & velocity variation for case 3. We can observe a rise in velocity & a drop in temperature & pressure of the fluid across the expansion fan indicating a region of pressure & temperature difference. Due to lower SGS criterion we can observe a smother resolution at the tip of expansion fan near the convex curve more prominently for case 3 than when compared to case 2 & case 1.