Week 5: Prandtl Meyer Shock problem

Introduction:

When a fluid flowing at supersonic velocity passes over a convex corner of the surface, it produces a discontinuity in flow

properties which is called the expansion fan. The fluid is flowing at so high speed that it tries to not follow the surface curve

and break off. It is not possible for the flow to separate directly so it expands suddenly and this causes shockwaves of

different Mach numbers spreading out like a fan. Expansion fans are isentropic processes that generate continuous and

smooth changes in the flow, causing the flow’s total properties to be conserved. The shock wave created is a type of

propagating disturbance that moves faster than the local speed of sound in the medium. The shockwave due supersonic

vehicle is created because the pressure variation created at a local point is more strong and faster than its neighbouring

points.

The use of this problem in the CFD domain is for validation purposes. The results from inviscid compressible fluid theory are

used for comparison of results from our CFD code.

For these kinds of supersonic flows, the pressure values are not defined at boundaries. We rather provide zero gradients or

Neumann boundary condition at the inlet boundary. The pressure values are extrapolated to boundaries from internal

regions. The region is initialized with proper velocity and pressure to avoid solution instability when the simulation starts. For

this case, the region is initialized with air at 50000 Pa and 680 m/s. So, when the simulation would start, the solver would

have a better initial guess. Also, at the outlet, the pressure condition is defined as a zero gradient boundary.

 

Q1. Shock Flow Boundary Conditions:

For solving the steady-state flow appropriate boundary conditions are needed. It is one of the required components of the

mathematical model. On the other hand, for solving transient flow, the appropriate initial condition is also required. 

Types of Boundary Conditions

Both ordinary and partial differential equations need boundary conditions to be solved. Different types of boundary conditions

can be imposed on the boundary of the domain. The choice of the boundary condition is very important as a bad imposition

of boundary condition may lead to the divergence of the solution or may also lead to the convergence of an incorrect

solution.

a) Dirichlet Boundary Condition:

In computational fluid mechanics, the classical Dirichlet boundary condition consists of the value of velocity and/or pressure

to be taken by a certain set of nodes.

i) Slip boundary condition: the velocity normal to the boundary is set to zero, while the velocity parallel to the boundary is let

free

ii) No-slip boundary condition: both the velocity normal to the boundary and the velocity parallel to the boundary are set

equal to zero.

At least one homogeneous boundary condition of the pressure/velocity has to be imposed as a reference for open domains.

b) Neumann Boundary Condition

When imposed on an ordinary or a partial differential equation, the Neumann boundary condition specifies the values that the

derivative of a solution is going to take on the boundary of the domain.

Constraints on the derivative of velocity can be seen in the application of a symmetry plane $\frac{\partial u}{\partial n}=0$

 

Since this condition is always applied in addition to a Dirichlet boundary condition on the velocity normal to the boundary, it is

naturally satisfied.

c) Robin Boundary Condition

The Robin boundary condition consists of a linear combination of the values of the field and its derivatives on the boundary.

Thus, it can also be said to be the linear combination of the Dirichlet and Neumann boundary conditions.

d) Mixed Boundary Condition

It consists of applying different types of boundary conditions in different parts of the domain. The mixed boundary condition

differs from the Robin condition because the latter consists of different types of boundary conditions applied to the same

region of the boundary, while the mixed condition implies different types of boundary conditions applied to different parts of

the boundary.

e) Cauchy Boundary Condition

The Cauchy boundary condition is a condition on both the unknown field and its derivatives. It differs from the Robin

condition because the Cauchy condition implies the imposition of two constraints (1 Dirichlet + 1 Neumann), while the Robin

condition implies only one constraint on the linear combination of the unknown function and its derivatives.

Why is Neumann Boundary condition used for supersonic Outlets?

Because the values of all the variables (pressure, velocity, temperature, Mach no) are not available before solving the

problem. So nothing is specified at a supersonic outlet and all variables are extrapolated from the domain interior. Only the

Normal derivative (Normal for the outlet boundary) of variables are specified as 0 (or Zero Gradient).

However, if any pressure value or velocity is specified at the outlet, the solution will not converge and become unstable.

 

Q2. What is a shock wave?

In physics, a shock wave is a type of propagating disturbance that moves faster than the local speed of sound in the medium.

Like an ordinary wave, a shock wave carries energy and can propagate through a medium but is characterized by an abrupt,

nearly discontinuous, change in pressure, temperature, and density of the medium.

When the speed of the moving object or source exceeds the speed of sound in the medium then the wavefronts lag behind

