Week 5: Prandtl Meyer Shock problem

Aim: Prandtl Meyer Shock problem

Objective: 

1. Shock flow boundary conditions

Do a literature search on what BC's are typically used for shock flow problems

2. What is a shock wave?

In your own words, describe the physics behind shock waves

3. Effect of SGS parameter on shock location

In the Prandtl Meyer shock problem, look at the effect of SGS temperature value on cell count and the shock location

It is expected to show the mesh generation using ParaView.

Also, create animation and upload it on Youtube and provide the Youtube link for animation in your report.

Introduction and Theory

A centered, Prandtl Meyer expansion wave is illustrated below. The type of flow highlighted in the below figure is inviscid, supersonic flow moving over the surface. A supersonic flow (`M>1`) is expanded around the expansion corner which creates an infinite number of weak. The leading edge makes an angle `mu_1` with respect to the upstream flow direction. The trailing edge makes an angle `mu_2` with respect to the downstream flow direction. `mu_1` and `mu_2` are Mach angles which are defined by 

Equation 1

`mu_1=sin^-1(frac{1}{M_1})`

Equation 2

`mu_2=sin^-1(frac{1}{M_2})`

`M_1` and `M_2` are upstream and the downstream Mach number.

As observed from the figure the Mach number increases after the expansion corner whereas temperature, pressure, density decreases after the expansion corner.

 

 

 

The flow around the expansion corner is assumed to be isentropic in nature. The flow in front of the expansion corner is uniform and parallel to the wall with Mach number `M_1` and the flow behind the expansion corner is uniform and parallel to the wall with Mach number `M_2`. The only exception to the flow is at the corner where the flow is discontinuous in nature. For the calorically perfect gas, there is an exact, analytical solution behind the expansion corner which is given by 

Equation 3

`f_2=f_1+theta`

where `f` is the Prandtl Meyer function and `theta` is the flow deflection angle. For a calorically perfect gas, the Prandtl Meyer function depends on `M` and `gamma` which is given by 

Equation 4

`f=sqrt(frac{gamma+1}{gamma-1})tan^-1sqrt(frac{gamma-1}{gamma+1}(M^2-1))-tan^-1sqrt((M^2-1))`

The inlet Mach number (`M_1`) is known, now from Equation 4 `f_1` can be calculated. The value of `theta` is already known, therefore `f_2` can be calculated by Equation 3. Now using `f_2`, `M_2` can be calculated by Equation 4. Thus by using the values of `M_1`, and `M_2` the pressure, temperature, density can be calculated behind the expansion corner.

Equation 5

`p_2=p_1{frac{1+[frac{gamma-1}{2}]M_1^2}{1+[frac{gamma-1}{2}]M_2^2}}^frac{gamma}{gamma-1}`

Equation 6

`T_2=T_1{frac{1+[frac{gamma-1}{2}]M_1^2}{1+[frac{gamma-1}{2}]M_2^2}}`

Equation 7

`rho_2=frac{p_2}{RT_2}`

Therefore the exact or analytical solution behind the expansion corner is determined by using these equations.

Shock wave

Shock waves are very thin regions in the gas flow field typically on the order of a few molecular mean free paths where flow properties are discontinuous across its thickness. A shock is a traveling wave, i.e. it moves relative to the fluid faster than the local speed of the sound.

Normal shock wave 

One special type of shock is the normal shock where the shock wavefront is perpendicular to the free stream flow and is observed in some internal and external supersonic flow applications.

Oblique shock

An oblique shock wave is a shock wave that unlike a normal shock is inclined with respect to the incident upstream flow direction. It will occur when a supersonic flow encounters a corner that effectively turns the flow into itself and compresses.

Let us now Mathematically describe the flow through an oblique shock wave. Consider an oblique shock wave oriented at an angle `beta` to the incoming supersonic flow which is at a Mach number of `M_1`, the equivalent velocity `V_1`. This velocity can be represented as a superposition of two components, a velocity component that is normal to the shock `u_1` and the velocity component parallel to the shock `Vt_1`. The relationship between `V_1` and the velocity component is shown below

