Week 5: Prandtl Meyer Shock problem

1. AIM:

• Simulate supersonic flow over a sharp corner to study Prandtl Meyer Expansion fan.
• Explain the phenomena of expansion fan.
• Expand upon the understanding of supersonic flows and shockwaves.
• Discuss boundary conditions as a special case for flows involving shocks.
• Propound on the AMR algorithm and SGS value for Converge adaptive meshing.
• Post process in Paraview.

2. GEOMETRY, BOUNDARY CONDITIONS AND INITIALIZATION:

• Geometry is a sharp corner of about 5° angle at the corner. 
• Inlet Boundary conditions: Dirichlet velocity of 680 m/s with temperature of inlet air of 286K. This makes the mach number of inlet flow as 2.0. Mach number is a ratio of velocity of an object or flow vs speed of sound in that medium. At any flow abbove mach 0.3 the compressibility effects of the air start to dominate and hence a density based solver has to be used to accurately account for it.
• Outlet Boundary conditions: At the outlet Neumann boundary condition was applied for both velocity and pressure.
• Wall: The top and bottom walls were given a slip condition.
• Timestepping and Cycles: Simulation was run for 25000 cycles and in some cases upto 50000 cycles with time minimum time step of 1E-09 and max of 1 second. The convection, diffusion and mach CFLs were kept at standard values of 1, 2 and 50 respectively.  The solution was continously monitored and CFL values were appropriately changed to accelerate the convergence rate.

What is a shock? And why does the expansion fan phenomena occur anyway?

Shocks are fundamentally a discontinuity or an abrupt change in flow field properties in a domain. Shock flows are most commonly known and encountered around the supersonic aircrafts, almost everyone has seen and heard of a distinct loud sound of an airplane flying over the speed of sound and has also seen a condensation cloud forming around the aircraft also called as vapour cone or shock cone. The cloud formation is because every high pressure shock wave front there is a normal or expansion zone, which is where air cools down and the moisure condenses to form that signature cone shape.

What happens when a body is moving faster than the speed of sound itself in the domain, the information which is communicated through domain at speed of sound is actually not fast enough. This causes the sound waves or disturbances in the medium to 'pile up' and at the region where the waves 'add up' there is a huge pressure discontinuity which moves into the domain as a shock wave. For a supersonic bogey, it can be imagined as a cone emanating beside the moving body as shown in figure above. This is also pretty common for large firearm bullets that travel over the speed of sound, here is a schlieren imaging of bullet going supersonic. You can see multiple oblique shock waves wherever the flow goes supersonic as it travels around the bullet disturbed by it.

The shock waves could be of multiple types, specifically when considering aerodynamic applications shockwaves are categorised into 3 main types:

1) Oblique shock waves

2) Normal shock waves

3) Bow shocks or Curved shock fronts

This can be made clear by a small image below:

The expansion fans are slightly different. Shocks are generated when there are multiple overlapping mach waves, so technically a supersonic body can move through domain without shock provided the mach waves are spaced appart sufficiently. The expansion fan however is generated anytime the fluid encounters the expansion especially a sharp corner. There are infinite mach waves in this region which help the fluid turn. 

The flow actually accelerates and velocity increases, the static pressure, temperature and density decreases. It might be counter intuitive since under sub sonic flow when fluid encounters expansion, it actually slows down. This velocity increase is because the flow was high compressed before the sharp expansion corner, the density then is much higher. As the flow encounters expansion area, the flow itself is expanded; that is its density drops and pressure drops; the velocity increases. This process is isentropic so the stagnation properties would remain the same. Some formulas to calculate the basic quantities of expansion fan around the sharp corner.

where `nu` is actually Prandtl-meyer function with M1 and M2, mach numbers before and after the expansion fan. It is defined as

The rest of the quantities can then be calculated as 

What are boundary conditions and different types of it?

Any CFD code first discretizes the domain in some finite control volumes, then it needs some initialization (read: initial guess value) and the boundary conditions (read: conditions at the boundaries of domains). It formulates the complete set of equations in a matrix form and uses different iterative techniques to solve them.

So with proper initial guess value, the boundary condition is also very important. This is not just to achieve convergence but also to reproduce physically relevant results. 

There are mostly following boundary conditions used in conventional CFD codes:

1) Dirichlet boundary conditions:

This type of boundary condition basically specifies the value of transport variables and scalars which has to be take by solution directly at the boundary. For instance: providing a temperature value directly at the face so that the convection of the heat can be studied at constant face temperature.

For above ODE with 'a' and 'b' boundaries, Dirichlet condition is specified as:

2) Neumann Boundary conditions:

Unlike Dirichlet boundary conditions where the value is directly specified on the face, the Neumann boundary specifies a value of normal gradient of particular scalar or variable at the face (in most cases zero normal gradient). For instance: providing zero normal gradient at the wall without providing the temperature directly to allow for the code to calculate and project the values on the wall based on domain informations. 

This is why we applied Neumann conditions at the outlet of our domain here, as you can see in images below the value of variables at the outlet face are not constant. Infact they are varying from bottom to top with expansion fan spreading out upto the outlet face. Below image demonstrates the variation with a Neumann boundary at outlet and a Dirichlet pressure boundary of 500000Pa. One can note that although the guess for outlet boundary was almost in ball park, the results lack verisimilitude.

Neumann Outlet condition                              With dirichlet outlet                                          AMR on outlet due to gradient

For the same ODE mentioned above with same boundaries where Dirichlet condition was specified, if we specified Neumann conditions it would look something like:

3) Mixed boundary conditions:

Although very rare but there can be ocassions where boundary conditions are the face are not exactly known but not totally unknown either. So it wouldn't be ideal to use any of Dirichlet or Neumann conditions alone, in such cases we can apply one condition to some portion of the boundary and another to some grouped portion of the same boundary depending on which one is suitable where. Such conditions are mixed.