the source forming a cone-shaped region with a source at the vertex. The edge of the cone forms a supersonic wavefront

with an unusually large amplitude called a shock wave. A sonic boom is heard when the shock waves reach an observer. The

occurrence of shock waves can be characterized by the instantaneous change in pressure, velocity and temperature in a fluid

flow. The region between the vehicle and the shock wave known as the shock layer will be a region of high pressure, density

and temperature than the free-stream flow conditions. When a fluid streamline crosses the standing shock wave, an abrupt

increase in the pressure, temperature and density of the fluid flow occurs with a decrease in velocity of the flow.

 

Prandtl Meyer Expansion Fan:

A supersonic expansion fan, technically known as Prandtl–Meyer expansion fan, a two-dimensional simple wave, is a centred

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 backwards to meet at a point.

Each wave in the expansion fan turns the flow gradually. It is physically impossible for the flow to turn through a single

"shock" wave because this would violate the second law of thermodynamics. Across the expansion fan, the flow accelerates

(velocity increases) and the Mach number increases, while the static pressure, temperature and density decrease. Since the

process is isentropic, the stagnation properties (e.g. the total pressure and total temperature) remain constant across the

fan. 

To understand the Prandtl Meyer shock wave, there is a need to understand what oblique waves are. The normal shock waves

are straight in which the flow before and after the wave is normal to the shock. It is considered as a special case in the

general family of oblique shock waves that occur in supersonic flow. In general, oblique shock waves are straight but inclined

at an angle to the upstream flow and produce a change in the flow direction. An oblique shock generally occurs, when a

supersonic flow is ‘turned into itself”.

Another class of two-dimensional waves occurring in supersonic flow shows the opposite effects of oblique shock. Such types

of waves are known as expansion waves. When the supersonic flow is “turned away from itself”, an expansion wave is

formed. Here, the flow is allowed to pass over a surface which is inclined at an angle θ to the horizontal and all the flow

streamlines are deflected downwards. The change in flow direction takes place across an expansion fan centred at point ‘A'.

The flow streamlines are smoothly curved till the downstream flow becomes parallel to the wall surface behind the point ‘A'.

Here, the flow properties change smoothly through the expansion fan except at point ‘A'. An infinitely strong oblique

expansion wave may be called as a Mach wave. An expansion wave emanating from a sharp convex corner is known as a

centred expansion which is commonly known as a Prandtl-Meyer expansion wave.

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 and viscosity.

$\text{Mach Number at Particular Temperature}=\frac{\text{Velocity of Air}}{\text{Speed of sound}}$

 

Classification of flow based on Mach number –

Subsonic: M < 0.8

Transonic: M belongs to 0.8 – 1.3

Supersonic: M belongs to 1.3 – 5.0

Hypersonic: M belongs to 5.0 – 10.0

High-Hypersonic: M belongs 10.0 – 25.0

Re-entry Speeds: M > 25.0

 

Geometry Setup:

1. Import Geometry:

Imported the .stl file geometry.

 

Enabling geometry bounding box.

Geometry>Options>Checking Geometry bounding box

 

2. Scaling down the geometry

Converge assumes the dimensions in metres i.e in below-shown image l_x=65000m, l_y=45152.5, l_z=8151m which is too

large to handle. So, need to scale down.

Geometry>Transform>Scale>Entity Type-Entire Surface>Scale Factor-Uniform=0.001>Apply

 

 Desired Scaled-down Geometry:

 

3. Creating Boundaries:

Geometry>Boundary>Flag> '+' create a new boundary>Create Multiple boundaries>OK

 

4. Assigning Boundaries:

In converge every part of geometry is assumed to be a triangle. Selecting triangles and grouped to assign to specific

boundaries.

 

i) Inlet:

$△$ as an entity selection-By Angle ${20}^0$

> Using cursor pick option> Selecting the triangles as shown in below image> Assigning them to inlet Boundary> Apply

 

ii) Front 2D:

$△$ as an entity selection-By Angle ${20}^0$

> Using cursor pick option> Selecting the triangles as shown in below image> Assigning them to Front 2D Boundary> Apply

 

iii) Back 2D:

$△$ as an entity selection-By Angle ${20}^0$

> Using cursor pick option> Selecting the triangles as shown in below image> Assigning them to Back 2D Boundary> Apply

 

iv) Top and Bottom Wall:

$△$ as an entity selection-By Angle ${20}^0$

> Using cursor pick option> Selecting the triangles as shown in below image> Assigning them to Top and Bottom Wall

Boundary> Apply

 

v) Outlet:

