PRANDTL MEYER 2D SUPERSONIC FLOW SIMULATION USING CONVERGE FLOW SOLVER

OBJECTIVE

The objective of this project is to set up a Prandtl Meyer Shock problem and solve using Converge Flow Solver with different inlet boundary and initial condition:

• Subsonic at 100 m/s
• Supersonic at 680 m/s

The overview of this project is to understand the following:

• What is Shock Wave?
• About shock flow boundary conditions
• Effect of SGS parameter on shock location
• Effect of SGS temperature value on cell count and shock location

Shockwaves

Shockwave or shock is a type of propagating disturbance that moves faster than the local speed of sound in the medium. Like an ordinary wave, a shockwave carries energy but is characterized by a discontinuous (or instantaneous) change in pressure (static), temperature and density by a large amount. Because, shockwave does no work, no heat addition while total enthalpy and total temperature are constant. In the case of shockwaves, the flow process is irreversible and entropy increases. Since it is non-isentropic, the total pressure downstream of the shock is always less than the total pressure upstream of the shock.

Types of shockwaves:

1. Normal Shockwave – Shock wave is perpendicular to the flow direction. It occurs in front of supersonic upstream flow and remains attached with the body. Here, upstream flow is always supersonic and downstream flow is subsonic. In general, one-dimensional analysis is performed due to discontinuous change in fluid properties along the direction of flow.