`V_1=sqrt(u_1^2+Vt_1^2`

The two velocity components are related to the oblique shock wave angle `beta` which is given by 

`beta=tan^-1(frac{u_1}{Vt_1})`

Applying continuity equation 

`rho_1u_1=rho_2u_2`

The tangential component of the momentum equation with a steady flow and nobody force  assumption

`(-rho_1u_1)Vt_1+(rho_2u_2)Vt_2=0`

which gives,

`Vt_1=Vt_2`

Therefore the shock oblique parallel component is the same upstream and downstream of an oblique shock.

The normal momentum equation is given by 

`p_1+rho_1u_1^2=p_2+rho_2u_2^2`

The energy equation is given as 

`h_1+frac{u_1^2}{2}=h_2+frac{u_2^2}{2}`

It can be concluded that the governing equation of a flow-through an oblique shock wave is dependent on the shock normal velocity `u_1` and `u_2`

For an oblique shock, the variation of the flow properties is given by the following relationship

`frac{T_2}{T_1}=frac{a_2^2}{a_1^2}=1+frac{2(gamma-1)(gammaM_1^2sin^2beta+1)}{(gamma+1)^2(M_1^2sin^2beta)}(M_1^2sin^2beta-1)`

`frac{p_2-p_1}{p_1}=frac{2gamma(M_1^2sin^2beta-1)}{gamma+1}`

`frac{rho_2}{rho_1}=frac{(gamma+1)M_1^2sin^2beta}{(gamma-1)M_1^2sin^2beta+2`

If `beta=frac{pi}{2}` the exact normal shock wave relation will be recovered. This indicates that the normal shock waves are just the unique case of oblique shock waves.

The Mach number downstream `M1_2` is given by the following equations 

`M_2^2sin^2(beta-theta)=frac{1+frac{gamma-1}{2}M_1^2sin^2beta}{gammaM_1^2sin^2beta-frac{gamma-1}{2}`

From the velocity triangle upstream and downstream of the shock 

`tanbeta=frac{u_1}{V_t}`

`tan(beta-theta)=frac{u_2}{V_t}`

Therefore combining this equation with the continuity equation and the density ratio the following relationships can be obtained 

`theta-beta-Mach` Number Relationship as 

`tantheta=2cotbetafrac{M_1^2sin^2beta-1}{M_1^2(gamma+cos2beta)+2}`

Literature review on the Boundary Conditions

Boundary conditions: Boundary conditions are the necessary constraints that are useful for solving ordinary and partial differential equations. The wrong selection of boundary conditions may lead to divergence or convergence to a wrong solution. The types of BC are listed below

Dirichlet BC: The Dirichlet BC specifies the value that the unknown functions need to take along the boundary of the domain.

Example: If `u(x,t)` represents the displacements of the vibrating string and its ends are fixed at `x=0` and `x=h` such that `u(0,t)=0` and `u(h,t)=0`, the condition is known as Dirichlet BC.

Neumann BC: The Neumann BC specifies the value of the derivative that the unknown function is going to take along the direction normal to the boundary which is written as `frac{deltau}{deltan}` which is defined as,

`frac{deltau}{deltan}``=(frac{deltau}{deltax_1},frac{deltau}{deltax_2},...).n`

Example: If `u(x,t)` represents the temperature distribution in a rod of length L. Assuming that the rods are perfectly insulated at x=0, x=L, the heat flux at these points are zero, so it follows that the appropriate BC as per Fourier Law of Heat Conduction is 

`frac{deltau(0,t)}{deltax}=0` and `frac{deltau(L,t)}{deltax}=0` is the Neumann BC.

Robin BC: The Robin BC is a weighted combination of Dirichlet BC and the Neumann BC which is expressed below