$△$ as an entity selection-By Angle ${20}^0$

> Using cursor pick option> Selecting the triangles as shown in below image> Assigning them to Front 2D Boundary> Apply

 

The geometry with assigned boundaries: 

 

4. The orientation of Normals.

Every geometry has a normal vector, it points perpendicular to it.

Normal Toggle - Enables Normal

To run the CFD simulation, the normal should point towards the fluid.

 

5. Diagnosis

This option will reflect the errors associated with geometry if there are any.

In this case, geometry is not having any errors.

 

Case Setup:

1. Application Type

 Selecting an application type as "Time Based".

 

2. Materials:

Selecting Air in the drop-down menu of predefined mixtures and keeping Gas Simulation enabled.

Confirming the permission to overwrite the gas.dat, mech.dat and therm.dat by loading predefined mixture Air.

 

a) Gas Simulation:

Keeping the default values for gas simulation > OK

 

b) Global Transport Properties:

Keeping all the default values for global transport properties > OK

 

3. Simulation Parameters:

 

a) Run Parameters:

Under Solver tab change the solver to Steady-state solver.

Simulation mode is Full hydrodynamic and Gas flow solver is Compressible.

Under Misc. tab checking the Steady-state monitor option and disabling use shared memory option  > OK.

 

b) Simulation Time Parameters:

Set the start time as 0, end time as 25000 cycles, Initial and minimum time step as 1e-9 and Maximum time step as 1

second.

Keep all the other CFL values default > OK.

Using initial and maximum time-step as 1e-9 because converge uses a dynamic time-step algorithm. At every time-step

converge tries to see if the time-step can be increased for next coming time-step.

So, if initial conditions and mesh conditions are good then time is going to increase automatically. 

 

c) Solver Parameters:

Changing the Navier-stokes solver type as Density-based.

Keeping all other settings as default > OK.

 

4. Regions and Initializations:

Grouping the boundaries to form volumetric zones called region. This defines a volumetric region with initial conditions like

initial pressure, initial temperature, velocity and species concentration 

Selecting Add to create a new region, then selecting Air for Species initial conditions and keep the other values default > OK.

Initial Conditions are:

Velocity = 680 m/s in x direction,

Temperature = 300 K,

Pressure = 50000 Pa and Species as Air

 

5. Boundary:

i) Inlet:

Boundary Type - INFLOW

Pressure Boundary Condition - Zero normal gradient (NE)

The velocity boundary condition - Specified value of 680 m/s in the x-direction.

Temperature boundary condition - Specified value of 286.1 K.

Air is selected for species boundary condition.

 

ii) Front 2D:

Boundary Type - TWO D

 

iii) Back 2D:

Boundary Type - TWO D

 

iv) Top & Bottom Wall:

Boundary type - WALL

Velocity boundary condition - Slip

Temperature boundary condition - Zero normal gradient

 

v) Outlet:

Boundary Type - OUTFLOW

Pressure boundary conditions - Zero Normal gradient.

Velocity boundary conditions - Zero Normal gradient

For Species Backflow Air is selected.

 

6. Turbulence Modeling:

Keep the Turbulence model as RNG k-epsilon

Keep all the settings default > OK.

 

6. Grid Control

Enabling Adaptive Mesh Refinement option.

 

i) Base Grid:

Base Grid size = 0.8 m

 

ii) Adaptive Mesh Refinement:

Adaptive Mesh Refinement is a method of adapting the accuracy of a solution within certain sensitive or turbulent regions of

simulation, dynamically and during the time the solution is being calculated. When solutions are calculated numerically, they

are often limited to pre-determined quantified grids as in the Cartesian plane which constitutes the computational grid, or

'mesh'.

The available region is Region 0 which is changed to active regions.

Minimum cells - 1

Maximum cells - 200000

Embedded type - Sub-Grid-Scale (SGS)

Temperature scheme is used for AMR with Embed type -SGS (Sub Grid-Scale)

Sub-grid criterion - 0.05 K 

Max embedding level as 2

Timing-control type - Sequential

The sub-grid criterion is varied for different cases.

How the Refinement of grid size is done using Embedded sizing is:

The data structure the Converge Software uses is the Octri-data structure.

Refinement Formula = Grid Size = $\left(\frac{N_{basesize}}{2^{Embed-Level}}\right)$

                                                = $\frac{0.8}{2^2}$

                                                = $\frac{0.8}{4}$ 

                                               = 0.2 m.

 

7. Output/Post-Processing

i) Post-Variable Selection

Using default settings

 

 

ii) Output Files