[Reference: https://www.grc.nasa.gov/www/k-12/airplane/normal.html]

Applications:

• Occurs at most of the supersonic inlet
• Generated in shock tubes for testing prototype
• Since gas turbine engines operate under subsonic conditions, a normal shock is introduced in the inlet compression system.

2. Oblique Shockwave – Shock wave is inclined to the upstream flow direction and also changes the direction of downstream flow in such a way that it turns toward itself. It may remain attached to the body at the corners (concave or convex) or pointed bodies (wedge or cones). Here, upstream flow is supersonic while downstream flow can be either supersonic (most often, also called as a weak shock) or subsonic (under some conditions, also called as a strong shock). Here, two-dimensional analysis is performed due to discontinuous change in fluid properties which is a function of both x and y.

[Reference: https://www.grc.nasa.gov/www/k-12/airplane/oblique.html]

Applications:

• To increase air pressure without using any rotating machinery going into the engine. Here, strong shocks are observed.
• Generated by the nose and by the leading edge of the wing and tail of a supersonic aircraft. Here, weak shocks are observed.

3. Bow Shockwave – These shockwaves are curved and formed a few distances in front of the body. Thus, detached shock generated for blunt bodies as well as for sharp bodies at low Mach numbers.

Applications:

• Space return vehicles such as space shuttles
• Bullets
• Formed at the bow (front) of the ship

Prandtl-Meyer Shockwave

Prandtl-Meyer Shockwave is also known as expansion wave or centred expansion. It shows the opposite effect of the oblique shock. An expansion wave is formed when the supersonic flow is turned away from itself. Here, the flow direction changes smoothly as well as flow properties across the expansion fan centred at point ‘A’.

The expansion wave comes from the sharp convex corner inclined at an angle, θ from horizontal along which flow direction changes. Thus, also called an angle of deflection. The expansion fan comprises an infinite number of the weak shockwave (supersonic) originating from point ‘A’ also called a Mach wave. Expansion fan is bounded by forward (M₁) and rearward (M₂) Mach line defined by Mach angle, i.e. μ₁ and μ₂.

`μ_Mn=sin^(-1)⁡(1/M_n )`

The expansion process is considered as isentropic regardless of shockwaves being non-isentropic. Because expansion takes place through a series of Mach waves such that there is no entropy change. At downstream, static properties of flow decreases and Mach number increases, but total pressure remains constant at inlet and outlet.

Mach Number

Mach number is a dimensionless quantity defined as the ratio of the velocity of flow to the local speed of sound. In aerodynamic and fluid dynamic applications, Mach number and Reynolds number are the important parameters related to compressibility (springiness) and viscosity (stickiness).

`'Mach Number at particular temperature'= (Velocity of air)/(Speed of sound)`

Classification of flow based on Mach number –

1. Subsonic: M < 0.8
2. Transonic: M belongs to 0.8 – 1.3
3. Supersonic: M belongs to 1.3 – 5.0
4. Hypersonic: M belongs to 5.0 – 10.0
5. High-Hypersonic: M belongs 10.0 – 25.0
6. Re-entry Speeds: M > 25.0

Boundary Conditions

Boundary conditions are constraints necessary for the solution of the boundary value problem.

For solving the steady-state problem, appropriate boundary conditions are needed. On the other hand, for solving a transient problem, the appropriate initial condition is also required. Boundary conditions are generally classified into three types explained below:

1. Dirichlet boundary condition

The prescribed value for the variable, h(x,y,z,t) is specified at the boundary of the domain. The predefined value may be either constant or vary in space or in time.

2. Neumann boundary condition

The gradient of the variable () is specified at the boundary of the domain where n denotes direction x, y, and z. If the gradient of the variable is zero, then there is no flow across that boundary.

3. Robin boundary condition

It is the linear combination of Dirichlet and Neumann boundary condition applied to a complete boundary given by:

`ah+b* (∂h)/(∂n) = constant`

where,

a and b are real parameters also called as ‘impedance condition’

h is Dirichlet boundary condition

`(∂h)/(∂n)`is Neumann boundary condition

4. Mixed boundary condition

Applying different types of boundary conditions in different parts of the boundary. While in robin boundary condition, single constraint consisting of different boundary conditions applied to complete part of a boundary.

5. Cauchy boundary condition

It consists of two boundary conditions (‘1’ Dirichlet and ‘1’ Neumann) i.e. two constraints applied to complete boundary. While, Robin boundary condition is a single constraint, a linear combination of Dirichlet and Neumann boundary condition applied to complete boundary.

Illustration for Shockwave

To understand the concept of a shockwave. Let’s take an example:

Speed of sound at 286.1 K = 339 m/s

If (case 1),

Inlet Velocity of air at 286.1 K = 100 m/s

Then, Mach Number for 100 m/s air at 286.1 `K = (100/339) = 0.295`

If (case 2),

Inlet Velocity of air at 286.1 K = 680 m/s

Then, Mach Number for 680 m/s air at 286.1 `K = (680/339) = 2.006`

For understanding Prandtl-Meyer shock problem, the subsonic and supersonic flow is simulated into a duct having a sharp corner in the way. In the case of supersonic flow, we must be able to observe an infinite number of weak Mach lines as well as a change in properties in the downstream flow as the per concept discussed in the previous section. While in the case of subsonic flow, shockwave is not generated.

Case Setup

Three standard steps for performing any simulation:

1. Pre-Processing
2. Solver
3. Post-Processing

Pre-Processing

Geometry Import

The geometry can be simply added to CONVERGE interface by dragging-and-dropping the compatible CAD file.

Further, imported geometry was checked for any errors using the diagnosis tool. As the geometry was clean, we proceed to define the boundary. Every other entity in converge is represented by triangles (or discretised into minimum possible triangles), which are flagged to different boundaries.

Boundary Creation

Geometry > Boundary > Flag > Add

For the Prandtl Meyer shock simulation, the complete geometry is defined through six boundary names:

1. Inlet
2. Outlet
3. Front and Back Wall (2D simulation)
4. Top and Bottom Wall (Slip boundaries)

Each boundary is assigned with their respective boundary ID number and defined using selection or grouping of triangles.

Case Setup (Some of the settings are only being mentioned/explained here)

Application:

• General Flow

Material:

• Air (predefined mixture)
• Gas simulation
• No reaction mechanism
• Species (Air comprises of O2 and N2)

Simulation Parameters:

• Steady-state solver
• Start-End time (0.0 – 25000 cycles)
• Initial time step (1.0e-7 s)
• Minimum time step (1.0e-7 s)
• Maximum time step (1.0 s)

Boundary:

Boundary Name	Boundary Type	Boundary Condition
		Neumann (Zero Normal Gradient)	Dirichlet
Inlet	INFLOW	Pressure	Velocity, 100 or 680 m/s
			Temperature, 286.1 K
Outlet	OUTFLOW	Pressure	Temperature, 300 K
		Velocity
Top Wall	WALL	Temperature	Stationery and Slip
Bottom Wall
Front Wall	TWO_D
Back Wall

As per the previously defined boundaries, we are required to provide boundary conditions for the surfaces bounding the entire domain. The boundary condition relevant to a specific boundary has been tabulated above.

In the case of shockwave problem, the flow field in front of the shockwave is called upstream and flow behind shockwave is called downstream. In downstream, the fluid properties changes across the expansion fan. Thus, the fluid properties are unknown and required to define Neumann boundary condition at the outlet.

Regions and initialization:

Initial condition: either 100 m/s (subsonic) or 680 m/s (supersonic)

Based on the case considered to be simulated, the initial condition is defined. If the initial condition is considered to be subsonic, then same inlet velocity is defined in boundary condition and vice-versa. For this subsonic case, we will not be able to generate shockwave.

Turbulence Modelling: Reynolds Averaged Navier-Stokes (RANS) – Default values considered

Base Grid: dx = dy = dz = 0.8

It gives the control to define base grid size, i.e., dx, dy, and dz which is considered to be the input for the discretization of geometry into cells or to create a mesh. Cell count estimation helps to predict cell counts after discretization. Here, we consider the base grid size as 0.8 for the baseline configuration. This base grid is valid for geometry having bounding box dimension as follows:

Length along x-axis = 65 m

Length along y-axis = 45.1525 m

Length along z-axis = 8.151 m

In case the geometry is bigger or smaller with respect to dimensions considered, then we must ensure that we either scale down or scale up the model.

Adaptive Mesh Refinement:

Maximum Cells = 200000

Active Regions = Region 0

AMR type = either based on Temperature or Velocity

Sub-grid criterion = 0.05 K

Embed scale = 2

Based on based grid size control, the cell size has been fixed for the entire domain. But in order to locate the Prandtl-Meyer shockwave, we require much finer mesh as across the expansion fan the discontinuous change in fluid properties is observed i.e. behind the Mach waves in downstream flow.

For this adaptive mesh refinement technique is used that can be either based on temperature or velocity as both parameters undergo a large amount of change. So, the curvature of either temperature or velocity can be monitored for local mesh refinement. Curvature is the second-order derivative of a parameter with respect to space. Whenever the value of curvature at any location exceeds the predefined value, the mesh gets refined or can be said AMR algorithm gets activated. Here, the predefined value of temperature is called a sub-grid scale criterion (SGS). Based on SGS (i.e. 0.05, 0.04, and 0.03), three cases for Prandtl-Meyer shock is simulated and studied its effect on cell count and location of shock. Embedded scale defines the level of refinement based on the base grid for a smoother transition between finer cell to base grid cell given by:

Finer Cells Grid Size `=(Base Grid Cell Size)/2^(Embedded Level) =0.08/2^2 =0.02`

Post Variable Selection: Default

Output files:

Time interval to writing 3D output data files – 500 cycle

The time interval for writing text output – 1 cycle

The time interval for writing restarting output – 1000 cycle

Solver

Using Cygwin, a dos-based solver is used to solve the problem with the help of numerical technique.

The input files directory for channel flow is to be set as a present working directory, through which solver is executed using the below-mentioned command line.

mpiexec.exe -n 4 converge.exe </dev/null> logfile &

This command-line helps to solve using parallel processing with better usage of system cores and also saves ‘logfile’, easy to understand and helpful for future reference.

Post-Processing

Effect of Sub-Grid Scale on Mesh (Cell Count)

The effect of sub-grid scale on cell count has been observed based on temperatures i.e. 0.05, 0.04, and 0.03 Kelvin. The three cases considered with their respective observed data has been tabulated below.

Case No.	Sub-Grid Scale Temperature (K)	Base Grid Size (m)	Based on Cycle History
			Initial Cell Count	Max Cell Count
1	0.05	0.08	4429	7359
2	0.04			13623
3	0.03			18289

The below-shown plot shows the cell count based on the last cycle for three sub-grid scale cases.

From the above-shown plot and tabulated data, it can be inferred that as the sub-grid scale based on temperature is decreased, the cell count increases as AMR gets activated.

The below figure shows the mesh at particular cycle history based on sub-grid scale definition. The effect on a mesh can be observed in the Prandtl-Meyer shock generation region as the AMR algorithm gets activated for different sub-grid scale temperature.

Case 1: SGS 0.05 K

Case 2: SGS 0.04 K

Case 3: SGS 0.03 K

Based on mesh refinement w.r.t. sub-grid scale in Prandtl-Meyer shock generation region, it can be clearly said with lower sub-grid scale temperature more accurately shock can be captured originated from the sharp convex corner and the other solutions related to flow characteristics also get converged. If analysed mesh for case 1, it can be seen that it has not been able to capture the complete expansion fan.

Now, let’s see the effect on flow properties downstream due to Prandtl-Meyer shock with respect to upstream for case 3 i.e. SGS 0.03 K and justify the simulation results with reference to theory of Prandtl-Meyer shock.

MACH NUMBER		VELOCITY
M1 (supersonic) < M2 (supersonic)		V1 (upstream) < V2 (downstream)
DENSITY		PRESSURE
ρ1 (upstream) > ρ2 (downstream)		P1 (upstream) > P2 (downstream)
	TEMPERATURE
	T1 (upstream) > T2 (downstream)

Some of the calculations based on simulation results –

Mach Number of upstream flow (given), M₁ = 2.006

Mach Number of downstream flow (unknown), M₂ = 2.27672

Using Mach angle relation,

`μ=sin^(-1)⁡(1⁄M)`

Forward Mach Angle = 29.90 degree

Backward Mach Angle = 26.05 degree

Using the Prandtl-Meyer function,

`ϑ(M)=√((γ+1)/(γ-1)) tan^(-1)⁡√((γ-1)/(γ+1) (M^2-1) )-tan^(-1)⁡√(M^2-1)`

`θ=ϑ(M_2 )-ϑ(M_1 )`

ν₁ = 26.545

ν₂ = 33.697

Angle of convex corner with horizontal based on Prandtl-Meyer function = ν₂ – ν₁ = 7.152 degree

The geometrical angle of the convex corner with horizontal = 5.352 degree

Subsonic Flow-through Duct

As discussed earlier, we considered one more case for having a subsonic inlet flow-through duct. The contours for this don’t show anything similar to as seen for the supersonic case. Rather, it can be much helpful for locating the point emanating the shockwave in supersonic flow. In subsonic flow, it appears to be in the initial phase but not developed while in supersonic flow we have clearly seen the developed shockwave. Since all the contours for different parameters doesn’t show much except for locating convex corner point, only one contour has been shown here.

Animation

https://youtu.be/qAOHN3L1YKk

Conclusions

• Prandtl-Meyer Shock shows the opposite effect of an oblique shockwave.
  
  
• Adaptive Mesh Refinement is an important tool to capture shockwave effectively using a sub-grid criterion based on the curvature of a parameter such as a velocity or a temperature.
  
  
• Sub-grid criterion value defined must low enough to capture converge and accurate solution.
  
  
• Due to Prandtl-Meyer shock, flow undergoes an abrupt change in reduced flow properties such as density, temperature, and pressure (static) across the expansion fan.
  
  
• Simulating subsonic flow could help to locate the shockwave emanating point in the duct.