Prandtl-Meyer Shock problem
Introduction
In this challenge we are setting up a simulation to capture expansion shock wave, when supersonic flow turns around a convex corner.
In our case we have a domain, in which we set inlet velocity as a 680m/s when the temperature is 286.1K, meaning that we have Mach = 2 at the inlet as a boundary condition.
The turning angle in our case is (180-174.648 = 5.352 degrees)
Supersonic Boundary Condition
At the outlet as a boundary condition we use Neumann Boundary condition, which is imposed on an ODE or PDE, the condition specifies the values in which the derivative of a solution is applied within the boundary of the domain.
We are using Neumann Boundary condition at the supersonic outlet because we have Velocity inlet BC, when the exit flow is fully developed, as an outflow Boundary condition assumes a zero normal gradient for all flow variables except pressure. The solver extrapolates the required information from interior.
https://en.wikipedia.org/wiki/Neumann_boundary_condition
http://www.bakker.org/dartmouth06/engs150/06-bound.pdf
Shock Wave
In physics Shock is a type of a propagating disturbance that moves faster than speed of sound. Shock wave carries energy and propagates through a medium but it is characterized as an unexpected, nearly discontinious change in pressure, density and temperature of the medium.
we know several types of shock waves.
1) Normal shock wave - at `90^circ` (perpendicular) to the shock medium\'s flow direction.
2) Oblique - at an angle to the direction of the flow
3) bow - occurs upstream of the front when the upstream flow velocity exceeds Mach 1.
Prandtl-Meyer Expansion Fan Theory
A supersonic expansion fan occurs when flow turns around the convex corner.
when flow turns at an angle it forms an expansion fan which consists of an infinite number of mach waves centered at the corner.
According to the second law of thermodynamics,which says that the total entropy of the isolated system can never decrease over time and it\'s constant if and only if all processes are reversible, it is impossible for the flow to turn through a single shock wave.
As a supersonic flow turns, normal component of the velocity increases `W_2>W_1`, while the tangential component ramains the same `v_1=v_2`. the corresponding change is Entropy `Deltas= s_2-s_1`
If `W_2>W_1`, `Deltas<0`. It\'s impossible, meaning it\'s impossible supersonic flow to turn thought only one shock wave. only one scenario is acceptable where we have turn through infinite shock waves, thus `Deltas->0`.
Across an expansion fan flow velocity increases, in other words Mach number increases, whereas static pressure temperature and density decreases.
First Mach line is at an angle `mu_1 = arcsin(1/M_1)` and the last mach line is at an angle `mu_2 = arcsin(1/M_2)`
`T_2/T_1 = ((1+(gamma-1)/2*M_1^2)/(1+(gamma-1)/2*M_2^2))`
`P_2/P_1 = ((1+(gamma-1)/2*M_1^2)/(1+(gamma-1)/2*M_2^2))^(gamma/(gamma-1))`
`rho_2/rho_1 = ((1+(gamma-1)/2*M_1^2)/(1+(gamma-1)/2*M_2^2))^(1/(gamma-1))`
Mach number after the turn `M_2`is related to initial Mach number `M_1` and turn angle `theta`.
`theta = nu(M_2)-nu(M_1)`
`nu(M)` is a Prandtl-Meyer function, which determines the turning angle through which given Mach number reaches a particular Mach number.
`nu(M) = sqrt((gamma+1)/(gamma-1))*arctan(sqrt((gamma-1)/(gamma+1)*(M^2-1)))-arctan(sqrt(M^2-1))`
In our case our inflow velocity is 680m/s, in other words Mach 2 (`M_1 = 2`) and the turning angle is 5.352degrees, in other words 0.09341rads.
code in matlab provided below
after substituting our values to the equations derived above we get
`nu(M_1) = 0.4604 => nu(M_2) = nu(M_1)+theta = 0.5538`
`M_2 = 2.19997`
`T_2 = 261.68K`
https://en.wikipedia.org/wiki/Prandtl–Meyer_expansion_fan
SET UP
we\'ve already talked about inlet and outlet boundary conditions.
Fluid: Air (Redlich-kwong real gas)
velocity baced BC at Inlet: 680m/s
Temperature at the Inlet : 286.1K
Outlet: Zero normal gradient (Neumann Boundary condition)
Top and Bottom Walls: slip condition
Front and Back Walls: 2D
Simulation time parameters.
Start time: 0 cyc
End time: 25000 cyc
Initial time-step: 1e-7s
Minimum time-step: 1e-7s
Maximum time-step: 1s
Turbulence Modeling: RNG `k-epsilon`
Geometry Dimensions
Mesh
Base-size: 0.8m
AMR: Sub grid scale (SGS) Temperature sub-grid criterion 0.05K(as an initial value)
we are simulating our problem with multiple values of SGS to see the effect of the temperature sub-grid citerion which is the second derivative of temperature within space (e.g `(partial^2T)/(partialx^2)`).
Second derivative of a function denotes to the curvature of the function, we are interested in where the temperature varies within the Domain, we know that when supersonic shock occurs velocity and temperature due to shock wave tend to vary dramatically, so we refine mesh to these specific places in the octree structure mesh to capture expansion fan (shock waves) more properly.
Max embedding level means that minimum value of the mesh element will be `(base)/2^((Embed- Level)) = 0.8/2^2= 0.2m`
Mesh Generation
POST PROCESSING
SGS 0.05K
SGS 0.03K
SGS 0.01K
It\'s clear that if we decrease SGS it means that mesh will be refined when the solver will detect even a little value of the second derivative of temperature, meaning we will capture shock waves more properly (it\'s noticable, when we decrease SGS lines of the shock waves become sharper).
If we increase the value of the SGS we will see that shock waves will become blur.
SGS 0.1K
SGS 0.2K
Velocity plots
Temperature Plots
Total cell count
In the image above it\'s clear that like we setted after 5000 iterations mesh is refined, in order to capture shock lines more properly.
100m/s initial condition case
Instead of setting 100m/s as an initial condition, we see that it did not effected the solution. At the very beggining of the simulation 100m/s was swept out of the domain and the simulation became the same when as an initial condition 680m/s was set. (just scaling colors are different but if you look at the values you will see that this animation shows the same as the animation of the simulaiton in which initial condition is 680m/s.)
Analytical and numerical
As we predicted with the analytical calculations Mach number `M_2 =2.19997` and temperature `T_2 = 161.68K` are pretty close to what we got from numerical analysis.
Conclusion
we successfully simulated prandtl-meyer expansion wave problem and learnt that SGS value is pretty important while predicting the position and the shape of the shock waves. As we were setting lower values of SGS we were getting sharper and hopefully more appropriate shapes of the shock wave lines.