Prandtl Meyer Shock problem
Prandtl-Meyer Shock Problem
Overview
The current report is divided into three sections. In the section I, a brief theoretical background of shocks and Prandtl-Meyer expansion fan is provided. Section II begins with an introduction to boundary conditions which are used in Mathematics to solve differential equations. Additionally section II also discusses some typical boundary conditions used in the field of Computational Fluid Dynamics (CFD). Section II ends with a short discussion on the inlet and outlet boundary conditions used for the case of steady inviscid supersonic flow. Section III consists of details regarding the simulation of Prandtl-Meyer expansion fan: geometry, materials, solver, initial and boundary conditions, grid, results and explanations.
I. Physics of shocks and expansion fans
A shock wave or shock is a disturbance travelling through a medium at a speed greater than the local speed of sound in that medium. Shocks are accompanied by abrupt changes in pressure, temperature and density of the medium. These changes occur across such small thickness that shocks appear as discontinuities.
In order to understand shocks, let us consider a scenario where an object (source of disturbance) is moving through a quiescent fluid at the local speed of sound in the fluid. The figure shown below displays the spherical wavefronts emanating from the object.
It can be seen that the object keeps up with the spherical wavefronts that it creates since the radius of each spherical wavefront is given by `aDeltat` where `Deltat` is the time elapsed and the distance moved by the source within the same duration is also `aDeltat`.
Let us now consider the scenario where the speed of the object exceeds the local speed of sound in the medium. In other words, the object is moving at supersonic speed.
As shown in the figure, at supersonic speed the spherical wavefronts are unable to keep up with the object (source of disturbance) since within the same time duration `Deltat` the wavefronts travel a distance of `aDeltat` whereas the object moves through a (greater) distance of `vDeltat`. At any time instant the different wavefronts are encompassed by a cone in three dimensions, called the Mach cone. Due to the grouping (constructive interference) of pressure wavefronts at the surface of the cone, shock waves are established. Any region of the fluid which passes through the mach cone, experiences shock and undergoes abrupt changes in pressure, temperature and density.
The half angle `theta` of the mach cone is known as Mach cone angle.
`sin theta=(aDeltat)/(vDeltat)=a/v=1/M` (Mach number)
`theta=arcsin (1/M)`
Like other waves, a shock wave transmits energy through the medium. Across the shock waves, energy is conserved but entropy is created. Hence shocks are highly irreversible processes. The second law of thermodynamics prevents the existence of shocks across which the fluid undergoes expansion. However there are scenarios in which supersonic flows can experience isentropic expansions and compressions. One such scenario where the supersonic flow experiences isentropic expansion at a sharp corner is Prandtl-Meyer expansion fan. Given below is an image showing the Prandtl-Meyer expansion fan.
image taken from https://en.wikipedia.org/wiki/Prandtl%E2%80%93Meyer_expansion_fan
The isentropic expansion experienced by the supersonic flow in Prandtl-Meyer expansion occurs across infinite number of Mach waves (all diverging from the sharp corner), each wave turning the flow by an infinitesimal angle.
Prandtl-Meyer function determines the angle through which a sonic flow (M=1) must turn to reach a particular Mach number (M). Mathematically,
`nu(M)=sqrt((gamma+1)/(gamma-1))arctan (sqrt((gamma+1)/(gamma-1)(M^2-1)))-arctan (sqrt(M^2-1))`
Also the turn angle `theta` is given by:
`theta=nu(M_2)-nu(M_1)`
II. Boundary conditions
Boundary conditions are essential for obtaining unique solution of differential equations (ODEs or PDEs). They are what separate the flow over an automobile from the flow over a supersonic jet (say), otherwise for both scenarios the governing PDEs are the same.
In this section, various types of boundary conditions used in Mathematics for solving differential equations are discussed. Additionally, the typical boundary conditions applicable to computational fluid dynamics problems are listed. A special attention is given to supersonic flow problems.
1. Boundary conditions used for solving differential equations in mathematics
The different types of boundary conditions that are used for solving differential equations in mathematics can be grouped as follows:
i. Dirichlet boundary condition: Specifies the values of a solution at the boundary at which it is imposed.
ii. Neumann boundary condition: Specifies the values of the normal derivative of a solution at the boundary at which it is imposed.
iii. Cauchy boundary condition: Specifies the values of both, a solution and its normal derivative at the boundary at which it is imposed.
iv. Mixed boundary condition: At the boundary at which it is imposed, it specifies the values of a solution over one part and the values of its normal derivative over the other part.
v. Robin boundary condition: Specifies the linear combination of the values of a solution and the values ofits normal derivative at the boundary at which it is imposed.
2. Typical Boundary conditions applicable to Computational Fluid Dynamics problems
Given below are the typical boundary conditions applied to CFD problems:
i. Wall: The wall type boundary condition can be imposed in the following ways:
(a) Solid wall: A solid wall can have slip or no-slip boundary condition depending on whether the flow is inviscid or viscous respectively.
- Slip: v=0 (No penetration boundary condition) : applicable to inviscid flows
- No-slip: u=v=w=0 : applicable to viscous flows
(b) Moving wall: For a moving wall the fluid in contact attains the same velocity as the wall.
`u=u_W; v=v_W; w=w_W`
(c) Porous wall: For a porous wall generally the pressure gradient through it is provided.
ii. Far-field: The far-field type boundary condition can be imposed in the following ways:
(a) Inlet/Inflow: Generally freestream conditions are applied under Inlet/Inflow type boundary condition.
`u=U_(oo); p=P_(oo)`
(b) Outlet/Outflow: The following boundary conditions are generally applied under Outlet/Outflow type boundary condition:
- Freestream
- Prescribed pressure
- Pressure extrapolated from interior (Neumann)
iii. Symmetry: The symmetry type boundary condition implies the following:
- No flow across the boundary
- No scalar flux across the boundary
iv. Periodic: The periodic type boundary condition implies that the flux of all variables at the periodic boundaries are repeated.
3. Case of steady invscid supersonic flow
Any flow can be determined by the so-called characteristic curves or wavefronts, which are the curves along which the flow gradients are indeterminate. They determine the direction of propagation of information through the fluid.
For hyperbolic PDEs, the number of boundary conditions that must be supplied at a boundary is equal to the number of characteristics that are directed from the exterior of the domain towards that boundary. The case of steady inviscid supersonic flow is governed by hyperbolic PDEs. Attached below are the pictures showing characteristics at the supersonic-inlet and supersonic-outlet boundaries for a one dimensional supersonic flow.
Supersonic-inlet
All three characteristics are directed from the exterior of the domain towards the supersonic-inlet boundary and no information is transmitted from the domain to the boundary. In other words, all three characteristics must be defined at the supersonic-inlet boundary.
Supersonic-outlet
None of the characteristics are directed from the exterior of the domain towards the supersonic-outlet boundary and no information is transmitted from the boundary to the domain. In other words, none of the three characteristics need to be defined at the supersonic-outlet boundary. Information from the solution should be used to update the boundary values. This makes Neumann type boundary condition, described in section II.1, an obvious choice for the supersonic-outlet boundary. For this reason while simulating flows with supersonic-outets, the zero gradient boundary conditions are set for pressure and velocity.
III. Simulation of Prandtl-Meyer expansion fan
In the current section, the Prandtl-Meyer expansion fan is simulated using CONVERGE CFD for the following cases:
Case 1: Initial velocity inside the domain=640m/s, SGS parameter=0.05K
Case 2: Initial velocity inside the domain=640m/s, SGS parameter=0.1K
Case 3: Initial velocity inside the domain=640m/s, SGS parameter=0.2K
Case 4: Initial velocity inside the domain=100m/s, SGS parameter=0.05K
The aim of the simulation is to study the following:
(a) Effect of sub-grid scale (SGS) parameter on the expansion fan and cell count.
(b) Effect of initial conditions (specifically velocity) on the results.
1. Geometry
The figure given below shows the geometry of the problem.
2. Material
Air (Predefined mixture) is used as the working fluid for performing gas simulation. Accordingy, the following global properties are utilized:
Turbulent Prandtl number: 0.9
Turbulent Schmidt number: 0.78
3. Solver
The Navier Stokes equations are solved using CONVERGE solver with the following settings:
i. Solver: Steady-state solver, Density-based
ii. Navier-Stokes solver scheme: PISO
iii. Simulation mode: Full hydrodynamic
iv. Gas flow solver: Compressible
v. Initial time-step: 1e-7s
vi. Minimum time-step: 1e-7s
vii. Maximum time-step: 1s
For performing steady-state simulation, the steady-state monitor is set up using average velocity at the outlet.
4. Initial conditions
The following settings are applied as the initial conditions:
Velocity: X-component: 680m/s; Y-component: 0; Z-component: 0 (For Cases 1, 2 and 3)
Velocity: X-component: 100m/s; Y-component: 0; Z-component: 0 (For Case 4)
Temperature: 286.1K
Pressure: 101325Pa
Species: Mass fractions: Oxygen: 0.23; Nitrogen: 0.77
5. Boundary conditions
The following boundary conditions are applied:
i. Inlet: INFLOW type boundary condition is applied with the following settings:
Pressure: 101325Pa
Velocity: X-component: 680m/s; Y-component: 0; Z-component: 0
Temperature: 286.1K
Species: Mass fractions: Oxygen: 0.23; Nitrogen: 0.77
ii. Outlet: OUTFLOW type boundary condition is applied with the following settings:
Pressure: Zero normal gradient
Velocity: Zero normal gradient
Temperature backflow: Zero normal gradient
Species backflow: Mass fractions: Zero normal gradient
iii. Slip_Walls: WALL type boundary condition is applied at the top and bottom walls with the following settings:
Wall motion type: Stationary
Surface movement: FIXED
Velocity: Slip
Temperature: Zero normal gradient
iv Front_Wall: TWO_D type boundary condition is imposed on the front wall. This ensures that the problem is solved in two dimensions.
v. Back_Wall: TWO_D type boundary condition is imposed on the back wall. This ensures that the problem is solved in two dimensions.
6. Grid
Computational grid with the following features is used for all cases:
i. Base grid: dx=dy=dz=0.8m
ii. Adaptive Mesh Refinement (AMR):
Minimum cells: 1
Maximum cells: 200000
Temperature Sub-Grid Scale (SGS) criterion: 0.05K (Case 1 and Case 4); 0.1K (Case 2); 0.2K (Case 3)
Maximum embedding level: 2
Timing control type: SEQUENTIAL
Start: 5000 cycles
End: 999999999 cycles
The figure shown below shows the base grid for all cases.
At steady state, the computational grids for the different cases are given below.
Case 1: Initial velocity inside the domain=640m/s, SGS parameter=0.05K
Case 2: Initial velocity inside the domain=640m/s, SGS parameter=0.1K
Case 3: Initial velocity inside the domain=640m/s, SGS parameter=0.2K
Case 4: Initial velocity inside the domain=100m/s, SGS parameter=0.05K
Discussion
It is evident that with an increase in the SGS criterion for temperature, the number of embedded cells decrease. This is because with an increase in the SGS criterion, the condition required for embedding new cells gets relaxed. Also the computational grids for Case 1 and Case 4 are same, since in both cases an SGS criterion of 0.05K is applied. At this point it can be easily predicted that the computational grid for Case 1 (and Case 4) is best suited to capture the expansion fan as it can resolve the smallest temperature gradients among all grids.
7. Results
7.1. Temperature contours
Section 7.1.1 shows the steady state temperature contours whereas section 7.1.2 shows the evolution of temperature contours with simulation cycles using animations.
7.1.1. Steady state
Case 1: Initial velocity inside the domain=640m/s, SGS parameter=0.05K
Case 2: Initial velocity inside the domain=640m/s, SGS parameter=0.1K
Case 3: Initial velocity inside the domain=640m/s, SGS parameter=0.2K
Case 4: Initial velocity inside the domain=100m/s, SGS parameter=0.05K
Discussion
The plots show that across the expansion fan the air temperature is reduced. For all four cases, the expansion fan is located near the sharp corner but the tip of the fan is best resolved for Case 1 (and Case 4), followed by Case 2 and Case 3 respectively. In other words, with an increase in the SGS criterion the resolution of the tip of the expansion fan decreases. Also the temperature contours for Case 4 are almost the same as those for Case 1.
7.1.2. Animations
Case 1: Initial velocity inside the domain=640m/s, SGS parameter=0.05K
Case 2: Initial velocity inside the domain=640m/s, SGS parameter=0.1K
Case 3: Initial velocity inside the domain=640m/s, SGS parameter=0.2K
Case 4: Initial velocity inside the domain=100m/s, SGS parameter=0.05K
Discussion
The animations show that the expansion fan is better visualized for each case after 5s, when the AMR feature is enabled.
7.2. Pressure contours
Section 7.2.1 shows the steady state pressure contours whereas section 7.2.2 shows the evolution of pressure contours with simulation cycles using animations.
7.2.1. Steady state
Case 1: Initial velocity inside the domain=640m/s, SGS parameter=0.05K
Case 2: Initial velocity inside the domain=640m/s, SGS parameter=0.1K
Case 3: Initial velocity inside the domain=640m/s, SGS parameter=0.2K
Case 4: Initial velocity inside the domain=100m/s, SGS parameter=0.05K
Discussion
The plots show that across the expansion fan the air pressure is reduced. With an increase in the SGS criterion, the resolution of expansion fan decreases and it becomes more diffused. Also the pressure contours for Case 4 are almost the same as those for Case 1.
7.2.2. Animations
Case 1: Initial velocity inside the domain=640m/s, SGS parameter=0.05K
Case 2: Initial velocity inside the domain=640m/s, SGS parameter=0.1K
Case 3: Initial velocity inside the domain=640m/s, SGS parameter=0.2K
Case 4: Initial velocity inside the domain=100m/s, SGS parameter=0.05K
Discussion
The animations show that the expansion fan is better visualized for each case after 5s, when the AMR feature is enabled.
7.3. Velocity contours
Section 7.3.1 shows the steady state velocity contours whereas section 7.3.2 shows the evolution of velocity contours with simulation cycles using animations.
7.3.1. Steady state
Case 1: Initial velocity inside the domain=640m/s, SGS parameter=0.05K
Case 2: Initial velocity inside the domain=640m/s, SGS parameter=0.1K
Case 3: Initial velocity inside the domain=640m/s, SGS parameter=0.2K
Case 4: Initial velocity inside the domain=100m/s, SGS parameter=0.05K
Discussion
The plots show that across the expansion fan the flow is accelerated. With an increase in the SGS criterion, the resolution of expansion fan decreases and it becomes more diffused. Also the velocity contours for Case 4 are almost the same as those for Case 1.
7.3.2. Animations
Case 1: Initial velocity inside the domain=640m/s, SGS parameter=0.05K
Case 2: Initial velocity inside the domain=640m/s, SGS parameter=0.1K
Case 3: Initial velocity inside the domain=640m/s, SGS parameter=0.2K
Case 4: Initial velocity inside the domain=100m/s, SGS parameter=0.05K
Discussion
The animations show that the expansion fan is better visualized for each case after 5s, when the AMR feature is enabled.
7.4. Density contours
Section 7.4.1 shows the steady state density contours whereas section 7.4.2 shows the evolution of density contours with simulation cycles using animations.
7.4.1. Steady state
Case 1: Initial velocity inside the domain=640m/s, SGS parameter=0.05K
Case 2: Initial velocity inside the domain=640m/s, SGS parameter=0.1K
Case 3: Initial velocity inside the domain=640m/s, SGS parameter=0.2K
Case 4: Initial velocity inside the domain=100m/s, SGS parameter=0.05K
Discussion
The plots show that across the expansion fan the density decreases. With an increase in the SGS criterion, the resolution of expansion fan decreases and it becomes more diffused. Also the density contours for Case 4 are almost the same as those for Case 1.
7.4.2. Animations
Case 1: Initial velocity inside the domain=640m/s, SGS parameter=0.05K
Case 2: Initial velocity inside the domain=640m/s, SGS parameter=0.1K
Case 3: Initial velocity inside the domain=640m/s, SGS parameter=0.2K
Case 4: Initial velocity inside the domain=100m/s, SGS parameter=0.05K
Discussion
The animations show that the expansion fan is better visualized for each case after 5s, when the AMR feature is enabled.
7.5. Cell count
The total cell count for all cases are plotted against the simulation cycles.
Discussion
Since AMR is enabled, the total cell count changes with simulation cycles. Also the highest cell count is generated for Case 1 (and Case 4), as expected, followed by Case 2 and Case 3 respectively. In other words, the total cell count decreases with an increase in the SGS criterion.
7.6. Mach number
The plots given below show the inlet and outlet Mach numbers at different simulation cycles, for all cases.
Case 1: Initial velocity inside the domain=640m/s, SGS parameter=0.05K
Case 2: Initial velocity inside the domain=640m/s, SGS parameter=0.1K
Case 3: Initial velocity inside the domain=640m/s, SGS parameter=0.2K
Case 4: Initial velocity inside the domain=100m/s, SGS parameter=0.05K
Discussion
The steady state Mach numbers at the outlet predicted by the simulations are as follows:
Case 1: 2.1338135
Case 2: 2.1342893
Case 3: 2.1343781
Case 4: 2.133992
The Mach number plots confirm that steady state is achieved for all cases. Additionally, the increase in Mach number at the exit suggests that air expands as it passes over the sharp corner.
Even though Case 4 starts with a different initial condition for velocity (100m/s), it achieves the same steady state as the other cases. This agrees with the idea explained earlier in section I.3. In case of supersonic flows the characteristics are directed from the inlet boundary to the domain and from the domain to the outlet boundary. Consequently the initial conditions affect the solution at the outlet boundary initially. However as the simulation progresses, the effect of inlet boundary conditions propagates through the domain and eventually reaches the outlet boundary. Hence if the initial conditions are not totally wrong or unphysical, their effect eventually dies down and the solution at steady state is ultimately determined by the inlet boundary conditions.
7.7. Comparison of Exit Mach number with theoretical value
As discussed previously, in section I.1, the Prandtl-Meyer function can be used to determine the turn angle for the flow across the expansion fan, provided that the inlet and outlet Mach numbers are known. For the present problem, the turn angle and the inlet Mach number are known. Based on these the Prandtl-Meyer function based on the outlet Mach number can be calculated, which can then be used to extract the analytical value of the outlet Mach number. This value can be compared with the values obtained from simulations.
`M_1=2.00074`
`nu(M_1)=sqrt((gamma+1)/(gamma-1))arctan (sqrt((gamma+1)/(gamma-1)(M_1^2-1)))-arctan (sqrt(M_1^2-1))=0.460769678148232`
Turn angle, `theta=0.093410021 rad`(based on the geometry)
Therefore:
`nu(M_2)=nu(M_1)+theta`
`sqrt((gamma+1)/(gamma-1))arctan (sqrt((gamma+1)/(gamma-1)(M_1^2-1)))-arctan (sqrt(M_1^2-1))=0.554179699148232`
Solving the above equation for `M_2`, we get:
`M_2=2.2008`
The table given below shows the comparison between the outlet Mach number obtained analytically and the outlet Mach numbers obtained for different simulations.
The table shows that the simulation results are within 5% deviation of the analytical results.
7.8. Conclusions
i. The AMR feature is better able to capture the expansion fan as compared to the base grid.
ii. Although the tip of the expansion fan is located near the sharp corner for all cases, its resolution increases with a decrease in the SGS criterion. For greater values of SGS criterion, the tip becomes diffused.
iii. The gradients - across the expansion fan - are better captured by simulations having small values of SGS criterion.
iv. With an increase in SGS criterion, the total cell count decreases.
v. In case of steady inviscid supersonic flows, the steady state solution is unaffected by the initial conditions provided that they are not unphysical.