Setting the time interval for writing 3D output data files and interval for writing restarting output to 100 and keep all the

other values default > OK

 

8. Export all the input files

Export all the input files into a separate folder

Files>Export>Export Input Files>setting up the desired location>OK

 

9. Copy and paste the mpiexec.exe file to the input files folder 

 

 

Run Simulation

Cygwin - Cygwin is a collection of GNU and open-source tools that provide functionality similar to Linux distribution on

windows.

To run the simulation, open Cygwin and navigate to the folder in which the input files are exported.

 

To run the simulation command with executable to be entered

$ mpiexec.exe -n 2 converge-intelmpi.exe restricted logfile.txt &

Taking output mpiexec.exe and sending it to  further Taking data from   and sending to logfile.txt

 

Results:

Base Grid:

 

Case 1:

Velocity = 680 m/s

SGS = 0.05 K

Computational Time:

 

Meshing:

SGS=0.05 K

 

Animation: 

From the mesh profile, the advantage of adaptive mesh refinement parameter is visible. It has embedded the mesh size by 2

levels with the SGS of 0.05 K. This refinement helps in capturing characteristics of the flow through the domain specifically at

the sharp corner, to have a better analysis of the shock wave. 

The mesh refinement with SGS of 0.05 K, it fails to refine the complete expansion fan region, which will affect the analysis of

the flow properties in those unrefined regions. 

 

Pressure Contour:

 

 Total Pressure:

 

Static Pressure:

 

 

Velocity Contour:

 

Velocity Plot:

 

 

Temperature Contour:

 

Temperature Plot: 

 

Density Contour:

 

Density Plot:

Mach No.:

Above plots and contours representing the variation of pressure, velocity, temperature, density, and Mach Number (Flow

velocity/speed of sound) through the domain length. As the flow propagates through the domain at a speed higher than the

speed of sound and with an abrupt increase of domain area, centred expansion waves or centred expansion fan is generated.

Hence, with the generation of the Prandtl-Meyer expansion fan, the Mach number and velocity increases whereas the flow

properties such as the density, temperature, and pressure decreases.

Mass Flow Rate:

Cell Count:

 

Case 2:

Velocity = 680 m/s

SGS = 0.04 K

Computational Time:

 

 

Meshing:

SGS=0.04 K

 

From the mesh profile, the advantage of adaptive mesh refinement parameter is visible. It has embedded the mesh size by 2

levels with the SGS of 0.04 K. This refinement helps in capturing characteristics of the flow through the domain specifically at

the sharp corner, to have a better analysis of the shock wave.

When compared to the mesh refinement with SGS of 0.05 K, it refines the expansion fan region better, for the analysis of the

flow properties in those regions. 

 

Pressure Contour:

 

 

Pressure Plot:

Total Pressure:

 

Static Pressure:

 

 

Velocity Contour:

 

Velocity Plot:

 

 

Temperature Contour:

  

Temperature Plot: 

 

Density Contour:

 

Density Plot:

Mach No.:

Above plots and contours representing the variation of pressure, velocity, temperature, density, and Mach Number (Flow

velocity/speed of sound) through the domain length. As the flow propagates through the domain at a speed higher than the

speed of sound and with an abrupt increase of domain area, centred expansion waves or centred expansion fan is generated.

Hence, with the generation of the Prandtl-Meyer expansion fan, the Mach number and velocity increases whereas the flow

properties such as the density, temperature, and pressure decreases.

 

Mass Flow Rate:

Cell Count:

 

 

Case 3:

Velocity = 680 m/s

SGS = 0.03 K

Computational Time: 

 

 

Meshing:

SGS=0.03 K

  

From the mesh profile, the advantage of adaptive mesh refinement parameter is visible. It has embedded the mesh size by 2

levels with the SGS of 0.03 K. This refinement helps in capturing characteristics of the flow through the domain specifically at

the sharp corner, to have a better analysis of the shock wave.

When compared to the mesh refinement with SGS of 0.04 K, it refines the expansion fan region near to complete refinement,

for the analysis of the flow properties in those regions. 

 

Pressure Contour:

  

Pressure Plot:

Total Pressure:

 

Static Pressure:

 

 

Velocity Contour:

 

Velocity Animation:

 

 

Velocity Plot:

 

 

Temperature Contour:

 

Temperature Animation:

 

 

Temperature Plot: 

 

Density Contour:

 

Density Plot:

Mach No.:

Above plots and contours representing the variation of Mach Number (Flow velocity/speed of sound), density, temperature

and pressure through the domain length. As the flow propagates through the domain at a speed higher than the speed of