`(alphau+betafrac{deltau}{deltan})deltaD=f(x)_(deltaD)`

where `alpha` and `beta` represent weights.

Example: If `u(x,t)` represents the temperature distribution in a rod of length L. The rods are poorly insulated at their end therefore the BC is given by 

`u(0,t)+frac{deltau(0,t)}{deltax}=0`

`u(L,t)+frac{deltau(L,t)}{deltax}=u_0`

Mixed BC: Mixed BC refers to the condition where Dirichlet is prescribed in some parts of the boundary and Neumann is prescribed in some parts of the boundary.

Example: If u(x,t) represents the temperature distribution in a rod of length L. The one end of the rod is immersed in a water bath with constant temperature and the other end of the rod is connected to the heater with a constant heat transfer rate then Mixed BC takes the form of  

`u(0,t)=A`

`frac{deltau(L,t)}{deltax}=B`

Outlet BC is Neumann

Explanation

In the present challenge the specified value of temperature, pressure, and velocity is unknown only the value of the normal gradient is known which is zero. Therefore, at the outlet Neumann BC is specified.

Case-setup

The 2d_supersonic_flow. STL file is loaded into the Converge-CFD set-up. 

The diagnosis dock is selected from the "View" option and the "Find" option is selected. It is found that there are no ("Intersections(0)","Nonmanifold Problems(0)","Open Edges(0)", "Overlapping Tris(0)","Normal Orientation(0)","Isolated tris (0)".

The "Normal Toggle" is selected from the ribbon, upon selecting the normal toggle it was observed that the normal is pointing outside the volume, technically the normal should point inside the volume where the fluid is flowing. Therefore in the geometry dock, the "Transform" option is clicked. The "Normal" tab is selected from the geometry dock and one of the triangles containing the normal pointing outward is selected. The "Apply" option is selected for removing the normals.

The normal direction of 2D boundaries must be aligned in Z-direction otherwise during post-processing, it will give an error "the Normal direction of TWO_D boundaries must be aligned in Z-direction".

 

The "Case Setup Dock" is selected from the "View" and "Begin Case Setup" is executed.

The "Time" based "Application" type is selected. The "Material" is selected and the predefined mixture is selected as air. The "Reactant mechanism" is checked off because it is not a combustion problem. The species is selected and the "Apply" is clicked. Under "Gas simulation" the Equation of state is selected as "Redlich Kwong", the critical temperature is 133K and the critical pressure is 3770000Pa. The Turbulent Prandtl number is 0.9 and the Turbulent Schmidt number is 0.78 under Global Transport parameters. The `O_2` and the `N_2` are selected as species. Under "Run parameters" the "Steady-state solver" is selected. The "Temporal type" is locked for the steady-state. The "simulation mode" is selected as "Full Hydrodynamic" because the geometry is simple so while creating the mesh inside the geometry it solves the NS equation as well, but if the geometry is complex "No hydrodynamic solver" is selected as such if there is an error it will point out immediately while creating the mesh, otherwise if the hydrodynamic solver is selected then it will be very tedious to identify the error in the simulation case set up.

Under "Simulation time parameters" end time is chosen as 25000 cycles (steady-state). The initial and minimum time step is 1e-09 s. The maximum time step is 1s, and the maximum convection CFL limit is 1 and the rest are default values.

The "Pressure" "Equation" is selected under the "Solver parameter"  and the "Preconditioner" is selected to "None". The "Maximum convection CFL limit (final stage)" is set to 0.5. The "Density-based solver" is selected it is suitable for even low Mach numbers to avoid extrapolating temperature calculations and the PISO coupling scheme is selected.

Base-grid

 

Adaptive Mesh Refinement

Adaptive mesh refinement is a technique that is used to refine the grids automatically based on fluctuating and moving conditions such as temperature or velocity. The feature is used to refine the grid in the flow region of interest such as flame propagation or high velocity without slowing down the simulation with a globally refined grid. There are two types of AMR, mainly sub-grid scale-based and value-based.  The AMR type and criteria are chosen in such a way that the embedding is added where the flow field is most unresolved in this case being expansion corner.

 

 The curvature in temperature is monitored, the curvature means `frac{delta^2T}{deltax^2}`, in AMR dialogue box subgrid criterion is set to 0.05, 0.04, 0.03. Whenever `frac{delta^2T}{deltax^2}` is `gt` than 0.05, 0.04, 0.03 K, the AMR is going to refine the mesh near the shock region.

Boundary selection

The boundary conditions for the different named selections are tabulated below

Results  

Case-1

Inlet velocity 678m/sec

Mach number = 1.994

Type of flow: Supersonic 

Grid size: 0.8m

Subgrid scale: 0.05K

Mesh

 

Contour plots

Velocity

Temperature

Pressure

Density

Animation file

Velocity

 

Temperature

Pressure

Density

Line plots

Velocity 

Temperature

Mach number

 

 

Pressure

Density

Total cell count

Subgrid scale: 0.04K

Contour plots

Velocity 

Temperature

Pressure

Density

Animation file

Velocity 

Temperature

Pressure

Density

Line plots

Velocity 

Temperature

Mach number

Pressure

Density

Total cell count

 Subgrid scale: 0.03K

Contour plots

Velocity

Temperature

Pressure

Density

Animation file

Velocity

Temperature

Pressure

Density

Line plots

Velocity

Temperature

Mach number

Pressure

Density

Total cell count

Case-2

Inlet velocity 100m/sec

Mach number = 0.294

Type of flow: Subsonic 

Grid size: 0.8m

Subgrid scale: 0.05K

Contour plots

Velocity

Temperature

Pressure

Density

Animation file

Velocity

Temperature

Pressure

Density

Line plots

Velocity

Temperature

Mach number

Pressure

 

Density

Total cells

Subgrid scale: 0.04K

Contour plots

Velocity

Temperature

Pressure

Density

Animation file

Velocity

Temperature

Pressure

Density

Line plots

Velocity

Temperature

Mach number

Pressure

Density

Total cell count 

 

Subgrid scale: 0.03K

Contour plots

Velocity

Temperature

Pressure

Density

Animation file

Velocity

Temperature

Pressure

Density

Line plots

Velocity

Temperature

Mach number

Pressure

 

Density

 

 Total cell count

 Summary of all cases 

Conclusion

• The steady-state simulation of Prandtl Meyer shock waves has been performed successfully with Converge-CFD and Cygwin command prompt terminal.
• The expansion waves were visualized for supersonic and subsonic flow. In the case of Supersonic flow, the infinite Mach waves were produced at the expansion corner showing the discontinuity the flow field exhibits. In the case of subsonic flow, the waves were not clear and visible.
• For the supersonic flow, the pressure, temperature, density decreases downstream, and the Mach number increase downstream.
• The Adaptive Mesh Refinement technique was used to refine the mesh grids near the expansion corner. The temperature type AMR was used and the subgrid criterion was set to 0.05K, 0.04K, 0.03K. The Converge typically refines the grid near the expansion corner where the flow field is not resolved. When the `frac{delta^2T}{deltax^2}gt`0.05, 0.03, 0.03 the mesh is refined without globally refining the grids thus saving a lot of time.
• It was also found that when the subgrid criterion decreases from 0.05K to 0.03K the total cell count increases.
• The shock location is more visible when the mesh is refined near the expansion corner i,e for SGS 0.03K otherwise for the higher SGS the shock location is blurred. 

References

Computational Fluid Dynamics-THE BASICS WITH APPLICATIONS BY JOHN D. ANDERSON, JR

Oblique Shock Wave: [https://www.youtube.com/watch?v=QCMwo4pAVGM&t=13s]

Adaptive Mesh Refinement: Converge_3.0_Manual.PDF