4) Robin boundary conditions:

A lot of times the robin boundary conditons are confused with mixed boundary conditions. In mixed boundary conditions, both dirichlet and neumann condition is applied to portions of faces which are disparate to each other. Robin boundary condition can be termed as a true hybrid of dirichlet and neumann boundary conditions.

Mathematically in simplest form, Robin conditons can look like:

here first term is Dirichlet term and second term in Neumann term clearly. 'a' and 'b' are the coefficients and could take a form of function essentially blending between Dirichlet-Neumann conditions.

5) Cauchy boundary conditions:

This is similiar to Robin in terms of both Dirichlet and Neumann conditions are applied to the boundary. But unlike a blend of linear function 'g' as mentioned before, in Cauchy boundary conditions we apply both Neumann and Dirichlet conditions throughout the face, this sets it apart from Mixed conditions as well. 

REGION INITIALIZATION:

Agood initial guess for a region will help the solver to converge on the solution much faster than a very bad initial guess, the region thus was initialized with 50000 Pa pressure and 300K of air in the domain.

We made two comparision, one with initialization of 100m/s and another with 680m/s same as inlet. As seen in graph below both. The graph on right represents an initialization of 680m/s, the jump in velocity is small. The graph on the left however, is with a poor guess value of 100m/s, the jump is massive in few initial iterations itself. However it must be noted that both soon converged to near abour 710m/s, but a good initial guess value would itself be around 710m/s. This would have put us in a position much closer to the solution and converge much faster. 

In extreme cases, there might be a possibility that solution might not even converge or product un-realistic results due to misfit initial conditions of domain.

Mesh was configured for base size of 0.8 in all three directions. Adaptive mesh refinement was used with different sub grid scales to study the effect of new grids and understand the Converge cartesian automatic meshing algorithm which will be explained later. All in all without AMR we have about 4500 cells in the domain.

3. RESULTS

As we used slip boundary condition on the wall, as expected there is no boundary layer on the walls. The results are what we expected to be, there is a fan of infinite mach waves emanating from the sharp corner helping the flow turn and expands as it moves through the domain. To demonstrate that we have a screenshot of highly refined density and velocity contour. As expected the density drops further down the flow direction and the velocity thus increases. The compressibility effects are highly dominant hence in the case of supersonic flow or infact any flow above M>0.3. 

As expected the pressure and temperature contour remains almost constant throught the expansion fan but it decreases drastically across the domain. Since the process is isentropic ideally, the Total pressure actually stays constant as expected.

Adaptive Mesh Refinement:

As we have already seen, Converge requires minimal input from the user's end for its own mesh generation. It automatically generates the best possible cartesian mesh for its domain. Ideally we would want as coarse mesh as possible with refinements at the places only and only where it is absolutely needed. The AMR algortihm of converge does just that, it adds embedding determined by algorithm to the portion of domain where the flow field variables are least resolved or so to speak where sub-grid field is largest. In simple words, it will check for the curvature gradients of the variable in the spatial domain and compares it to user defined tolerance. 

In converge, this subgrid field is defined as a difference of actual field and resolved field,

The sub-grid for any scalar can be expressed as an infinite series (Bedford and Yeo
(1993) and Pomraning (2000), as is given by

Since it is not possible to evaluate the entire series, only the first term (the second-order
term) in the series is used to approximate the scale of the sub-grid :

A cell is embedded if the absolute value of the sub-grid field is above a user-specified value.
Conversely, a cell is released (i.e., the embedding is removed) if the absolute value of the
sub-grid is below 1/5th of the user-specified value.

Back to Prandtl Meyer Simulation:

The above simulation was run with 3 different sub-grid scale values of 0.1, 0.05 and 0.01. The embedding level was set to 2, which means that a single control volme of 0.8m cell size we specified will be broken down upto 16 elements of 0.2m cell size. This embedding was given to the temperature as a scalar field for the adaptive algorithm. So basically the `phi`mentioned above becomes Temperature in this case.

Here are the three meshes and their refined grids at 25000 cycles.

Following shows clear difference between the velocity contours of the various SGS. From left to right we have more and more refinement of SGS from 0.1 to 0.01. As one can see with finer resolution a very sharp contour is produced as results are much more refined compared to very diffused coarse grid solution of the highest SGS criteria on the left.

One can see similiar case of refinement in Pressure contours from coarse grids to highly refined grids.

Temperature contour

Density contour

LINE PLOTS & DISCUSSIONS:

SubGrid-Scale criterion was varied from 0.1 to 0.01 and as one can see there is a variation of cell count throughout the simulation as the algorithm proceeds embedding multiple cells based on the SGS criteria. The smaller the criteria was, more number of cells were formed in the domain - hence more resolution but slower simulation as the timestep would decrease to meet various CFL constraints. This leads to the trading off of time-accuracy of simulation by the user. 

Cell count for SGS of 0.1

Cell count for SGS of 0.05

Cell count for SGS of 0.01

For reference, following plots are enclosed for temperature, pressure and mass flow rate through the inlet and outlet boundaries. For this particular instance, we have taken SGS of 0.01 and about 50000 cycles

Pressure plot:

Mass flow rate plot:

Temperature plot:

.

Videos of how the solution psuedo-transits through the time shows minor oscillations before finally converging to an expansion fan, the solution was stopped at 50000 iterations however it has still not fully reached the steady required state. But even under such conditions, the fan and related quantities can be seen changes as expected.

What if we run same simulation but with Sub-sonic velocity?

As expected, if the flow is subsonic, there is infact a decrease in velocity at certain areas of expansion. There is almost no change in density as expected throughout the domain.