Prandtl Meyer Shock Problem
AIM:
- To set up the case of Prandtl Meyer shock problem in Converge studio.
- Process the inputs from pre-processing done in converge studio using Cygwin terminal.
- Post process the results in Paraview.
- To do a literature search on what boundary conditions are typically used for shock flow problems.
- To look at the effect of change in SGS temperature difference value on cell count and shock location.
INTRODUCTION AND THEORY:
- Physics behind shock waves:
When the speed of source equals the speed of sound in that medium, the wave fronts cannot escape the source. The resulting pile of waves forms a large amplitude of “sound barrier” that makes sustained flight at this speed difficult and risky.
The below picture shows the motion of a source in a medium whose velocity is equal to the velocity of the sound in that medium. The wavefronts are forming a barrier of sound as the source moves in the medium and that barrier is perpendicular to the direction of motion of the source :
What is a shock wave ?
When the speed of a source exceeds the speed of sound, the wave fronts lag behind the source in a cone-shaped region with the source at its vertex. The edge of the cone forms a supersonic (means the speed of the moving body is more than the speed of sound in that medium) wave front with an unusually large amplitude called a “shock wave”. When a shock wave reaches an observer a “sonic boom” is look heard.
The image below shows the formation of shock waves and we can also observe the lagging of the wavefronts behind its source :
Unlike ordinary sound waves, the speed of a shock wave varies with its amplitude. The speed of a shock wave is always greater than the speed of sound in the fluid and decreases as the amplitude of the wave decreases. When the shock wave speed equals the normal speed, the shock wave dies and is reduced to an ordinary sound wave.
Mach number:
The ratio of the speed of a moving object (v) to the speed of sound (c) in a fluid is known as the Mach number (Ma) in honor of ERNST MACH (1838–1916), the Moravian physicist, psychologist, and philosopher who studied sound and ballistics.
Ma = v/c
The Mach number is a dimensionless measure of speed common in aerodynamics. Mach 0.5 is half the speed of sound, Mach 2 is twice the speed of sound, and so on. Speeds less than the speed of sound have a Mach number between zero and one and are described as subsonic. Those greater than the speed of sound have Mach numbers greater than one are a described as supersonic. Speeds approximately equal to the speed of sound have Mach numbers approximately equal to one and are described as transonic.
The shock wave from a supersonic object is a cone composed of overlapping spherical wavefronts. As any one of these wavefronts forms, it propagates radially outward at speed c and acquires a radius ct. At the same time the source, traveling at speed v moves forward vt. These two displacements form the leg and hypotenuse, respectively, of a right triangle and can be used to determine the Mach angle u at the vertex of the shock cone.
Sinu = c/v
When an object travels slower than sound, the ratio in this equation is greater than one, and the equation does not have a real solution. This makes absolute sense as there is no shock wave to speak of at subsonic speeds. Traveling at the speed of sound makes the ratio equal one and results in a Mach angle of ninety degrees. At transonic speeds the shock wave is a wall of high pressure moving with the object, perpendicular to its velocity. Above the speed of sound, the ratio is less than one and the Mach angle is less than ninety degrees. The faster the object moves, the narrower the cone of high pressure behind it becomes. Measuring the vertex angle is thus a simple way to determine the speed of a supersonic object.
- Boundary conditions:
Boundary conditions are those mathematical equations which are associated with the condition or state of a particular differential equation or any function whose derivative (or the function itself) is involved in the differential equation at the extremes of the domain.
Such boundary conditions are predominantly of two types although there are other boundary conditions which can be used. They are: 1) Dirichlet boundary conditions and 2) Neumann boundary conditions.
Dirichlet boundary conditions are those boundary conditions which specify the state of the function involved in the differential equation, which means that the value of the function is known at all parts of the boundary. The main motive behind the use of boundary condition is to help in solving the differential equation. Consider an example of a channel flow in which the fluid flowing is air. In this case, let us say that the flow of the fluid is due to the difference in the pressure at the inlet and the outlet of the channel. So, while defining the boundary conditions we use can use the value of pressure at inlet and outlet conditions as the pressure is constant at inlet and outlet conditions respectively.
But, if we consider a situation where the air which is entering the channel in the above example is being heated at a constant rate and the air going out is being heated at a constant rate then, we cannot use the Dirichlet boundary condition. Here, we have to use the boundary condition which goes something like this: The pressure gradient with respect to time is equal to some constant at the inlet and the pressure gradient with respect to time at the outlet of the channel is constant. These boundary conditions in which the value of the derivative of the dependent variable function is known is called as Neumann Boundary Condition.
Why typically Neumann boundary conditions are used for shock wave problems?
For Shock flow problems, typically Neumann boundary conditions are used because, the shock waves, by nature are surrounded by the constant change in physical properties like pressure, temperature and velocity. Hence they need a condition at the boundary where the change in pressure, temperature or velocity can be defined at the boundary instead of defining the pressure or temperature or velocity at the inlet and outlet.
In the Prandtl Meyer shock expansion fan problem, at the outlet, we have considered Neumann boundary conditions because of the changes in the physical quantities at the outlet like pressure and velocity.
Prandtl Meyer expansion fan:
Expansion fans are an important phenomenon that occur when the supersonic waves encounter a convex corner or an expansion corner. The direction of the flow changes smoothly across the fan and the mach number increases. Each wave in the expansion fan changes the direction of the flow gradually.
The formation of expansion fan is shown in the figure below:
Where ,
M1 is the Mach number at inlet and M1>1 (Supersonic wave);
M2 is the Mach number at outlet with M2>M1;
B= Turn angle;
A1= arcsin(1/M1), corresponding to first Mach line;
A2= arcsin(1/M2), corresponding to last Mach line.
Setting up the case of Prandtl Meyer shock problem in Converge studio. Let it be case set up 1.
GEOMETRY:
The geometry of the Prandtl meyer expansion fan is shown below.
This geometry is drawn in any CAD software and is imported as an stl file. This type of file is a file that converts a file into a triangulated form of the CAD file format. This type of file is required to open in Converge studio.
The dimensions of the geometry are as follows:
l_x= 65m; l_y = 45.1525m; l_z= 8.151m.
After the geometry is imported, it is checked for any errors in the diagnosis dock.
It checks for the following problems in the geometry:
1) Intersections.
2) Nonmanifold edges.
3) Open edges.
4) Overlapping Tris.
5) Normal orientation.
6) Isolated tris.
The diagnosis dock is shown in the picture below.
We can click the find option marked in the above picture of diagnosis dock to check for the problems in the geometry.
As it can be seen in the above picture that there are no problems in the geometry.
Next step is to set up the case. In this process we input all the necessary parameters for before processing the solution for the simulation.
CASE SETUP:
i) Application type: In the application type we select the time-based type as in our case the simulation is not for an IC engine. In the version prior to converge studio 3.0, the application type can be selected as General flow.
ii) Materials: We check Gas simulation, reaction mechanism. We select Air as the pre-defined mixture
The values for the global transport parameters are taken as follows:
1) Prandtl Number: 0.9
2) Schmidt Number: 0.78
iii) Simulation parameters:
In simulation parameters we check run parameters, simulation time parameters and solver parameters.
Under run parameters set the simulation type to steady state and run the simulation on full hydro mode. Also in the miscellaneous section, check the steady state monitor. in the steady state monitor section add the quantity which you want to be monitored. In this case, Avg. Velocity is monitored at the inlet boundary.
Under the simulation time parameters the conditions are set as following:
1) start time = 0 cycles
2) End time = 25000 cycles
3) Initial time step = 1e-7 s
4) Minimum time step = 1e-7 s
5) Maximum time step = 1 s
6) Maximum convetion cfl limit = 1.0
In the solver parameters, the maximum convection cfl limit (final stage) is set to 0.5.
Navier-stokes solver scheme: PISO
iv) Boundary Conditions:
The boundary conditions are as follows:
Inlet: The boundary type INFLOW is selected.
Pressure boundary condition type is selected as “Zero normal gradient (NE)”
Velocity boundary condition type is set as “Specified value (DI)” and a value of 100 m/s is provided.
Temperature boundary condition type is set as “Specified value (DI)” and a value of 286.1 K is provided.
Species: Air with mass fractions of oxygen : 0.23 and Nitrogen: 0.77
Outlet: The boundary type OUTFLOW is selected.
Pressure boundary condition type is selected as “Zero normal gradient (NE)”
Velocity boundary condition type is selected as “Zero normal gradient (NE)”
Temperature boundary condition type is set as “Specified value (DI)” and a value of 300 K is provided.
Species: Air with mass fractions of oxygen : 0.23 and Nitrogen: 0.77
Front: The boundary type TWO_D is select
Back: The boundary type TWO_D is selected.
Slip: The boundary type WALL is selected.
Velocity boundary conditions:
Wall motion type: Stationary
Type of wall boundary: Slip boundary.
V) Regions and Initialization:
Initial velocity: 100 m/s.
vi) Mesh:
Base grid size: dx=dy=dz=0.8 m
In Grid control, enable Adaptive mesh refinement.
Under adaptive mesh refinement (AMR) set the following as mentioned below:
Maximum Cells: 200,000
Place the available regions to active regions.
Check the temperature tab.
Set the maximum embedding level as ‘2’.
Sub-grid criterion as 0.05 k.
Timing control type as ‘sequential’.
Now, the input files are exported to a folder and are processed in the software Cygwin terminal 64. Here the solver solves the steady state Navier stokes equations and the results thus obtained are post processed in the Software called Paraview.
The above case setup is run after changing the initial velocity at the inlet boundary from 100 m/s to 680 m/s.
Let us call this as case set up 2.
Again the case set up 1 is considered and the sub grid scale criterion is changed from 0.05k to the following values:
- 0.1k. Let us call this as Case set up 3.
- 0.01k. Let us call this as Case set up 4.
RESULTS AND DISCUSSION:
Case set up 1 (Initial velocity at the inlet = 100 m/s, Sub grid criterion = 0.05k):
The mesh size dx=dy=dz=0.8m is shown below:
The above mesh does not have any adaptive mesh refinement as it is a sub-sonic flow. The mesh size of the above mesh is dx=dy=dz=0.8m
The Velocity contour is shown below:
We cannot see any abrupt increase in the velocity magnitude at the convex corner as the speed of air is less than the speed of sound in the medium.
The Pressure contour is shown below:
There is no formation of expansion fan as the flow is subsonic.
The Temperature contour is shown below:
The temperature contour above, does not show a formation of expansion wave as the speed of air at the inlet is 100 m/s which is sub-sonic. The mach number is less than 1.
Plot of Average velocity at inlet and outlet boundaries vs the simulation cycles for mesh size dx=dy=dz= 0.8m is shown below:
Average velocity at the inlet is 100 m/s.
Average velocity at the outlet varies upto first 1000 cycles of simulation then attains a constant value of 100 m/s.
This behavior of the velocity at the outlet indicates that no expansion wave is formed as the flow is sub-sonic.
Plot of Total pressure at inlet and outlet boundaries vs the simulation cycles is shown below:
From the above plot of total pressure we can observe that the total pressure at inlet overlaps with that of outlet. This indicates that there is no formation of expansion fan and the flow is sub-sonic.
Plot of Average temperature at inlet and outlet boundaries vs the simulation cycles is shown below:
From the above plot it can be observed that the average temperature at the inlet is 286.1k and at outlet, initiallly it varies upto 1500 simulation cycles and then reaches a constant value of about 271k.
Plot of Mass flow rate at inlet and outlet boundaries vs the simulation cycles is shown below:
From the above plot we can observe that the mass flow rate reaches a value of 0 kg/s which shows that steady state is reached.
Animations for Temperature, Velocity and Pressure for case set up 1 are as follows:
Link to the temperature animation is:
Link to the velocity animation is:
Link to the pressure animation is:
The above three animations for case set up 1 in which the initial velocity at the inlet is 100m/s shows that there is no formation of expansion fan at the convex corner. This indicates that the flow is sub-sonic in nature.
Case set up 2 (Initial velocity at the inlet = 680 m/s, Sub grid criterion = 0.05k):
The mesh size with base grid dx=dy=dz=0.8m is shown below:
The above mesh size shows the base grid size of dx=dy=dz=0.8m.
We can observe that some of the cells are further divided into finer cells and some of them are further divided into finer cells. This uneven, dynamic distribution of cells which changes according to the detailness of a particular quantity under study and at the particular time of requirement during the simulation is called Adaptive mesh refinement.
In Adaptive mesh refinement, the finest cells can be calculated as:
Finest cells = Base Grid size/(2^(embed level))
For case setup 2 it would be finest cells= 0.8/(2^2) = 0.2 m.
The velocity contour is shown below:
In the above velocity contour we can observe that the velocity of the air at the inlet is 680m/s and it continues to be so, until the air reaches the convex corner. After that the supersonic air flow changes its direction gradually at the convex corner and as a result of this we see that there is sudden increase in velocity of the air and forms mach lines. The first mach line has the lowest velocity among all the mach lines in expansion fan and the last mach line has the highest velocity. Then afterwards, the velocity gradually increases until it reaches the outlet.
Pressure contour is shown below:
In the above pressure contour we can observe that the pressure of air at the inlet is about 880000 Pa. It maintains the same value until it reaches the convex corner. At the convex corner, as the air gradually changes the direction of flow, we can see abrupt decrease in pressure and this region forms a shape of a fan. After this the pressure of air contiues to decrease until it reaches the outlet.
Temperature contour is shown below:
In the above temperature contour plot we can observe that the temperature of air at the inlet is around 286.1k and then it continues to remain same (with a little gradient of ~0.2k) until it reaches the same until it reaches the convex corner. At the convex corner as it gradually changes the direction of flow, the temperature of air abruptly decreases and if we observe the shape of the region in which this change occurs it looks like a fan. This is the Prandtl meyer expansion fan. After the air comes out of this region, the temperature of the air again decreases abruptly.
Plot of Average velocity at inlet and outlet boundaries vs the simulation cycles for base grid size dx=dy=dz= 0.8m is shown below:
In the above plot we can observe the change in average velocities at the inlet and outlet boundaries from the start to the end of this steady state simulation. At inlet throughout the simulation we can observe that the average velocity is 680m/s.This happens because we have maintained the velocity at the inlet throught the simultion to be 680m/s. At outlet we can observe that at the initial cycles of simulation, there is a slight increase in the velocity and then it attains a constant value. This is because initially the error in calculation of the velocity of the solver is more and then it gradually makes the error minimum dynamically.
Plot of Total pressure at inlet and outlet boundaries vs the simulation cycles is shown below:
In the above plot, we can observe that the total pressure at the inlet and the outlet boundaries increases and just above around 5000 cycles, the slope of the above plot suddenly decreases. The above nature of change does not necessarily show the exact change in pressure at each point of time during the simulation as this is a steady state simulation where we focus on the end result and not the path which the solver has taken to reach it.
Plot of Static pressure at inlet and outlet boundaries vs the simulation cycles is shown below:
In the above plot, we can observe that the static pressure at the inlet and the outlet boundaries increases and just above around 5000 cycles, the slope of the above plot suddenly decreases. The above nature of change does not necessarily show the exact change in pressure at each point of time during the simulation as this is a steady state simulation where we focus on the end result and not the path which the solver has taken to reach it.
Plot of Average temperature at inlet and outlet boundaries vs the simulation cycles is shown below:
In the above plot, we can observe that initially the average temperature of air at inlet and outlet boundaries is 286.1k and 300k respectively. The at the inlet continues to remain constant as we have specified the temperature at the inlet to be 286.1k and we are not increasing or decreasing the temperature of air at the inlet throughout the simulation. Where as, on the other hand, due to the abrupt increase in the temperature at the region of expansion fan, the solver had to adapt to the sudden changes in the average temperature inside the geometry which inturn changes the average temperature at the outlet. Due to this we can observe the changes in the above plot in case of the air average temperature at the outlet.
Plot of Mass flow rate at inlet and outlet boundaries vs the simulation cycles is shown below:
In the above plot we can observe that the mass flow rate at the inlet and at the outlet decreases and increases at the same rate. The above plot does not necessarily show the actual nature of air mass flow at each point of time during the simulation but shows the path followed by the solver to reach convergence as this is a steady state simulation.
Animations for Temperature, Velocity and Pressure for case set up 2 are as follows:
Link to the temperature animation is:
Link to the velocity animation is:
Link to the Pressure animation is:
In the above animations we can observe that how the temperature pressure and velocity of the air changes throughout the simulation. We can clearly see the formation of expansion fan and how the solver refines the data where more details about the temperature, pressure and velocity are required. These animations does not show the actual nature of the fluid flow in real time but shows the path followed by the solver to achieve convergence as this is a steady state simulation.
The effect of change in sub grid scale temperature difference on the total cell count in a simulation:
If in between any two sub-grids there exists a temperature difference that exceeds the value of sub grid scale temperature difference, the Adaptive mesh refinnement applies to those sub-grid cells and the cells get refined into more finer cells and hence the cell count increases.
So, if we increase the value of sub grid scale temperature difference, the refinement criteria of cells goes up and more temperature difference between two sub-grids is required to refine grid cells and this decreases the amount of cells which are refined hence the amount of cells are comparatively less.
Total cell count of the simulation for this case and by changing the Sub gridscale temperature to 0.1k and then to 0.01k VS the simulation cycles is shown below:
This plot shows the effect of sub grid scale temperature difference on cell count.
As we can see in the plot below that the total cell count is higher in the cases where the sub grid scale temperature difference has a comparatively higher value.
Also the change in the number of cells in each case set up can be seen by observing the meshes in each case setup as shown below:
For the case setup having SGC as 0.1k:
As we can see from the above mesh that fewer number of cells are further divided into finer cells compared to the other two cases.
For the case setup having SGC as 0.05k:
In this case more number of cells are further divided into finer cells when compared to the first case and hence it has more number of cells for simulation.
For the case setup having SGC as 0.01k:
In this case, when compared to other two cases has the largest number of cells which are divided into finer cells hence it has the highest number of cells for simulation.
Thus we can conclude that as the sub-grid scale temperature difference increases the number of cell count of simulation decreases.
The effect of change in Sub grid scale temperature difference on the location of expansion fan:
Below shows the temperature contours of the prandtl meyer shock problem with initial velocity at the inlet equal to 680 m/s but at different sub grid scale temperature differences of 0.1k, 0.05k and 0.01k:
Animations for the temperature contour plots with different SGS parameters:
Link to the temperature animation for SGS parameter 0.1k is:
Link to the temperature animation for SGS parameter 0.05k is:
Link to the temperature animation for SGS parameter 0.01k is:
As we can see from the above three contour plots and the animations of the temperature for the three different SGS temperature difference values of 0.1k, 0.05k and 0.01k.
We can see that the location of the expansion fan, the shock, has been narrowed down as we decrease the SGS parameter. Its nature becomes more clear and well-defined as we decrease the SGS parameter.
Conclusion:
In the Prandtl Meyer shock problem, we have observed how the properties of air like pressure, temperature and velocity change suddenly at the convex corner.
We observed the nature of sub-sonic and supersonic flows near the convex corner.
For more accurate analysis of the expansion fan, the sub grid scale temperature difference must be reasonably low. But this increases the computation time of simulation.