sound and with an abrupt increase of domain area, centred expansion waves or centred expansion fan is generated.

Hence, with the generation of the Prandtl-Meyer expansion fan, the Mach number increases whereas the flow properties such

as the density, temperature, and pressure decreases.

 

Mass Flow Rate:

Cell Count:

Case 4:

Velocity = 680 m/s

SGS = 0.01 K

Computational Time: 

 

Meshing:

SGS=0.01 K

 

Animation:

 

From the mesh profile, the advantage of adaptive mesh refinement parameter is visible. It has embedded the mesh size by 2

levels with the SGS of 0.03 K. This refinement helps in capturing characteristics of the flow through the domain specifically at

the sharp corner, to have a better analysis of the shock wave.

When compared to the mesh refinement with SGS of all the above cases, it gives the expansion fan region complete

refinement, for the analysis of the flow properties in those regions. 

 

Pressure Contour:

  

 Pressure Plot:

Total Pressure:

 

Static Pressure:

 

 

Velocity Contour:

  

 

Velocity Plot:

 

 

Temperature Contour:

 

Temperature Plot: 

 

Density Contour:

 

Density Plot:

Mach No.:

Above plots and contours representing the variation of Mach Number (Flow velocity/speed of sound), density, temperature

and pressure through the domain length. As the flow propagates through the domain at a speed higher than the speed of

sound and with an abrupt increase of domain area, centred expansion waves or centred expansion fan is generated.

Hence, with the generation of the Prandtl-Meyer expansion fan, the Mach number increases whereas the flow properties such

as the density, temperature, and pressure decreases.

We can see the more refined solution and better looking contour towards the end of the fan and along its edges, it has

happened due to a stricter SGS criterion. It is understood that as the SGS will be reduced we can expect a more refined

solution with an increase in mesh sizes and better result. 

 

Mass Flow Rate:

Cell Count:

 

Case 5:

Velocity - 100 m/s

SGS - 0.05 K 

 

 

Meshing:

As the flow is subsonic, there will be no much change in the temperature. So, the Adaptive mesh refinement has not

appeared on the surface.

The base-grid size is 0.8m.

In the case of subsonic flow no matter how small the SGC value is we cannot observe any adaptive mesh refinement. This

indicates that the curvature of temperature is not changing as the flow progresses.

 

Pressure Contour:

 

 

Pressure Plot:

Total Pressure:

 

 

Static Pressure:

 

 

Velocity Contour:

The Inlet velocity given is 100m/s. Due to the viscous force and resistances, the maximum velocity recorded at the end of

the simulation is 108.4 m/s.

The min. velocity recorded inside is 91 m/s.  Because this is a subsonic problem, the shock wave is not produced. So that the

AMR does not come into the picture for this simulation.

 

 

Velocity Plot:

 

 

Temperature Contour:

The above contour shows how the temperature is behaving inside the duct when flow happens inside.

 

 

Temperature Plot: 

The average inlet temperature is 286.1 K given from the inlet boundary condition.

The temperature at the outlet is gradually decreasing and settles to around 279 K due to the decrease of the pressure inside

the duct.

 

Density Contour:

 

Density Plot:

 

Mass Flow Rate:

Cell Count:

A total of 4300 cells has been generated, and due to the subsonic flow phenomenon, there is no AMR. So, there is no

disturbance in the curve of the cell-count.

 

Subsonic Vs Supersonic Flow Comparison

In supersonic flow, the Prandtl-Meyer expansion fan is visible and the mesh refinement is also visible in the shock region

Whereas on the other hand, in case of a subsonic flow, the Prandtl-Meyer expansion fan is not at all visible and even the

mesh refinement is not visible. This is because there are no sudden fluctuations in pressure, velocity, temperature or density.

Also, the average Mach number increases from the inlet to the outlet in case of supersonic flow and decreases in case of

subsonic flow.

 

 Q3. Effect of SGS parameter on shock location and cell count.

SGS is a sub-grid scale, which is a mesh refiner in the adaptive mesh refinement that depending on the gradients among the

parameters. SGS refines the mesh grid by the variations in the curve of the parameters. This help in capturing the details of

the flow in the gradient region. 

Here gradient refers to the second-order PDE of the parameter $\frac{{\partial }^{2}T}{\partial x^{}}$

If we specify the value like temperature to be 0.01 degree in subgrid criteria(SGC), then if the rise in temperature in the

region is above the average temperature of the region the mesh gets finer depending on the scale is set. This adaptive

technique provides detailed insight into the turbulent regions.